Peano 4
Loading...
Searching...
No Matches
refinement.h
Go to the documentation of this file.
1#ifndef EXASEIS_REFINEMENT_HEADER
2#define EXASEIS_REFINEMENT_HEADER
3
4// #include "kernels/GaussLegendreBasis.h"
6
7namespace Refinement {
8 template <class Shortcuts>
10 public:
11 virtual bool eval(
12 const double* luh,
15 double t,
16 const int level,
18 ) = 0;
19 };
20 template <class Shortcuts>
21 class trackVelocity: public RefinementCriterion<Shortcuts> {
22 public:
23 trackVelocity(int a_basisSize, double a_refine_threshold, double a_coarsen_threshold):
24 basisSize(a_basisSize),
25 refine_threshold_2(a_refine_threshold * a_refine_threshold),
26 coarsen_threshold_2(a_coarsen_threshold * a_coarsen_threshold){};
27
28 bool eval(
29 const double* luh,
32 double t,
33 const int level,
35 ) override {
36 Shortcuts s;
37 kernels::idx4 idx(basisSize, basisSize, basisSize, s.SizeVariables + s.SizeParameters);
38 double maxVelocity_2 = 0;
39 double minVelocity_2 = std::numeric_limits<double>::max();
40 for (int i = 0; i < basisSize; i++) {
41 for (int j = 0; j < basisSize; j++) {
42 for (int k = 0; k < basisSize; k++) {
43 double velocity_2 = luh[idx(i, j, k, s.v + 0)] * luh[idx(i, j, k, s.v + 0)]
44 + luh[idx(i, j, k, s.v + 1)] * luh[idx(i, j, k, s.v + 1)]
45 + luh[idx(i, j, k, s.v + 2)] * luh[idx(i, j, k, s.v + 2)];
46
47 maxVelocity_2 = std::max(maxVelocity_2, velocity_2);
48 minVelocity_2 = std::min(minVelocity_2, velocity_2);
49 }
50 }
51 }
52
54
55 if (refine_threshold_2 < maxVelocity_2) {
57 return false;
58 }
59
60 if (coarsen_threshold_2 > minVelocity_2) {
62 return false;
63 }
64
65 return false;
66 };
67
68 private:
72 };
73
74 template <class Shortcuts>
75 class StaticAMR: public RefinementCriterion<Shortcuts> {
76 public:
77 bool eval(
78 const double* luh,
81 double t,
82 const int level,
84 ) override {
85 // static AMR
86 if (!tarch::la::equals(t, 0.0)) {
88 return true;
89 }
90 return false;
91 }
92 };
93
94 template <class Shortcuts>
96 public:
97 RefineBetweenPositions(int level, double a_left[3], double a_right[3]):
98 max_level(level) {
99 rectangle_left[0] = a_left[0];
100 rectangle_left[1] = a_left[1];
101 rectangle_left[2] = a_left[2];
102
103 rectangle_right[0] = a_right[0];
104 rectangle_right[1] = a_right[1];
105 rectangle_right[2] = a_right[2];
106 }
107
108 bool eval(
109 const double* luh,
112 double t,
113 const int level,
115 ) override {
116 double left_vertex[3];
117 double right_vertex[3];
118 bool elementBetweenPoints = true;
119
120 for (int i = 0; i < Dimensions; i++) {
121 left_vertex[i] = center[i] - dx[i] * 0.5;
122 right_vertex[i] = center[i] + dx[i] * 0.5;
123 elementBetweenPoints &= ((left_vertex[i] < rectangle_right[i]) && (right_vertex[i] > rectangle_left[i]));
124 }
125
126 if (elementBetweenPoints && (max_level > level)) { // solver->getMaximumAdaptiveMeshLevel()-1
128 return true;
129 }
130 return false;
131 }
132
133 private:
134 double rectangle_left[3];
137 };
138
139 template <class Shortcuts>
141 public:
142 RefineDownToPositionCustomCoordinates(int a_basisSize, int level, double a_position[3], int a_top, double* a_FL):
143 basisSize(a_basisSize),
144 max_level(level),
145 top(a_top) {
146 position[0] = a_position[0];
147 position[1] = a_position[1];
148 position[2] = a_position[2];
149 FL = new double[a_basisSize];
150 std::copy_n(a_FL, a_basisSize, FL);
151 }
152
153 bool eval(
154 const double* luh,
157 double t,
158 const int level,
160 ) {
161 double left_vertex = std::numeric_limits<double>::min();
162
163 Shortcuts s;
164 kernels::idx4 idx(basisSize, basisSize, basisSize, s.SizeVariables + s.SizeParameters);
165 for (int j = 0; j < basisSize; j++) {
166 for (int k = 0; k < basisSize; k++) {
167 double vertex = 0;
168 for (int i = 0; i < basisSize; i++) {
169 if (top == 2) {
170 vertex += luh[idx(i, j, k, s.curve_grid + top)] * FL[i];
171 } else if (top == 1) {
172 vertex += luh[idx(j, i, k, s.curve_grid + top)] * FL[i];
173 } else if (top == 0) {
174 vertex += luh[idx(j, k, i, s.curve_grid + top)] * FL[i];
175 }
176 }
177 left_vertex = std::max(left_vertex, vertex);
178 }
179 }
180
181 bool elementAbovePointSource = (left_vertex < position[top]);
182 if (elementAbovePointSource) { // solver->getMaximumAdaptiveMeshLevel()-1
184 return true;
185 }
186 return false;
187 }
188
189 private:
190 double position[3];
191 double* FL;
194 int top;
195 };
196
197 template <class Shortcuts>
198 class RefineDownToPosition: public RefinementCriterion<Shortcuts> {
199 public:
200 RefineDownToPosition(int a_basisSize, int level, double a_position[3], int a_top):
201 basisSize(a_basisSize),
202 max_level(level),
203 top(a_top) {
204 position[0] = a_position[0];
205 position[1] = a_position[1];
206 position[2] = a_position[2];
207 }
208
209 bool eval(
210 const double* luh,
213 double t,
214 const int level,
216 ) {
217 double left_vertex[Dimensions];
218 for (int d = 0; d < Dimensions; d++) {
219 left_vertex[d] = center[d] - dx[d] * 0.5;
220 }
221
222 bool elementAbovePointSource = (left_vertex[top] < position[top]);
223 if (elementAbovePointSource) { // solver->getMaximumAdaptiveMeshLevel()-1
225 return true;
226 }
227 return false;
228 }
229
230 private:
231 double position[3];
234 int top;
235 };
236
237 template <class Shortcuts>
239 public:
240 RefinePositionCustomCoordinates(double a_position[3], int a_basisSize, double* a_FL, double* a_FR) {
241 position[0] = a_position[0];
242 position[1] = a_position[1];
243 position[2] = a_position[2];
244
245 basisSize = a_basisSize;
246 FR = new double[a_basisSize];
247 std::copy_n(a_FR, a_basisSize, FR);
248
249 FL = new double[a_basisSize];
250 std::copy_n(a_FL, a_basisSize, FL);
251 };
252
253 inline bool eval(
254 const double* luh,
257 double t,
258 const int level,
260 ) {
261 double left_vertex[3];
262 double right_vertex[3];
263 Shortcuts s;
264 bool positionInElement = true;
265 kernels::idx4 idx(basisSize, basisSize, basisSize, s.SizeVariables + s.SizeParameters);
266 for (int i = 0; i < basisSize; i++) {
267 left_vertex[0] += luh[idx(0, 0, i, s.curve_grid + 0)] * FL[i];
268 right_vertex[0] += luh[idx(basisSize - 1, basisSize - 1, i, s.curve_grid + 0)] * FR[i];
269
270 left_vertex[1] += luh[idx(0, i, 0, s.curve_grid + 1)] * FL[i];
271 right_vertex[1] += luh[idx(basisSize - 1, i, basisSize - 1, s.curve_grid + 1)] * FR[i];
272
273 left_vertex[2] += luh[idx(i, 0, 0, s.curve_grid + 2)] * FL[i];
274 right_vertex[2] += luh[idx(i, basisSize - 1, basisSize - 1, s.curve_grid + 2)] * FR[i];
275 }
276
277 for (int i = 0; i < Dimensions; i++) {
278 positionInElement &= ((left_vertex[i] <= position[i]) && (right_vertex[i] >= position[i]));
279 }
280
281 if (positionInElement) {
283 return true;
284 }
285 return false;
286 }
287
288 private:
289 double position[3];
290 double* FL;
291 double* FR;
293 };
294
295 template <class Shortcuts>
296 class RefinePosition: public RefinementCriterion<Shortcuts> {
297 public:
298 RefinePosition(double a_position[3]) {
299 position[0] = a_position[0];
300 position[1] = a_position[1];
301 position[2] = a_position[2];
302 };
303
304 inline bool eval(
305 const double* luh,
308 double t,
309 const int level,
311 ) {
312 double left_vertex[3];
313 double right_vertex[3];
314
315 bool positionInElement = true;
316
317 for (int i = 0; i < Dimensions; i++) {
318 left_vertex[i] = center[i] - dx[i] * 0.5;
319 right_vertex[i] = center[i] + dx[i] * 0.5;
320 }
321
322 for (int i = 0; i < Dimensions; i++) {
323 positionInElement &= ((left_vertex[i] <= position[i]) && (right_vertex[i] >= position[i]));
324 }
325
326 if (positionInElement) {
328 return true;
329 std::cout << "refine ps" << std::endl;
330 }
331 return false;
332 }
333
334 private:
335 double position[3];
336 };
337
338 template <class Shortcuts>
339 class CoarseBoundaryLayer: public RefinementCriterion<Shortcuts> {
340 private:
341 std::vector<double> threshold_p[3];
342 std::vector<double> threshold_m[2];
343
345 int max_depth; // Maximal tracked depth
347 int top;
348 int left;
349 int front;
350
351 public:
353 int a_max_depth,
354 int a_coarsest_level,
355 double a_coarsest_size,
356 int a_top,
357 const double domainSize[3],
358 const double domainOffset[3]
359 ):
360 max_depth(a_max_depth),
361 coarsest_level(a_coarsest_level),
362 coarsest_size(a_coarsest_size),
363 top(a_top) {
364
365 left = (top + 1) % 3;
366 front = (top + 2) % 3;
367
368 // We first compute properties of the mesh for the coarsest mesh size
369
370 // The coorindates for the coars layers around the fine domain
371 double coarse_layer[3];
372
373 // We aim for 1/3 of the domain in each direction
374 /* coarse_layer[left ] = 1.0/3.0;
375 coarse_layer[top ] = 1.0/3.0;
376 coarse_layer[front] = 1.0/3.0;*/
377
378 coarse_layer[left] = 12.5 / 25.0;
379 coarse_layer[top] = 12.5 / 25.0;
380 coarse_layer[front] = 12.5 / 25.0;
381
382 double coarse_elts[3];
383 // Number of elements for the estimated layer
384 coarse_elts[left] = domainSize[left] / coarsest_size;
385 coarse_elts[top] = domainSize[top] / coarsest_size;
386 coarse_elts[front] = domainSize[front] / coarsest_size;
387
388 // Adjust positions to fit with coarsest mesh size,
389 // we round down elements inside the fine domain
390 coarse_layer[left] = std::floor(coarse_elts[left] * coarse_layer[left]) / coarse_elts[left];
391 coarse_layer[top] = std::floor(coarse_elts[top] * coarse_layer[top]) / coarse_elts[top];
392 coarse_layer[front] = std::floor(coarse_elts[front] * coarse_layer[front]) / coarse_elts[front];
393
394 // for each possible depth (0 beeing the coarsest layer) we store the refinement threshold
395 // left
396 threshold_p[0].resize(max_depth);
397 threshold_m[0].resize(max_depth);
398
399 // front
400 threshold_p[1].resize(max_depth);
401 threshold_m[1].resize(max_depth);
402
403 // top
404 threshold_p[2].resize(max_depth);
405
406 for (int depth = 0; depth < max_depth; depth++) {
407 // the total number of elements in each dimenion for the current depth
408 double num_elts[3];
409 num_elts[left] = coarse_elts[left] * pow(3.0, depth);
410 num_elts[top] = coarse_elts[top] * pow(3.0, depth);
411 num_elts[front] = coarse_elts[front] * pow(3.0, depth);
412
413 // the number of elements on one side of each layer for the current depth
414 int num_layer_elements[3];
415 num_layer_elements[left] = coarse_layer[left] * num_elts[left];
416 num_layer_elements[top] = coarse_layer[top] * num_elts[top];
417 num_layer_elements[front] = coarse_layer[front] * num_elts[front];
418
419 // For each depth level we add a single element to the step wise layer
420 for (int d = 0; d < depth; d++) {
421 num_layer_elements[left] += std::pow(3.0, d);
422 num_layer_elements[top] += std::pow(3.0, d);
423 num_layer_elements[front] += std::pow(3.0, d);
424 }
425
426 // width of the coarse layer for the current depth
427 double threshold_width[3];
428 threshold_width[0] = domainSize[left] / num_elts[left] * num_layer_elements[left];
429 threshold_width[1] = domainSize[front] / num_elts[front] * num_layer_elements[front];
430 threshold_width[2] = domainSize[top] / num_elts[top] * num_layer_elements[top];
431
432 // boundaries of the coarse domain
433 // index for level = 2 is 0
434 threshold_m[0][depth] = domainOffset[left] + threshold_width[0];
435 threshold_m[1][depth] = domainOffset[front] + threshold_width[1];
436
437 threshold_p[0][depth] = domainOffset[left] + domainSize[left] - threshold_width[0];
438 threshold_p[1][depth] = domainOffset[front] + domainSize[front] - threshold_width[1];
439 threshold_p[2][depth] = domainOffset[top] + domainSize[top] - threshold_width[2];
440 }
441 }
442
443 bool eval(
444 const double* luh,
447 double t,
448 const int level,
450 ) {
451 // depth of the current element
452 int depth = level - coarsest_level;
453
454 if (depth >= max_depth) {
455 return false;
456 }
457
458 // check if cell is in the boundary layer
459 bool cellInBoundaryLayer = (center[left] < threshold_m[0][depth] || center[left] > threshold_p[0][depth])
460 || (/*center[1] < boundary_y_left ||*/ center[top] > threshold_p[2][depth])
461 || (center[front] < threshold_m[1][depth] || center[front] > threshold_p[1][depth]);
462
463 if (cellInBoundaryLayer) {
465 return true;
466 }
467
468 return false;
469 }
470 };
471
472 template <class Shortcuts>
474 public:
475 RefineCubeAroundPosition(double a_position[3], double radius) {
476 for (int i = 0; i < Dimensions; i++) {
477 left_vertice[i] = a_position[i] - radius;
478 right_vertice[i] = a_position[i] + radius;
479 }
480 }
481
482 bool eval(
483 const double* luh,
486 double t,
487 const int level,
489 ) {
490 bool inRadius = true;
491 for (int i = 0; i < Dimensions; i++) {
492 inRadius &= center[i] < right_vertice[i];
493 inRadius &= center[i] > left_vertice[i];
494 }
495
496 if (inRadius) {
498 return true;
499 }
500
501 return false;
502 }
503
504 private:
505 double left_vertice[3];
506 double right_vertice[3];
507 };
508
509 template <class Shortcuts>
510 class RefineFilterCube: public RefinementCriterion<Shortcuts> {
511 public:
512 RefineFilterCube(double a_left_vertice[3], double a_right_vertice[3], int a_level_inside, int a_level_outside):
513 level_inside(a_level_inside),
514 level_outside(a_level_outside) {
515 for (int i = 0; i < Dimensions; i++) {
516 left_vertice[i] = a_left_vertice[i];
517 right_vertice[i] = a_right_vertice[i];
518 }
519 }
520
521 bool eval(
522 const double* luh,
525 double t,
526 const int level,
528 ) {
529 bool left_vertice_in = true;
530 bool right_vertice_in = true;
531 for (int i = 0; i < Dimensions; i++) {
532 left_vertice_in &= center[i] - dx[i] < right_vertice[i];
533 left_vertice_in &= center[i] - dx[i] > left_vertice[i];
534
535 right_vertice_in &= center[i] + dx[i] < right_vertice[i];
536 right_vertice_in &= center[i] + dx[i] > left_vertice[i];
537 }
538
539 if (left_vertice_in || right_vertice_in) {
540 if (level == level_inside) {
543 }
544 }
545 } else {
546 if (level == level_outside) {
549 }
550 }
551 }
552 return true;
553 }
554
555 private:
556 double left_vertice[3];
557 double right_vertice[3];
560 };
561
562} // namespace Refinement
563
564#endif
And from this we can write down f$ nabla phi_i nabla phi_i dx but since we are constructing matrix let s investigate the f$ j
we integrate over each cell and then take the sum across each of the cells We also consider the terms that enter the f$ k
we integrate over each cell and then take the sum across each of the cells We also consider the terms that enter the f$ so we are integrating with f$ phi_k phi_k dx
std::vector< double > threshold_m[2]
Definition refinement.h:342
std::vector< double > threshold_p[3]
Definition refinement.h:341
bool eval(const double *luh, const tarch::la::Vector< Dimensions, double > &center, const tarch::la::Vector< Dimensions, double > &dx, double t, const int level, ::exahype2::RefinementCommand &refine)
Definition refinement.h:443
CoarseBoundaryLayer(int a_max_depth, int a_coarsest_level, double a_coarsest_size, int a_top, const double domainSize[3], const double domainOffset[3])
Definition refinement.h:352
bool eval(const double *luh, const tarch::la::Vector< Dimensions, double > &center, const tarch::la::Vector< Dimensions, double > &dx, double t, const int level, ::exahype2::RefinementCommand &refine) override
Definition refinement.h:108
RefineBetweenPositions(int level, double a_left[3], double a_right[3])
Definition refinement.h:97
bool eval(const double *luh, const tarch::la::Vector< Dimensions, double > &center, const tarch::la::Vector< Dimensions, double > &dx, double t, const int level, ::exahype2::RefinementCommand &refine)
Definition refinement.h:482
RefineCubeAroundPosition(double a_position[3], double radius)
Definition refinement.h:475
RefineDownToPositionCustomCoordinates(int a_basisSize, int level, double a_position[3], int a_top, double *a_FL)
Definition refinement.h:142
bool eval(const double *luh, const tarch::la::Vector< Dimensions, double > &center, const tarch::la::Vector< Dimensions, double > &dx, double t, const int level, ::exahype2::RefinementCommand &refine)
Definition refinement.h:153
bool eval(const double *luh, const tarch::la::Vector< Dimensions, double > &center, const tarch::la::Vector< Dimensions, double > &dx, double t, const int level, ::exahype2::RefinementCommand &refine)
Definition refinement.h:209
RefineDownToPosition(int a_basisSize, int level, double a_position[3], int a_top)
Definition refinement.h:200
bool eval(const double *luh, const tarch::la::Vector< Dimensions, double > &center, const tarch::la::Vector< Dimensions, double > &dx, double t, const int level, ::exahype2::RefinementCommand &refine)
Definition refinement.h:521
RefineFilterCube(double a_left_vertice[3], double a_right_vertice[3], int a_level_inside, int a_level_outside)
Definition refinement.h:512
RefinePositionCustomCoordinates(double a_position[3], int a_basisSize, double *a_FL, double *a_FR)
Definition refinement.h:240
bool eval(const double *luh, const tarch::la::Vector< Dimensions, double > &center, const tarch::la::Vector< Dimensions, double > &dx, double t, const int level, ::exahype2::RefinementCommand &refine)
Definition refinement.h:253
bool eval(const double *luh, const tarch::la::Vector< Dimensions, double > &center, const tarch::la::Vector< Dimensions, double > &dx, double t, const int level, ::exahype2::RefinementCommand &refine)
Definition refinement.h:304
RefinePosition(double a_position[3])
Definition refinement.h:298
virtual bool eval(const double *luh, const tarch::la::Vector< Dimensions, double > &center, const tarch::la::Vector< Dimensions, double > &dx, double t, const int level, ::exahype2::RefinementCommand &refine)=0
bool eval(const double *luh, const tarch::la::Vector< Dimensions, double > &center, const tarch::la::Vector< Dimensions, double > &dx, double t, const int level, ::exahype2::RefinementCommand &refine) override
Definition refinement.h:77
bool eval(const double *luh, const tarch::la::Vector< Dimensions, double > &center, const tarch::la::Vector< Dimensions, double > &dx, double t, const int level, ::exahype2::RefinementCommand &refine) override
Definition refinement.h:28
trackVelocity(int a_basisSize, double a_refine_threshold, double a_coarsen_threshold)
Definition refinement.h:23
examples::exahype2::elastic::VariableShortcuts s
Definition loh.cpp:10
bool equals(const Matrix< Rows, Cols, Scalar > &lhs, const Matrix< Rows, Cols, Scalar > &rhs, const Scalar &tolerance=NUMERICAL_ZERO_DIFFERENCE)
Compares to matrices on equality by means of a numerical accuracy.
Simple vector class.
Definition Vector.h:134