Peano 4
Loading...
Searching...
No Matches
DiffuseInterface.h
Go to the documentation of this file.
1#ifndef TEMPLATE_DIFFUSEINTERFACE_HEADER
2#define TEMPLATE_DIFFUSEINTERFACE_HEADER
3
4#include "../../../../ExaHyPE/kernels/GaussLegendreBasis.h"
5#include "../../../../ExaHyPE/kernels/KernelUtils.h"
6
7enum Position { above, at, below };
8
9template <class Shortcuts, int basisSize>
11public:
13 Topography* a_topo,
14 DomainInformation* a_domainInfo,
15 SolverInformation* a_solverInfo,
16 double a_interface_factor,
17 double a_max_distance_to_direct_surface
18 ) {
19 topography = a_topo;
20 info = a_domainInfo;
21 isDG = a_solverInfo->isDG();
22
23 max_distance_to_direct_surface = a_max_distance_to_direct_surface;
24
25 dx = info->getDx();
26
27 interface_factor = a_interface_factor;
29
30 // number of cells in one dimension
31 domain_cells = int(std::ceil(info->domainSize[0] / dx));
32
33 integration_nodes.resize(basisSize);
34 for (int node = 0; node < basisSize; node++) {
35 if (isDG) {
36 // todo replace with solver nodes
37 integration_nodes[node] = kernels::legendre::nodes[basisSize - 1][node];
38 } else {
39 integration_nodes[node] = (node + 1) / (basisSize + 1);
40 }
41 }
42 initLimits();
43 };
44
46
47 // position 0: above
48 // 1: at
49 // 2: below topography
50 int findRelativePosition(double& distance, double y_min, double y_max, int i, int j) {
52 distance = topography_limits_min[idx_cells(j, i)] - y_max;
53 if (distance > 0) {
54 return Position::above;
55 }
56
57 distance = y_min - topography_limits_max[idx_cells(j, i)];
58 if (distance > 0) {
59 return Position::below;
60 }
61 distance = 0.0;
62 return Position::at;
63 }
64
65 void getAlphaPatch(const double* const x, double* alpha) {
67 kernels::idx4 idx(basisSize, basisSize, basisSize, Dimensions);
68
69 constexpr int center = basisSize / 2;
70
71 // find cell id on surface of patch
72 int cell_j = static_cast<int>(std::floor((x[idx(center, center, center, 2)] - info->domainOffset[2]) / dx));
73 int cell_i = static_cast<int>(std::floor((x[idx(center, center, center, 0)] - info->domainOffset[0]) / dx));
74
75 double patch_distance = 0;
76 // position 0: above 1: at 2: below topography
77 int position = findRelativePosition(
78 patch_distance, x[idx(0, 0, 0, 1)], x[idx(0, basisSize - 1, 0, 1)], cell_i, cell_j
79 );
80
81 // compute distance of patch to underlaying surface
82 if (patch_distance >= max_distance_to_direct_surface) {
83 if (position == Position::above) {
84 std::fill_n(alpha, basisSize * basisSize * basisSize, 0.0);
85 return;
86 }
87 if (position == Position::below) {
88 std::fill_n(alpha, basisSize * basisSize * basisSize, 1.0);
89 return;
90 }
91 }
92
93 double cell_dx = (x[idx(0, basisSize - 1, 0, 1)] - x[idx(0, 0, 0, 1)]);
94 // add distance from top center to edge as buffer
95 patch_distance += sqrt(3.0 / 2.0) * dx;
96
97 // minmal measuerd distance for each patch node
98 double minDistance[basisSize * basisSize * basisSize];
99 std::fill_n(minDistance, basisSize * basisSize * basisSize, patch_distance);
100 // boolean if the patch is done
101 bool done[basisSize * basisSize * basisSize];
102 std::fill_n(done, basisSize * basisSize * basisSize, false);
103
104 Position positions[basisSize * basisSize * basisSize];
105 getMinDistancePatch(patch_distance, x, minDistance, positions, done);
106 computeAlphaPatch(alpha, minDistance, positions);
107
108 // std::copy_n(minDistance,basisSize*basisSize*basisSize,alpha);
109 }
110
111 void computeAlphaPatch(double* alpha, double* minDistance, Position* positions) {
112 kernels::idx3 idx_alpha(basisSize, basisSize, basisSize);
113 for (int node_z = 0; node_z < basisSize; node_z++) {
114 for (int node_y = 0; node_y < basisSize; node_y++) {
115 for (int node_x = 0; node_x < basisSize; node_x++) {
116 Position position_node = positions[idx_alpha(node_z, node_y, node_x)];
117
118 double r = minDistance[idx_alpha(node_z, node_y, node_x)] / interval_size;
119 if (r < 0.5) {
120 if (position_node == Position ::above)
121 r = -r;
122 alpha[idx_alpha(node_z, node_y, node_x)] = (std::sin(r * M_PI) + 1) * 0.5;
123 } else {
124 alpha[idx_alpha(node_z, node_y, node_x)] = position_node == Position::below ? 1 : 0;
125 }
126#ifdef Asserts
128 alpha[idx_alpha(node_z, node_y, node_x)] >= 0 && alpha[idx_alpha(node_z, node_y, node_x)] <= 1,
129 alpha[idx_alpha(node_z, node_y, node_x)]
130 );
131#endif
132 }
133 }
134 }
135 }
136
138 std::vector<double>& patch_x, std::vector<double>& patch_z, double center_x, double center_z, double patch_radius
139 ) {
140
141 kernels::idx4 idx(basisSize, basisSize, basisSize, Dimensions);
142 kernels::idx3 idx_node(basisSize, basisSize, basisSize);
143
144 constexpr int center = basisSize / 2;
145
146 // Todo: Is this basis size a valid coice?
147 int num_circles = basisSize;
148 std::vector<double> radius;
149 radius.resize(num_circles);
150
151 int num_args[num_circles];
152 num_args[0] = 1;
153 num_args[1] = basisSize;
154
155 double delta_radius = patch_radius / static_cast<double>(num_circles - 1);
156 for (int i = 0; i < num_circles; i++) {
157 radius[i] = delta_radius * i;
158 }
159
160 int node_counter = 0;
161
162 patch_x.resize(num_args[0]);
163 patch_z.resize(num_args[0]);
164
165 patch_x[node_counter] = center_x;
166 patch_z[node_counter] = center_z;
167
168 node_counter++;
169
170 for (int j = 1; j < num_circles; j++) {
171 num_args[j] = static_cast<int>(std::ceil(radius[j] / radius[1])) * num_args[1];
172 patch_x.resize(patch_x.size() + num_args[j]);
173 patch_z.resize(patch_z.size() + num_args[j]);
174 for (int i = 0; i < num_args[j]; i++) {
175 double arg = 2 * M_PI * i / static_cast<double>(num_args[j]);
176 patch_x[node_counter] = std::cos(arg) * radius[j] + center_x;
177 patch_z[node_counter] = std::sin(arg) * radius[j] + center_z;
178
179 /*truncate on domain boundary*/
180 patch_x[node_counter] = std::min(
181 std::max(patch_x[node_counter], info->domainOffset[0]), info->domainOffset[0] + info->domainSize[0]
182 );
183 patch_z[node_counter] = std::min(
184 std::max(patch_z[node_counter], info->domainOffset[2]), info->domainOffset[2] + info->domainSize[2]
185 );
186 node_counter++;
187 }
188 }
189 }
190
192 double& min_distance,
193 Position& position,
194 const double* const x,
195 std::vector<double>& patch_x,
196 std::vector<double>& patch_y,
197 std::vector<double>& patch_z
198 ) {
199
200 double new_min_distance = std::numeric_limits<double>::max();
201 for (int i = 0; i < patch_x.size(); i++) {
202 double distance = std::sqrt(
203 (patch_x[i] - x[0]) * (patch_x[i] - x[0]) + (patch_y[i] - x[1]) * (patch_y[i] - x[1])
204 + (patch_z[i] - x[2]) * (patch_z[i] - x[2])
205 );
206 if (distance < new_min_distance) {
207 new_min_distance = distance;
208 position = (patch_y[i] - x[1]) > 0 ? Position::above : Position::below;
209 }
210 }
211 min_distance = new_min_distance;
212 }
213
215 double patch_distance, const double* const x, double* min_distance, Position* positions, bool* done
216 ) {
217
218 /*span patch around center of the cell*/
219 std::vector<double> patch_x;
220 std::vector<double> patch_z;
221
222 kernels::idx4 idx(basisSize, basisSize, basisSize, Dimensions);
223 kernels::idx3 idx_node(basisSize, basisSize, basisSize);
224
225 double center_x = (x[idx(0, 0, basisSize - 1, 0)] + x[idx(0, 0, 0, 0)]) / 2.0;
226 double center_z = (x[idx(basisSize - 1, 0, 0, 2)] + x[idx(0, 0, 0, 2)]) / 2.0;
227 spanPatchAroundCenter(patch_x, patch_z, center_x, center_z, patch_distance);
228
229 /*evaluate topography on patch*/
230 std::vector<double> patch_y;
231 topography->topography_onSet(patch_x, patch_z, patch_y);
232
233 /*compute distance to surface for each node*/
234
235 for (int node_z = 0; node_z < basisSize; node_z++) {
236 for (int node_y = 0; node_y < basisSize; node_y++) {
237 for (int node_x = 0; node_x < basisSize; node_x++) {
238 double new_min_distance = 0;
239 Position new_position;
241 new_min_distance, new_position, x + idx(node_z, node_y, node_x, 0), patch_x, patch_y, patch_z
242 );
243
244 // If new min distance doesn't change we are done
245 if (new_min_distance < min_distance[idx_node(node_z, node_y, node_x)]) {
246 positions[idx_node(node_z, node_y, node_x)] = new_position;
247 min_distance[idx_node(node_z, node_y, node_x)] = new_min_distance;
248 } else {
249 done[idx_node(node_z, node_y, node_x)] = true;
250 }
251 }
252 }
253 }
254
255 bool all_done = true;
256 for (int i = 0; i < basisSize * basisSize * basisSize; i++) {
257 if (!done[i]) {
258 all_done = false;
259 break;
260 }
261 }
262
263 /* if(!isDG){ */
264 /* for(int node_z = 0; node_z < basisSize; node_z++){ */
265 /* for(int node_y = 0; node_y < basisSize; node_y++){ */
266 /* for(int node_x = 0; node_x < basisSize; node_x++){ */
267 /* std::cout << min_distance[idx_node(node_z,node_y,node_x)] << std::endl; */
268 /* } */
269 /* } */
270 /* } */
271 /* std::cout << " "<< std::endl; */
272 /* } */
273
274 if (all_done) { // no closer nodes found -> we are done
275 return;
276 } else {
277 double patch_distance = std::numeric_limits<double>::max();
278
279 for (int node_z = 0; node_z < basisSize; node_z++) {
280 for (int node_y = 0; node_y < basisSize; node_y++) {
281 for (int node_x = 0; node_x < basisSize; node_x++) {
282
283 if (!done[idx_node(node_z, node_y, node_x)]) {
284 double node_distance = 0.0;
285 node_distance = std::sqrt(
286 (x[idx(node_z, node_y, node_x, 0)] - center_x
287 ) * (x[idx(node_z, node_y, node_x, 0)] - center_x)
288 + (x[idx(node_z, node_y, node_x, 2)] - center_z
289 ) * (x[idx(node_z, node_y, node_x, 2)] - center_z)
290 )
291 + min_distance[idx_node(node_z, node_y, node_x)];
292
293 patch_distance = std::min(node_distance, patch_distance);
294 }
295 }
296 }
297 }
298 getMinDistancePatch(patch_distance, x, min_distance, positions, done);
299 return;
300 }
301 }
302
303 double getAlpha(const double* const x) {
305 int cell_j = static_cast<int>(std::floor((x[2] - info->domainOffset[2]) / dx));
306 int cell_i = static_cast<int>(std::floor((x[0] - info->domainOffset[0]) / dx));
307
308 double min_dist_limit = topography_limits_min[idx2(cell_j, cell_i)] - x[1];
309 if (min_dist_limit > max_distance_to_direct_surface) {
310 return 0;
311 }
312
313 double max_dist_limit = x[1] - topography_limits_max[idx2(cell_j, cell_i)];
314 if (max_dist_limit > max_distance_to_direct_surface) {
315 return 1;
316 }
317
318 double min_distance = std::min(std::abs(max_dist_limit), std::abs(min_dist_limit));
319 min_distance = getMinDistance(min_distance, x[0], x[1], x[2]);
320
321 double alpha;
322 if (std::abs(min_distance) < interval_size * 0.5) {
323 double r = min_distance / interval_size;
324 alpha = (std::sin(r * M_PI) + 1) * 0.5;
325 } else {
326 alpha = min_distance < 0 ? 0 : 1;
327 }
328 return alpha;
329 };
330
331 double getMinDistance(double max_distance, double x, double y, double z) {
332
333 std::vector<double> patch_x;
334 std::vector<double> patch_z;
335 std::vector<double> patch_y;
336 std::vector<double> radius;
337
338 int num_circles = 5; // Todo: Is this a valid choice?
339 int num_args[num_circles];
340 num_args[0] = 1;
341 num_args[1] = 10; // Todo: Is this a valid coice?
342
343 radius.resize(num_circles);
344
345 double delta_radius = std::abs(max_distance) / static_cast<double>(num_circles - 1);
346 for (int i = 0; i < num_circles; i++) {
347 radius[i] = delta_radius * i;
348 }
349
350 int node_counter = 0;
351
352 patch_x.resize(1);
353 patch_z.resize(1);
354 patch_x[0] = x;
355 patch_z[0] = z;
356
357 node_counter++;
358
359 for (int j = 1; j < num_circles; j++) {
360 num_args[j] = static_cast<int>(std::ceil(radius[j] / radius[1])) * num_args[1];
361 patch_x.resize(patch_x.size() + num_args[j]);
362 patch_z.resize(patch_z.size() + num_args[j]);
363 for (int i = 0; i < num_args[j]; i++) {
364 double arg = 2 * M_PI * i / static_cast<double>(num_args[j]);
365 patch_x[node_counter] = std::cos(arg) * radius[j] + x;
366 patch_z[node_counter] = std::sin(arg) * radius[j] + z;
367
368 patch_x[node_counter] = std::min(
369 std::max(patch_x[node_counter], info->domainOffset[0]), info->domainOffset[0] + info->domainSize[0]
370 );
371 patch_z[node_counter] = std::min(
372 std::max(patch_z[node_counter], info->domainOffset[2]), info->domainOffset[2] + info->domainSize[2]
373 );
374
375 node_counter++;
376 }
377 }
378
379 topography->topography_onSet(patch_x, patch_z, patch_y);
380
381 double new_max_distance = max_distance;
382 double distance;
383 node_counter = 0;
384 for (int j = 0; j < num_circles; j++) {
385 for (int i = 0; i < num_args[j]; i++) {
386 distance = std::sqrt((patch_y[node_counter] - y) * (patch_y[node_counter] - y) + radius[j] * radius[j]);
387 if (std::abs(new_max_distance) > std::abs(distance)) {
388 new_max_distance = distance;
389 }
390 node_counter++;
391 }
392 }
393
394 if (std::abs(new_max_distance - max_distance) < 1.0e-10) { // no closer nodes found
395 return max_distance;
396 } else {
397 return getMinDistance(new_max_distance, x, y, z);
398 }
399 }
400
401 void initLimits() {
402 std::vector<double> x;
403 std::vector<double> y;
404 std::vector<double> z;
405
408
409 x.resize(basisSize);
410 z.resize(basisSize);
411
413 for (int j = 0; j < domain_cells; j++) {
414 for (int i = 0; i < domain_cells; i++) {
415 for (int node = 0; node < basisSize; node++) {
416 if (isDG) {
417 x[node] = kernels::legendre::nodes[basisSize - 1][node] * dx + i * dx + info->domainOffset[0];
418 z[node] = kernels::legendre::nodes[basisSize - 1][node] * dx + j * dx + info->domainOffset[2];
419 } else {
420 x[node] = (node + 1) / static_cast<double>(basisSize + 1) * dx + i * dx + info->domainOffset[0];
421 z[node] = (node + 1) / static_cast<double>(basisSize + 1) * dx + j * dx + info->domainOffset[2];
422 }
423 }
424
425 topography->topography_onPatch(x, z, y);
426
427 double min = std::numeric_limits<double>::max();
428 double max = -std::numeric_limits<double>::max();
429
430 for (int node = 0; node < basisSize * basisSize; node++) {
431 min = std::min(min, y[node]);
432 max = std::max(max, y[node]);
433 }
434 topography_limits_max[idx2(j, i)] = max;
435 topography_limits_min[idx2(j, i)] = min;
436 }
437 }
438 }
439
440private:
441 Topography* topography;
443
444 bool isDG;
445
446 double dx;
451
452 std::vector<double> topography_limits_max;
453 std::vector<double> topography_limits_min;
454 std::vector<double> integration_nodes;
455
457};
458
459#endif
#define assertion1(expr, param)
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
Position
@ at
@ above
@ below
Topography * topography
double getAlpha(const double *const x)
std::vector< double > integration_nodes
void computeAlphaPatch(double *alpha, double *minDistance, Position *positions)
int findRelativePosition(double &distance, double y_min, double y_max, int i, int j)
std::vector< double > topography_limits_min
void getMinDistancePatch(double patch_distance, const double *const x, double *min_distance, Position *positions, bool *done)
DomainInformation * info
double max_distance_to_direct_surface
double getMinDistance(double max_distance, double x, double y, double z)
void spanPatchAroundCenter(std::vector< double > &patch_x, std::vector< double > &patch_z, double center_x, double center_z, double patch_radius)
void findMinDistanceToNode(double &min_distance, Position &position, const double *const x, std::vector< double > &patch_x, std::vector< double > &patch_y, std::vector< double > &patch_z)
DiffuseInterface(Topography *a_topo, DomainInformation *a_domainInfo, SolverInformation *a_solverInfo, double a_interface_factor, double a_max_distance_to_direct_surface)
void getAlphaPatch(const double *const x, double *alpha)
std::vector< double > topography_limits_max
virtual bool isDG()=0