MeshFunction_Level.h 13.3 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
/******************************************************************************
 *
 * Extension of AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
 * Authors: Simon Praetorius et al.
 *
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 *
 *
 * See also license.opensource.txt in the distribution.
 * 
 ******************************************************************************/


#ifndef EXTENSIONS_MESH_FUNCTION_LEVEL_H
#define EXTENSIONS_MESH_FUNCTION_LEVEL_H
21

22
namespace AMDiS { namespace extensions {
23
24


25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
/** \brief
 * Abstract class that can be passed to RefinementLevel* as indicator where
 * to refine the mesh up to which level. It is an AbstractFunction that 
 * overloads the operator() method to return a level or a meshsize depending
 * on the coords/data passed to the operator.
 * You can switch between meshsize and level with the methods hToLevel(double) and 
 * levelToH(int)
 **/
template<typename T, typename T2>
class MeshRefinementFunction : public AbstractFunction<T2, T>
{
public:

  MeshRefinementFunction(Mesh* mesh_) :
    AbstractFunction<T2, T>(0),
    mesh(mesh_), 
    globalSize(0)
  {
      h0 = getMacroMeshSize(mesh);
      reduction = 1.0 / std::sqrt(2.0); // if dim==2
  }

  int getGlobalSize() { return globalSize; }

  double meshSize() { return h0; }
  
  virtual void init(ElInfo* elInfo) {}

  virtual T2 operator()(const T &value) const { return globalSize; }
  virtual T2 getLevel() const { return globalSize; }
  
  virtual double indicator(const T &value) const { return 1.0; }

protected:

  int hToLevel(double h) {
      int level = static_cast<int>(std::floor(std::log(h / h0) / std::log(reduction)));
      return level;
  }

  double levelToH(int level) {
      double h = std::pow(reduction,level)*h0;
      return h;
  }

  double getMacroMeshSize(Mesh* mesh) {
      FixVec<WorldVector<double>, VERTEX> coords = mesh->getMacroElement(0)->getCoord();
      double h = 0.0;
73
74
      for (int i = 0; i < coords.getSize(); ++i)
          for (int j = i + 1; j < coords.getSize(); ++j)
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
              h = std::max(h, norm(coords[i]-coords[j]));
      return h;
  }

protected:

  Mesh* mesh;

  int globalSize;

  double h0;
  double reduction;
};


90
91
92
93
94
95
96
97
98
99
100
101
/** \brief
 * Base-Implementation of a MeshRefinementFunction for refinement of phase-field
 * interfaces. operator() must be overwritten in subclasses. Class provides
 * basic parameters for the refinement like refinement-level values on interface
 * or in the interesting domain and defines the refinement range for the phase-
 * field variable, i.e. interface is defined in [0.05, 0.95]
 **/
template<typename T>
class PhaseFieldRefinementBase : public MeshRefinementFunction<T, int>
{
public:

102
  PhaseFieldRefinementBase(Mesh *mesh, std::string name = "mesh->refinement") : 
103
104
105
106
107
108
109
110
111
    MeshRefinementFunction<T, int>(mesh),
    lInner(10),
    lOuter(10),
    lInterface(14),
    minPhase(0.05),
    maxPhase(0.95),
    minOuterPhase(0.001),
    maxOuterPhase(0.999)
  {
112
113
114
      Parameters::get(name + "->level in inner domain",lInner);
      Parameters::get(name + "->level in outer domain",lOuter);
      Parameters::get(name + "->level on interface",lInterface);
115
116
117
118
119
120

      lInterface-= mesh->getMacroElementLevel();
      lInner-= mesh->getMacroElementLevel();
      lOuter-= mesh->getMacroElementLevel();

      int local_globalSize = 10;
121
      Parameters::get(name + "->initial level", local_globalSize);
122
123
      MeshRefinementFunction<T, int>::globalSize = local_globalSize;

124
125
126
127
      Parameters::get(name + "->min interface value",minPhase);
      Parameters::get(name + "->max interface value",maxPhase);
      Parameters::get(name + "->min outer interface value",minOuterPhase);
      Parameters::get(name + "->max outer interface value",maxOuterPhase);
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
  }

  double meshSize() { return levelToH(MeshRefinementFunction<T, int>::globalSize); }

protected:

  int lInner;
  int lOuter;
  int lInterface;
  
  double minPhase;
  double maxPhase;
  double minOuterPhase;
  double maxOuterPhase;
};


/** \brief
 * Implementation of PhaseFieldRefinementBase. Provides the method operator()
 * that defines the regions for interface refinement and domain refinement, etc.
 **/
class PhaseFieldRefinement : public PhaseFieldRefinementBase<double>
{
public:

153
154
155
    PhaseFieldRefinement(Mesh* mesh_, std::string name = "mesh->refinement") 
      : PhaseFieldRefinementBase<double>(mesh_, name) 
    { };
156
157
158
159
160
161
    
    int operator()(const double& phase) const {
      int result= lOuter;
      if (minPhase < phase && phase < maxPhase)
        result= lInterface;         // auf dem Interface
      else if (phase > minOuterPhase && phase <= minPhase)
162
        result= std::max(lOuter, static_cast<int>(std::floor((lOuter+lInterface)/2.0)));
163
      else if (phase < maxOuterPhase && phase >= maxPhase)
164
        result= std::max(lInner, static_cast<int>(std::floor((lInner+lInterface)/2.0)));
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
      else if (phase > (minPhase+maxPhase)/2.0)
        result= lInner;
      return result;
    }
};


/** \brief
 * Implementation of PhaseFieldRefinementBase. Provides the method operator()
 * that defines the regions for interface refinement and domain refinement, etc.
 * Additionally to the phase-field value a list of refinement points can be given
 * where the mesh should be refined up to the a given level, e.g. in corners of 
 * the macro-mesh.
 **/
class PhaseFieldCoordsRefinement : public PhaseFieldRefinementBase< std::pair<WorldVector<double>, double> >
{
public:
182
183
184
185
186
    PhaseFieldCoordsRefinement(Mesh* mesh_, std::vector<WorldVector<double> > points_, double radius_, std::string name = "mesh->refinement") 
      : PhaseFieldRefinementBase< std::pair<WorldVector<double>, double> >(mesh_, name),
	lPoints(14),
	points(points_),
	radius(radius_)
187
    {
188
        Parameters::get(name + "->level on points", lPoints);
189
190
191
192
193
194
195
196
197
        lPoints-= mesh->getMacroElementLevel();
    }
    
    int operator()(const std::pair<WorldVector<double>, double>& data) const {
      double phase = data.second;
      int result= lOuter;
      if (minPhase < phase && phase < maxPhase)
        result = lInterface;         // auf dem Interface
      else if ((phase > minOuterPhase && phase <= minPhase) || (phase < maxOuterPhase && phase >= maxPhase))
198
        result = std::max(lInner, static_cast<int>(std::floor((lOuter+lInterface)/2.0)));
199
200
201
202
203
204
205
206
207
      else if (phase > 0.5)
        result = lInner;

      WorldVector<double> x = data.first;
      double minDist = 1.e15;
      for (unsigned i = 0; i < points.size(); ++i) {
        minDist = std::min(minDist, norm(points[i]-x));
      }
      double lambda = std::max(0.0, 1.0 - minDist/radius);
208
      result = std::max(result, static_cast<int>(std::floor(lOuter+lambda*(lPoints-lOuter))));
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230

      return result;
    }
    
private:

    int lPoints;

    std::vector<WorldVector<double> > points;

    double radius;
};


/** \brief
 * Implementation of MeshRefinementFunction. Provides the method operator()
 * that defines a list of refinement points where the mesh should be refined 
 * up to the a given level, e.g. in corners of the macro-mesh.
 **/
class CoordsRefinement : public MeshRefinementFunction< WorldVector<double>, int >
{
public:
231
232
233
234
    CoordsRefinement(Mesh* mesh_, std::vector<WorldVector<double> > points_, double radius_, std::string name = "mesh->refinement") 
      : MeshRefinementFunction< WorldVector<double>, int >(mesh_),
	points(points_),
	radius(radius_)
235
236
    {
        lInner= 15; lPoints = 20;
237
238
        Parameters::get(name + "->level on points",lPoints);
        Parameters::get(name + "->level in inner domain",lInner);
239
240
241
242
243
244
245
246
247
248
249
250
251
        lPoints-= mesh->getMacroElementLevel();
        lInner-= mesh->getMacroElementLevel();

        globalSize= lInner;
    }
    
    int operator()(const WorldVector<double>& x) const {
      double minDist = 1.e15;
      for (unsigned i = 0; i < points.size(); ++i) {
        minDist = std::min(minDist, norm(points[i]-x));
      }
      int result = lInner;
      double lambda = std::max(0.0, 1.0 - minDist/radius);
252
      result = std::max(result, static_cast<int>(std::floor(lInner+lambda*(lPoints-lInner))));
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272

      return result;
    }
    double meshSize() { return levelToH(lInner); }
    
private:
    int lInner, lPoints;
    std::vector<WorldVector<double> > points;
    double radius;
};


/** \brief
 * Implementation of PhaseFieldRefinementBase. Provides the method operator()
 * that defines the level of refinement depending on a list of phase-fields, e.g.
 * a diffuse-domain representation and a Cahn-Hilliard field.
 **/
class PhaseFieldChRefinement : public PhaseFieldRefinementBase< std::vector<double> >
{
public:
273
274
275
  PhaseFieldChRefinement(Mesh* mesh_, std::string name = "mesh->refinement") 
    : PhaseFieldRefinementBase< std::vector<double> >(mesh_, name),
      lInterfaceDomain(14)
276
  {
277
    Parameters::get(name + "->level on domain interface",lInterfaceDomain);
278
279
    lInterfaceDomain-= mesh->getMacroElementLevel();
    
280
    Parameters::get(name + "->initial level", globalSize);
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
  }
  
  int operator()(const std::vector<double>& phases) const 
  {
    int result= lOuter;
    for (unsigned i = 0; i < phases.size(); ++i) {
      int localResult= lOuter;
      if (minPhase < phases[i] && phases[i] < maxPhase
        && (i < 1 || phases[i-1] > 0.5))
        localResult= (i==0 ? lInterfaceDomain : lInterface);         // auf dem Interface
      else if (phases[i] > 0.5
        && (i < 1 || phases[i-1] > 0.5))
        localResult= lInner;            // im Innern des Gebietes
      if (((phases[i] > minOuterPhase && phases[i] <= minPhase) || (phases[i] < maxOuterPhase && phases[i] >= maxPhase))
        && (i < 1 || phases[i-1] > 0.5))
296
        localResult= static_cast<int>(std::floor((lInner + lInterface) / 2.0));
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
      result= std::max(result, localResult);
    }
    return result;
  }
    
private:

  int lInterfaceDomain;
};


/** \brief
 * Base-Implementation of a MeshRefinementFunction for refinement of phase-field
 * interfaces. operator() must be overwritten in subclasses. Class provides
 * basic parameters for the refinement like refinement-level values on interface
 * or in the interesting domain and defines the refinement range for the phase-
 * field variable, i.e. interface is defined in [0.05, 0.95]
 **/
template<typename T>
class SignedDistRefinementBase : public MeshRefinementFunction<T, int>
{
public:

320
  SignedDistRefinementBase(Mesh *mesh, std::string name = "mesh->refinement") :
321
    MeshRefinementFunction<T, int>(mesh),
322
    name(name),
323
324
325
326
327
328
329
    lInner(10),
    lOuter(10),
    lInterface(14),
    interfaceWidth(0.2),
    fadeOutWidth(0.0),
    signInInnerDomain(-1.0)
  {
330
331
332
      Parameters::get(name + "->level in inner domain",lInner);
      Parameters::get(name + "->level in outer domain",lOuter);
      Parameters::get(name + "->level on interface",lInterface);
333
334
335
336
337
338

      lInterface-= mesh->getMacroElementLevel();
      lInner-= mesh->getMacroElementLevel();
      lOuter-= mesh->getMacroElementLevel();

      int local_globalSize = 10;
339
      Parameters::get(name + "->initial level", local_globalSize);
340
341
      MeshRefinementFunction<T, int>::globalSize = local_globalSize;

342
343
344
      Parameters::get(name + "->interface width",interfaceWidth);
      Parameters::get(name + "->fade out width",fadeOutWidth);
      Parameters::get(name + "->sign in inner domain",signInInnerDomain);
345
346
347
348
349
  }

  double meshSize() { return levelToH(MeshRefinementFunction<T, int>::globalSize); }

protected:
350
  std::string name;
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369

  int lInner;
  int lOuter;
  int lInterface;

  double interfaceWidth;
  double fadeOutWidth;
  double signInInnerDomain;
};


/** \brief
 * Implementation of PhaseFieldRefinementBase. Provides the method operator()
 * that defines the regions for interface refinement and domain refinement, etc.
 **/
class SignedDistRefinement : public SignedDistRefinementBase<double>
{
public:

370
    SignedDistRefinement(Mesh* mesh_, std::string name = "mesh->refinement") : SignedDistRefinementBase<double>(mesh_, name) {};
371
372
373

    int operator()(const double& dist) const {
      int result= lOuter;
374
      if (std::abs(dist) < interfaceWidth)
375
        result= lInterface;         // auf dem Interface
376
377
      else if (std::abs(dist) < fadeOutWidth+interfaceWidth && signInInnerDomain*dist < 0.0) {
	double lambda = std::abs(dist)/fadeOutWidth - interfaceWidth/fadeOutWidth;
378
        result= static_cast<int>(std::floor(lambda*lInner+(1.0-lambda)*lInterface));
379
380
      } else if (std::abs(dist) < fadeOutWidth+interfaceWidth && signInInnerDomain*dist > 0.0) {
	double lambda = std::abs(dist)/fadeOutWidth - interfaceWidth/fadeOutWidth;
381
        result= static_cast<int>(std::floor(lambda*lOuter+(1.0-lambda)*lInterface));
382
383
384
385
386
387
      } else if (signInInnerDomain*dist < 0.0)
        result= lInner;
      return result;
    }
};

388

389
390
391
class ESIndicator : public MeshRefinementFunction< std::vector<double>, int >
{
public:
392
393
  ESIndicator(Mesh* mesh_, std::string name = "mesh->refinement") 
    : MeshRefinementFunction< std::vector<double>, int >(mesh_)
394
  {
395
    Parameters::get(name + "->initial level", globalSize);
396
397
398
399
400
  }
  
  int operator()(const std::vector<double>& means) const 
  {
    double tol = means[0];
401
    return (std::abs(means[1]-means[2]) > tol ? 0 : -1);
402
403
404
405
  }
    
  double indicator(const std::vector<double>& means) const 
  {
406
    return std::abs(means[1]-means[2]);
407
408
  }
};
409

410
411
} }

412
#endif // EXTENSIONS_MESH_FUNCTION_LEVEL_H