Peano 4
Loading...
Searching...
No Matches
DiffuseInterface_interpolate.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
7template <class Shortcuts, int basisSize>
9public:
10 DiffuseInterface(Topography* a_topo, DomainInformation* a_info, int a_interface_factor) {
11 interface_factor = a_interface_factor;
12 topography = a_topo;
13 info = a_info;
14 isDG = a_info->isDG();
15
16 kernels::initGaussLegendreNodesAndWeights(std::set<int>()); // empty set as it is not used
17
18 info->getDx(dx);
19
20 interval_size = dx / basisSize * interface_factor;
21 max_patch_size = (dx * 5.0);
22
23 _nx = int(std::ceil(info->domainSize[0] / dx));
24 _nz = int(std::ceil(info->domainSize[2] / dx));
25
26 integration_nodes.resize(basisSize);
27
28 topo_x.resize(_nx);
29 topo_z.resize(_nz);
30 topo_y.resize(_nx * _nz);
31
32 normals_x.resize(_nx * _nz);
33 normals_y.resize(_nx * _nz);
34 normals_z.resize(_nx * _nz);
35
36 for (int node = 0; node < basisSize; node++) {
37 if (isDG) {
38 integration_nodes[node] = kernels::legendre::nodes[basisSize - 1][node];
39 } else {
40 integration_nodes[node] = (node + 1) / (basisSize + 1);
41 }
42 }
43 initCells();
44 };
45
46 void getAlphaPatch(const double* const x, double* alpha) {
47 kernels::idx2 idx2(_nz, _nx);
48 int j = static_cast<int>(std::floor((x[2] - info->domainOffset[2]) / dx));
49 int i = static_cast<int>(std::floor((x[0] - info->domainOffset[0]) / dx));
50
51 kernels::idx4 id_xyzn(basisSize, basisSize, basisSize, Dimensions);
52
53 double min_dist_limit = topography_limits_min[idx2(cell_j, cell_i)] - x[id_xyzn(0, basisSize - 1, 0, 1)];
54 if (min_dist_limit > max_patch_size) {
55 std::fill_n(alpha, 0.0, basisSize * basisSize * basisSize);
56 return;
57 }
58
59 double max_dist_limit = x[id_xyzn(0, 0, 0, 1)] - topography_limits_max[idx2(cell_j, cell_i)];
60 if (max_dist_limit > max_patch_size) {
61 std::fill_n(alpha, 0.0, basisSize * basisSize * basisSize);
62 return;
63 }
64
65 // identify position of topography relative to cell
66 for (int node_z = 0; node_z < basisSize; node_z++) {
67 for (int node_x = 0; node_x < basisSize; node_x++) {
68 topo_y[idy(i, node_z, j, node_x)];
69 }
70 }
71 }
72
73 void initCells() {
76
77 kernels::idx2 id_cells(_nz, _nx);
78 kernels::idx2 idx(_nx, basisSize);
79 kernels::idx2 idz(_nz, basisSize);
80 kernels::idx3 idy(_nz, basisSize, _nx, basisSize);
81
82 for (int i = 0; i < _nx; i++) {
83 for (int node = 0; node < basisSize; node++) {
84 if (isDG) {
85 topo_x[idx(i, node)] = kernels::legendre::nodes[basisSize - 1][node] * dx + i * dx + info->domainOffset[0];
86 } else {
87 topo_x[idx(i, node)] = (node + 1) / static_cast<double>(basisSize + 1) * dx + i * dx + info->domainOffset[0];
88 }
89 }
90 }
91
92 for (int i = 0; i < _nz; i++) {
93 for (int node = 0; node < basisSize; node++) {
94 if (isDG) {
95 topo_z[idz(i, node)] = kernels::legendre::nodes[basisSize - 1][node] * dx + i * dx + info->domainOffset[2];
96 } else {
97 topo_z[idz(i, node)] = (node + 1) / static_cast<double>(basisSize + 1) * dx + i * dx + info->domainOffset[2];
98 }
99 }
100 }
101
102 topography->topography_onPatch(topo_x, topo_z, topo_y);
103 for (int j = 0; j < _nz; j++) {
104 for (int i = 0; i < _nx; i++) {
105 double min = std::numeric_limits<double>::max();
106 double max = -std::numeric_limits<double>::max();
107 for (int node_z = 0; node_z < basisSize; node_z++) {
108 for (int node_x = 0; node_x < basisSize; node_x++) {
109 min = std::min(min, y[idy(j, node_z, i, node_x)]);
110 max = std::max(max, y[idy(j, node_z, i, node_x)]);
111 }
112 }
113 topography_limits_max[id_cells(j, i)] = max;
114 topography_limits_min[id_cells(j, i)] = min;
115 }
116 }
117
118 // init normals
119 double grad_x[basisSize];
120 for (int j = 0; j < _nz; j++) {
121 for (int i = 0; i < _nx; i++) {
122 kernels::idx2 idx_grad(basisSize, basisSize);
123 double grad_x[basisSize * basisSize];
124 double grad_z[basisSize * basisSize];
125 // compute gradient
126 for (int node_z = 0; node_z < basisSize; node_z++) {
127 for (int node_x = 0; node_x < basisSize; node_x++) {
128 for (int n = 0; n < basisSize; n++) {
129 grad_x[idx_grad(
130 node_z, node_x
131 )] += kernels::legendre::dudx[num_nodes - 1][node_x][n] * topo_y[idy(j, node_z, i, n)];
132 grad_z[idx_grad(
133 node_z, node_x
134 )] += kernels::legendre::dudx[num_nodes - 1][node_z][n] * topo_y[idy(j, n, i, node_x)];
135 }
136 }
137 }
138 }
139 }
140 }
141
142private:
143 Topography* topography;
145
146 bool isDG;
147 double dx;
148
149 double interval_size;
150 double interface_factor;
151
153
154 int _nx;
155 int _nz;
156
157 std::vector<double> topography_limits_max;
158 std::vector<double> topography_limits_min;
159
160 std::vector<double> integration_nodes;
161
162 std::vector<double> topo_x;
163 std::vector<double> topo_y;
164 std::vector<double> topo_z;
165
166 std::vector<double> normals_x;
167 std::vector<double> normals_y;
168 std::vector<double> normals_z;
169};
170
171#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
std::vector< double > normals_y
Topography * topography
std::vector< double > topo_x
std::vector< double > integration_nodes
std::vector< double > normals_x
std::vector< double > topo_z
std::vector< double > topography_limits_min
DiffuseInterface(Topography *a_topo, DomainInformation *a_info, int a_interface_factor)
DomainInformation * info
void getAlphaPatch(const double *const x, double *alpha)
std::vector< double > normals_z
std::vector< double > topography_limits_max
std::vector< double > topo_y