Peano 4
Loading...
Searching...
No Matches
DiffuseInterface_Convolution.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) {
11 topography = a_topo;
12 info = a_info;
13
14 int nx = std::pow(3, info->meshLevel);
15 int ny = std::pow(3, info->meshLevel);
16 // TODO: how do i get a good approximation for the real buffered nz ?
17 int nz = std::pow(3, info->meshLevel);
18
19 std::vector<double> alpha;
20 alpha.resize(nx * ny * nz * (basisSize) * (basisSize) * (basisSize));
21 kernels::initGaussLegendreNodesAndWeights(std::set<int>()); // empty set as it is not used
22
23 info->getDx(dx);
24
25 patch_size = (dx * 2.0) / basisSize;
26
27 // patch_size = (dx * 2.0 );
28 patch_dx = dx / 3.0;
29 patch_cells = int(std::ceil(patch_size / patch_dx));
30 domain_cells = int(std::ceil(info->domainSize[0] / dx));
31
32 integration_nodes.resize(basisSize);
33 integration_weights.resize(basisSize);
34
35 for (int node = 0; node < basisSize; node++) {
36 integration_nodes[node] = kernels::legendre::nodes[basisSize - 1][node] * patch_dx;
37 integration_weights[node] = kernels::legendre::weights[basisSize - 1][node] * patch_dx;
38 }
39
40 patch.resize((patch_cells * basisSize) * (patch_cells * basisSize) * (patch_cells * basisSize));
41
42 kernels::idx6 idx(patch_cells, patch_cells, patch_cells, basisSize, basisSize, basisSize);
43
45
46 for (int cell_z = 0; cell_z < patch_cells; cell_z++) {
47 double offset_z = (cell_z - patch_cells / 2.0) * patch_dx;
48 for (int cell_y = 0; cell_y < patch_cells; cell_y++) {
49 double offset_y = (cell_y - patch_cells / 2.0) * patch_dx;
50 for (int cell_x = 0; cell_x < patch_cells; cell_x++) {
51 double offset_x = (cell_x - patch_cells / 2.0) * patch_dx;
52 for (int node_z = 0; node_z < basisSize; node_z++) {
53 double z = offset_z + integration_nodes[node_z];
54 for (int node_y = 0; node_y < basisSize; node_y++) {
55 double y = offset_y + integration_nodes[node_y];
56 for (int node_x = 0; node_x < basisSize; node_x++) {
57 double x = offset_x + integration_nodes[node_x];
58 patch[idx(cell_z, cell_y, cell_x, node_z, node_y, node_x)] = testFunction(x, y, z);
59
60 volume_testfunction += patch[idx(cell_z, cell_y, cell_x, node_z, node_y, node_x)]
62 * integration_weights[node_x];
63 }
64 }
65 }
66 }
67 }
68 }
69 initLimits();
70 };
71
72 double getAlpha(const double* const x) {
73 std::vector<double> patch_x;
74 std::vector<double> patch_z;
75 std::vector<double> patch_y;
76 int cell_index = int(std::floor((x[2] - info->domainOffset[2]) / dx)) * domain_cells
77 + int(std::floor((x[0] - info->domainOffset[0]) / dx));
78
79 if (x[1] < topography_limits_min[cell_index] - dx) {
80 // std::cout << "." << std::endl;
81 return 0;
82 }
83
84 if (x[1] > topography_limits_max[cell_index] + dx) {
85 // std::cout << "." << std::endl;
86 return 1;
87 }
88
89 patch_x.resize(patch_cells * basisSize);
90 patch_z.resize(patch_cells * basisSize);
91
92 for (int cell = 0; cell < patch_cells; cell++) {
93 for (int node = 0; node < basisSize; node++) {
94 patch_z[cell * basisSize + node] = x[2] + (cell - patch_cells / 2.0) * patch_dx + integration_nodes[node];
95 patch_x[cell * basisSize + node] = x[0] + (cell - patch_cells / 2.0) * patch_dx + integration_nodes[node];
96 }
97 }
98
99 topography->topography_onPatch(patch_x, patch_z, patch_y);
100
101 std::vector<double> topography_heaviside;
102 gen_heaviside_topography(x, patch_y, topography_heaviside);
103 double alpha = convolution_with_patch(topography_heaviside);
104 return alpha;
105 };
106
107 double testFunction(double x, double y, double z) {
108 double r2 = x * x + y * y + z * z;
109 double r2_patch = patch_size * patch_size;
110 if (r2 > r2_patch) {
111 return 0;
112 } else {
113 return 1.0 / (r2 + 1);
114 // return std::exp(r2_patch/(r2-r2_patch));
115 }
116 }
117
119 const double* const x, std::vector<double>& patch_y, std::vector<double>& topography_heaviside
120 ) {
121 kernels::idx6 idx(patch_cells, patch_cells, patch_cells, basisSize, basisSize, basisSize);
122 kernels::idx2 idx2(patch_cells * basisSize, patch_cells * basisSize);
123
124 topography_heaviside.resize(patch_cells * basisSize * patch_cells * basisSize * patch_cells * basisSize);
125
126 for (int cell_z = 0; cell_z < patch_cells; cell_z++) {
127 for (int node_z = 0; node_z < basisSize; node_z++) {
128 for (int cell_x = 0; cell_x < patch_cells; cell_x++) {
129 for (int node_x = 0; node_x < basisSize; node_x++) {
130 for (int cell_y = 0; cell_y < patch_cells; cell_y++) {
131 double offset_y = x[1] + (cell_y - patch_cells / 2.0) * patch_dx;
132 for (int node_y = 0; node_y < basisSize; node_y++) {
133 double y = offset_y + integration_nodes[node_y];
134 topography_heaviside[idx(
135 cell_z, cell_y, cell_x, node_z, node_y, node_x
136 )] = patch_y[idx2(cell_z * basisSize + node_z, cell_x * basisSize + node_x)] < y ? 1 : 0;
137 }
138 }
139 }
140 }
141 }
142 }
143 }
144
145 double convolution_with_patch(std::vector<double>& topography_heaviside) {
146 double alpha = 0;
147 kernels::idx6 idx(patch_cells, patch_cells, patch_cells, basisSize, basisSize, basisSize);
148 for (int cell_z = 0; cell_z < patch_cells; cell_z++) {
149 for (int node_z = 0; node_z < basisSize; node_z++) {
150 for (int cell_x = 0; cell_x < patch_cells; cell_x++) {
151 for (int node_x = 0; node_x < basisSize; node_x++) {
152 for (int cell_y = 0; cell_y < patch_cells; cell_y++) {
153 for (int node_y = 0; node_y < basisSize; node_y++) {
154 double weight = integration_weights[node_z] * integration_weights[node_y] * integration_weights[node_x];
155
156 alpha += topography_heaviside[idx(cell_z, cell_y, cell_x, node_z, node_y, node_x)]
157 * patch[idx(cell_z, cell_y, cell_x, node_z, node_y, node_x)] * weight;
158 }
159 }
160 }
161 }
162 }
163 }
164 return alpha / volume_testfunction;
165 }
166
167 void initLimits() {
168 std::vector<double> x;
169 std::vector<double> y;
170 std::vector<double> z;
171
174
175 x.resize(basisSize);
176 z.resize(basisSize);
177
178 for (int i = 0; i < domain_cells; i++) {
179 for (int j = 0; j < domain_cells; j++) {
180 for (int node = 0; node < basisSize; node++) {
181 x[node] = kernels::legendre::nodes[basisSize - 1][node] * dx + i * dx + info->domainOffset[0];
182 z[node] = kernels::legendre::nodes[basisSize - 1][node] * dx + j * dx + info->domainOffset[2];
183 }
184
185 topography->topography_onPatch(x, z, y);
186
187 double min = std::numeric_limits<double>::max();
188 double max = std::numeric_limits<double>::min();
189 for (int node = 0; node < basisSize * basisSize; node++) {
190 min = std::min(min, y[node]);
191 max = std::max(max, y[node]);
192 }
195 }
196 }
197 }
198
199 int getArea(int x, int z) {
200 int cell_x = (x - info->domainOffset[0]) / info->domainSize[0];
201 int cell_z = (z - info->domainOffset[2]) / info->domainSize[2];
202 return
203 }
204
205private:
206 Topography* topography;
208
209 double dx;
210
212 double patch_dx;
213
215 int domain_cells;
216
217 std::vector<double> topography_limits_max;
218 std::vector<double> topography_limits_min;
219 std::vector<double> patch;
220 std::vector<double> integration_nodes;
221 std::vector<double> integration_weights;
222
223 double volume_testfunction;
224};
225
226#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 COMMENT THIS OUT. I DON 'T THINK IT 'S READY YET Now we are working on a cell-by-cell basis, let us examine the @f$ 1^{st} @f$ cell. We note that only @f$ \phi_0 @f$ and @f$ \phi_1 @f$ are non-zero here, so our original equation reduces to:\f{eqnarray *}{ u_1 \ \sum_{k=1, 2} \int_{c_k}(\nabla \phi_0, \ \phi_0) \+dx \=\ \int_{c_0}(f, \ \phi_0) dx \f} and we drop the @f$ \forall \phi @f$ requirement since only @f$ \phi_0 @f$ is non-zero here. In effect, this @f$ \int_{c_0}(\nabla \phi_0, \ \phi_0) dx \=\ \int_{c_0}(f, \ \phi_0) dx @f$ will become our matrix element, as we shall see shortly. --> The basis of functions that we use in this example are piecewise linear defined in each cell(which is an interval of width @f$ h @f$)
Topography * topography
std::vector< double > patch
double getAlpha(const double *const x)
std::vector< double > integration_nodes
std::vector< double > topography_limits_min
void gen_heaviside_topography(const double *const x, std::vector< double > &patch_y, std::vector< double > &topography_heaviside)
DomainInformation * info
double convolution_with_patch(std::vector< double > &topography_heaviside)
double testFunction(double x, double y, double z)
std::vector< double > integration_weights
DiffuseInterface(Topography *a_topo, DomainInformation *a_info)
std::vector< double > topography_limits_max