Peano 4
Loading...
Searching...
No Matches
ContextCurvilinear.h
Go to the documentation of this file.
1#ifndef EXASEIS_CONTEXT_CURVILINEAR_HEADER
2#define EXASEIS_CONTEXT_CURVILINEAR_HEADER
3
4// #include "kernels/GaussLegendreBasis.h"
5#include "Context.h"
7#include "curvi/coordinate.h"
9#include "curvi/geometry/block.h"
10#include "curvi/interface/interface.h"
11#include "curvi/kdTree/root.h"
13
14template <class Shortcuts, int basisSize>
15class ContextCurvilinear: public Context<Shortcuts, basisSize> {
16
17public:
19 std::string& scenario_string,
20 std::string& a_topography_string,
21 DomainInformation* a_domain_info,
23 ):
24 Context<Shortcuts, basisSize>(scenario_string, a_domain_info) {
25
26 domain_info = a_domain_info;
27 solver_info = a_solver_info;
28 topography_string = a_topography_string;
29 }
30
32 uint elements_l[3];
33 elements_l[Coordinate::X] = domain_info->elements[0];
34 elements_l[Coordinate::Y] = domain_info->elements[1];
35 elements_l[Coordinate::Z] = domain_info->elements[2];
36
37 double size[3];
38 size[Coordinate::X] = domain_info->domainSize[0];
39 size[Coordinate::Y] = domain_info->domainSize[1];
40 size[Coordinate::Z] = domain_info->domainSize[2];
41
42 double offset[3];
43 offset[Coordinate::X] = domain_info->domainOffset[0];
44 offset[Coordinate::Y] = domain_info->domainOffset[1];
45 offset[Coordinate::Z] = domain_info->domainOffset[2];
46
47 interface = new Curvi::Interface(topography_string, offset, size, elements_l, basisSize - 1);
48 }
49
50 Root* getRoot() { return interface->getRoot(); }
51
52 void initTree() { interface->initTree(); }
53
55 delete interface;
56 delete domain_info;
57 delete solver_info;
58 // delete topography;
59 }
60
61 virtual void initUnknownsPatch(
62 double* luh,
65 double t,
66 double dt
67 ) {
68
69 if (tarch::la::equals(t, 0.0)) {
70 Shortcuts s;
71 int numberOfData = s.SizeVariables + s.SizeParameters;
72 kernels::idx4 id_xyzf(basisSize, basisSize, basisSize, numberOfData);
73 kernels::idx3 id_xyz(basisSize, basisSize, basisSize);
74
75 int ez = std::round((center[2] - domain_info->domainOffset[2] - dx[2] / 2.0) / dx[2]);
76 int ey = std::round((center[1] - domain_info->domainOffset[1] - dx[1] / 2.0) / dx[1]);
77 int ex = std::round((center[0] - domain_info->domainOffset[0] - dx[0] / 2.0) / dx[0]);
78
79#ifndef QUICK_IDENTITY
80 Block element = interface->getElement(solver_info->nodes, ez, ey, ex);
81 int index_offset[3];
82 element.getIndexOffset(index_offset);
83
84 for (int k = 0; k < basisSize; k++) {
85 int e_k = index_offset[Coordinate::Z] + k;
86 for (int j = 0; j < basisSize; j++) {
87 int e_j = index_offset[Coordinate::Y] + j;
88 for (int i = 0; i < basisSize; i++) {
89 int e_i = index_offset[Coordinate::X] + i;
90 // x,y,z
91 luh[id_xyzf(k, j, i, s.curve_grid + 0)] = element(
92 2, Coordinate::Z, e_k, Coordinate::Y, e_j, Coordinate::X, e_i
93 );
94 luh[id_xyzf(k, j, i, s.curve_grid + 1)] = element(
95 1, Coordinate::Z, e_k, Coordinate::Y, e_j, Coordinate::X, e_i
96 );
97 luh[id_xyzf(k, j, i, s.curve_grid + 2)] = element(
98 0, Coordinate::Z, e_k, Coordinate::Y, e_j, Coordinate::X, e_i
99 );
100 }
101 }
102 }
103
104#else
105 for (int k = 0; k < basisSize; k++) {
106 for (int j = 0; j < basisSize; j++) {
107 for (int i = 0; i < basisSize; i++) {
108 luh[id_xyzf(k, j, i, s.curve_grid + 0)] = solver_info->nodes[i] * dx[0] + center[0] - dx[0] / 2;
109 luh[id_xyzf(k, j, i, s.curve_grid + 1)] = solver_info->nodes[j] * dx[1] + center[1] - dx[1] / 2;
110 luh[id_xyzf(k, j, i, s.curve_grid + 2)] = solver_info->nodes[k] * dx[2] + center[2] - dx[2] / 2;
111 }
112 }
113 }
114#endif
115
116 double derivatives[basisSize * basisSize * basisSize * 10];
117 std::fill_n(derivatives, basisSize * basisSize * basisSize * 10, 0.0);
118 kernels::idx4 id_der(basisSize, basisSize, basisSize, 10);
119
120 // compute metric derivatives//
121#ifdef QUICK_IDENTITY
122
123 for (int k = 0; k < basisSize; k++) {
124 for (int j = 0; j < basisSize; j++) {
125 for (int i = 0; i < basisSize; i++) {
126 derivatives[id_der(k, j, i, 0)] = 1.0;
127 derivatives[id_der(k, j, i, 1)] = 1.0;
128 derivatives[id_der(k, j, i, 5)] = 1.0;
129 derivatives[id_der(k, j, i, 9)] = 1.0;
130 }
131 }
132 }
133#else
135#endif
136
137 for (int k = 0; k < basisSize; k++) {
138 for (int j = 0; j < basisSize; j++) {
139 for (int i = 0; i < basisSize; i++) {
140 std::copy_n(derivatives + id_der(k, j, i, 0), 10, luh + id_xyzf(k, j, i, s.jacobian));
141 }
142 }
143 }
144
146 getElementCenter(luh, curv_center);
147 for (int k = 0; k < basisSize; k++) {
148 for (int j = 0; j < basisSize; j++) {
149 for (int i = 0; i < basisSize; i++) {
150 double coords[3] = {
151 luh[id_xyzf(k, j, i, s.curve_grid + 0)],
152 luh[id_xyzf(k, j, i, s.curve_grid + 1)],
153 luh[id_xyzf(k, j, i, s.curve_grid + 2)]};
154 this->scenario->initUnknownsPointwise(coords, curv_center, t, dt, luh + id_xyzf(k, j, i, 0));
155 }
156 }
157 }
158
159 for (int k = 0; k < basisSize; k++) {
160 for (int j = 0; j < basisSize; j++) {
161 for (int i = 0; i < basisSize; i++) {
162 for (int m = 0; m < numberOfData; m++) {
163 assertion5(std::isfinite(luh[id_xyzf(k, j, i, m)]), k, j, i, m, domain_info->meshLevel);
164 }
165 }
166 }
167 }
168 }
169 }
170
171 void initPointSourceLocation(double pointSourceLocation[][3]) override {
172 this->scenario->initPointSourceLocation(pointSourceLocation);
173
174 for (int p = 0; p < 1; p++) {
175 Coordinate dir_top = static_cast<Coordinate>(_TOP);
176 double coords[3];
177 // invert coordinates as curvi order is zyx
178 coords[0] = pointSourceLocation[p][2];
179 coords[1] = pointSourceLocation[p][1];
180 coords[2] = pointSourceLocation[p][0];
181
182 pointSourceLocation[p][_TOP] = this->interface->invertProjection(dir_top, coords);
183 }
184 }
185
187 Shortcuts s;
188 int numberOfData = s.SizeVariables + s.SizeParameters;
189 kernels::idx4 id_xyzf(basisSize, basisSize, basisSize, numberOfData);
190
191 dx[0] = 0;
192 dx[1] = 0;
193 dx[2] = 0;
194
195 for (int i = 0; i < basisSize; i++) {
196 for (int j = 0; j < basisSize; j++) {
197 dx[0] = std::max(
198 dx[0], luh[id_xyzf(i, j, basisSize - 1, s.curve_grid + 0)] - luh[id_xyzf(i, j, 0, s.curve_grid + 0)]
199 );
200 dx[1] = std::max(
201 dx[1], luh[id_xyzf(i, basisSize - 1, j, s.curve_grid + 1)] - luh[id_xyzf(i, 0, j, s.curve_grid + 1)]
202 );
203 dx[2] = std::max(
204 dx[2], luh[id_xyzf(basisSize - 1, i, j, s.curve_grid + 2)] - luh[id_xyzf(0, i, j, s.curve_grid + 2)]
205 );
206 }
207 }
208 }
209
210 void getElementCenter(const double* const luh, tarch::la::Vector<Dimensions, double>& center) {
211 int center_i1 = int(std::ceil(basisSize / 2.0));
212 int center_i2 = int(std::floor(basisSize / 2.0));
213
214 Shortcuts s;
215 int numberOfData = s.SizeVariables + s.SizeParameters;
216 kernels::idx4 id_xyzf(basisSize, basisSize, basisSize, numberOfData);
217
218 center[0] = (luh[id_xyzf(center_i1, center_i1, center_i1, s.curve_grid + 0)]
219 + luh[id_xyzf(center_i2, center_i2, center_i2, s.curve_grid + 0)])
220 / 2.0;
221 center[1] = (luh[id_xyzf(center_i1, center_i1, center_i1, s.curve_grid + 1)]
222 + luh[id_xyzf(center_i2, center_i2, center_i2, s.curve_grid + 1)])
223 / 2.0;
224 center[2] = (luh[id_xyzf(center_i1, center_i1, center_i1, s.curve_grid + 2)]
225 + luh[id_xyzf(center_i2, center_i2, center_i2, s.curve_grid + 2)])
226 / 2.0;
227 }
228
230
231 // Topography* topography;
232protected:
233 Curvi::Interface* interface;
234
237 std::string topography_string;
238};
239
240#endif
#define assertion5(expr, param0, param1, param2, param3, param4)
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
#define _TOP
virtual void initUnknownsPatch(double *luh, const tarch::la::Vector< Dimensions, double > &center, const tarch::la::Vector< Dimensions, double > &dx, double t, double dt)
void initPointSourceLocation(double pointSourceLocation[][3]) override
ContextCurvilinear(std::string &scenario_string, std::string &a_topography_string, DomainInformation *a_domain_info, SolverInformationADERDG< basisSize - 1 > *a_solver_info)
Curvi::Interface * interface
DomainInformation * getDomainInfo()
void getElementSize(const double *const luh, tarch::la::Vector< Dimensions, double > &dx)
void getElementCenter(const double *const luh, tarch::la::Vector< Dimensions, double > &center)
DomainInformation * domain_info
SolverInformationADERDG< basisSize - 1 > * solver_info
static void metricDerivatives(double dudx[][num_nodes], const double *const coordinates, const double *const dx, double *derivatives)
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