Peano
Loading...
Searching...
No Matches
Riemann.cpp
Go to the documentation of this file.
1#include "Riemann.h"
2
3#include "exahype2/VolumeIndex.h"
4#include "exahype2/enumerator/AoSLexicographicEnumerator.h"
5#include "exahype2/enumerator/FaceAoSLexicographicEnumerator.h"
6
7#include "peano4/utils/Loop.h"
8#include "peano4/datamanagement/FaceMarker.h"
9
10void exahype2::dg::interpolateRiemannSolution(
11 const peano4::datamanagement::FaceMarker& marker,
12 int order,
13 int unknowns,
14 const double* __restrict__ InterpolationMatrix1d,
15 const double* __restrict__ coarseGridFaceQ,
16 double* __restrict__ fineGridFaceQ
17) {
18 exahype2::enumerator::FaceAoSLexicographicEnumerator faceEnumerator(
19 marker.getSelectedFaceNumber(),
20 order+1,
21 1, // overlap
22 unknowns,
23 0 // auxiliaryVariables
24 );
25
26 #if Dimensions==2
27 const int totalNumberOfValues = (order+1) * 2 * unknowns;
28 const int InterpolationMatrixDimensions = (order+1); // no of rows/cols
29 #else
30 const int totalNumberOfValues = (order+1) * (order+1) * 2 * unknowns;
31 const int InterpolationMatrixDimensions = (order+1)*(order+1);
32 #endif
33 std::fill_n( fineGridFaceQ, totalNumberOfValues, 0.0 );
34
35 const int normal = marker.getSelectedFaceNumber()%Dimensions;
36 dfore(fineGridDoF,order+1,normal,0) {
37 dfore(coarseGridDoF,order+1,normal,0) {
38 double tensorProductWeight = 1.0;
39 for (int d=0; d<Dimensions; d++) {
40 if (d!=normal) {
41 int subMatrix = marker.getRelativePositionWithinFatherCell()(d);
42 assertion(subMatrix>=0);
43 assertion(subMatrix<=2);
44
45 int col = fineGridDoF(d);
46 int row = coarseGridDoF(d);
47 int linearisedMatrixIndex = subMatrix * InterpolationMatrixDimensions * InterpolationMatrixDimensions
48 + row * InterpolationMatrixDimensions
49 + col;
50 assertion(linearisedMatrixIndex>=0);
51 assertion2(linearisedMatrixIndex<3, linearisedMatrixIndex, InterpolationMatrixDimensions);
52 tensorProductWeight *= InterpolationMatrix1d[linearisedMatrixIndex];
53 }
54 }
55
56 tarch::la::Vector<Dimensions,int> shiftedFineGridDoF = fineGridDoF;
57 tarch::la::Vector<Dimensions,int> shiftedCoarseGridDoF = coarseGridDoF;
58 for (int leftRight=0; leftRight<2; leftRight++)
59 for (int unknown=0; unknown<unknowns; unknown++) {
60 shiftedFineGridDoF(normal) = leftRight;
61 shiftedCoarseGridDoF(normal) = leftRight;
62
63 fineGridFaceQ[ faceEnumerator(shiftedFineGridDoF,unknown) ] += tensorProductWeight * coarseGridFaceQ[ faceEnumerator(shiftedCoarseGridDoF,unknown) ];
64 }
65 }
66 }
67}
68
69void exahype2::dg::restrictAndAccumulateProjectedFacePolynomial(
70 const peano4::datamanagement::FaceMarker& marker,
71 int order,
72 int numberOfProjectedQuantities,
73 int unknowns,
74 int auxiliaryVariables,
75 const double* __restrict__ RestrictionMatrix1d,
76 const double* __restrict__ fineGridFaceQ,
77 double* __restrict__ coarseGridFaceQ
78) {
79 assertion(numberOfProjectedQuantities>=1);
80 // @todo Remove as it should work in principle for other values as well, but I haven't tested it
81 assertionEquals(numberOfProjectedQuantities,1);
82
83 exahype2::enumerator::FaceAoSLexicographicEnumerator faceEnumerator(
84 marker.getSelectedFaceNumber(),
85 order+1,
86 numberOfProjectedQuantities,
87 unknowns,
88 auxiliaryVariables
89 );
90
91 #if Dimensions==2
92 [[maybe_unused]] const int totalNumberOfValues = (order+1) * (numberOfProjectedQuantities*2) * (unknowns+auxiliaryVariables);
93 const int RestrictionMatrixDimensions = (order+1); // no of rows/cols
94 #else
95 [[maybe_unused]] const int totalNumberOfValues = (order+1) * (order+1) * (numberOfProjectedQuantities*2) * (unknowns+auxiliaryVariables);
96 const int RestrictionMatrixDimensions = (order+1)*(order+1);
97 #endif
98
99 const int normal = marker.getSelectedFaceNumber()%Dimensions;
100 const int leftRightFace = marker.getSelectedFaceNumber() < Dimensions ? 1 : 0;
101 dfore(fineGridDoF,order+1,normal,0) {
102 dfore(coarseGridDoF,order+1,normal,0) {
103 double tensorProductWeight = 1.0;
104 for (int d=0; d<Dimensions; d++) {
105 if (d!=normal) {
106 int subMatrix = marker.getRelativePositionWithinFatherCell()(d);
107 assertion(subMatrix>=0);
108 assertion(subMatrix<=2);
109
110 int row = fineGridDoF(d);
111 int col = coarseGridDoF(d);
112 int linearisedMatrixIndex = subMatrix * RestrictionMatrixDimensions * RestrictionMatrixDimensions
113 + row * RestrictionMatrixDimensions
114 + col;
115 assertion(linearisedMatrixIndex>=0);
116 assertion2(linearisedMatrixIndex<3, linearisedMatrixIndex, RestrictionMatrixDimensions);
117 tensorProductWeight *= RestrictionMatrix1d[linearisedMatrixIndex];
118 }
119 }
120
121 tarch::la::Vector<Dimensions,int> shiftedFineGridDoF = fineGridDoF;
122 tarch::la::Vector<Dimensions,int> shiftedCoarseGridDoF = coarseGridDoF;
123 for (int leftRight=leftRightFace * numberOfProjectedQuantities; leftRight<(leftRightFace+1)*numberOfProjectedQuantities; leftRight++)
124 for (int unknown=0; unknown<unknowns+auxiliaryVariables; unknown++) {
125 shiftedFineGridDoF(normal) = leftRight;
126 shiftedCoarseGridDoF(normal) = leftRight;
127
128 assertion( faceEnumerator(shiftedCoarseGridDoF,unknown)>=0 );
129 assertion( faceEnumerator(shiftedFineGridDoF,unknown)>=0 );
130
131 coarseGridFaceQ[ faceEnumerator(shiftedCoarseGridDoF,unknown) ] += tensorProductWeight * fineGridFaceQ[ faceEnumerator(shiftedFineGridDoF,unknown) ];
132 }
133 }
134 }
135}
136
137void exahype2::dg::clearSolutionProjection(
138 int order,
139 int unknowns,
140 int auxiliaryVariables,
141 int numberOfProjectedQuantities,
142 double* __restrict__ faceQ
143) {
144 #if Dimensions==2
145 const int numberOfDoFs = (order+1);
146 #else
147 const int numberOfDoFs = (order+1) * (order+1);
148 #endif
149 const int numberOfDoubles = numberOfDoFs * (unknowns + auxiliaryVariables) * 2 * numberOfProjectedQuantities;
150 std::fill_n( faceQ, numberOfDoubles, 0.0 );
151}
152
153void exahype2::dg::clearRiemannResult(
154 int order,
155 int unknowns,
156 double* __restrict__ faceQ
157) {
158 clearSolutionProjection(
159 order,
160 unknowns,
161 0, // auxiliaryVariables,
162 1, // numberOfProjectedQuantities
163 faceQ
164 );
165}
166
167void exahype2::dg::projectVolumetricDataOntoFaces(
168 [[maybe_unused]] const double* __restrict__ cellQ,
169 [[maybe_unused]] int order,
170 [[maybe_unused]] int unknowns,
171 [[maybe_unused]] int auxiliaryVariables,
172 [[maybe_unused]] const double* const __restrict__ BasisFunctionValuesLeft,
173 [[maybe_unused]] double* const __restrict__ faceQLeft,
174 [[maybe_unused]] double* const __restrict__ faceQRight,
175 [[maybe_unused]] double* const __restrict__ faceQBottom,
176 [[maybe_unused]] double* const __restrict__ faceQUp
177) {
178#if Dimensions==2
179 exahype2::enumerator::AoSLexicographicEnumerator cellEnumerator(1,order+1,0,unknowns,auxiliaryVariables);
180 exahype2::enumerator::FaceAoSLexicographicEnumerator leftRightFaceEnumerator(0,order+1,1,unknowns,auxiliaryVariables);
181 exahype2::enumerator::FaceAoSLexicographicEnumerator bottomTopFaceEnumerator(1,order+1,1,unknowns,auxiliaryVariables);
182
183 // Project left/right
184 for (int y=0; y<order+1; y++) {
185 const auto left = volumeIndex(1, y);
186 const auto right = volumeIndex(0, y);
187 for (int unknown=0; unknown<unknowns+auxiliaryVariables; unknown++) {
188 faceQLeft [ leftRightFaceEnumerator(left,unknown) ] = 0.0;
189 faceQRight[ leftRightFaceEnumerator(right,unknown) ] = 0.0;
190 }}
191
192 for (int y=0; y<order+1; y++) {
193 const auto left = volumeIndex(1, y);
194 const auto right = volumeIndex(0, y);
195 for (int x=0; x<order+1; x++) {
196 for (int unknown=0; unknown<unknowns+auxiliaryVariables; unknown++) {
197 faceQLeft [ leftRightFaceEnumerator(left,unknown) ] += BasisFunctionValuesLeft[x] * cellQ[ cellEnumerator(0,{x,y},unknown) ];
198 faceQRight[ leftRightFaceEnumerator(right,unknown) ] += BasisFunctionValuesLeft[order-x] * cellQ[ cellEnumerator(0,{x,y},unknown) ];
199 }}}
200
201 // Project up/down
202 for (int x=0; x<order+1; x++) {
203 const auto bottom = volumeIndex(x, 1);
204 const auto up = volumeIndex(x, 0);
205 for (int unknown=0; unknown<unknowns+auxiliaryVariables; unknown++) {
206 faceQBottom [ bottomTopFaceEnumerator(bottom,unknown) ] = 0.0;
207 faceQUp [ bottomTopFaceEnumerator(up,unknown) ] = 0.0;
208 }}
209
210 for (int x=0; x<order+1; x++) {
211 const auto bottom = volumeIndex(x, 1);
212 const auto up = volumeIndex(x, 0);
213 for (int y=0; y<order+1; y++) {
214 for (int unknown=0; unknown<unknowns+auxiliaryVariables; unknown++) {
215 faceQBottom [ bottomTopFaceEnumerator(bottom,unknown) ] += BasisFunctionValuesLeft[y] * cellQ[ cellEnumerator(0,{x,y},unknown) ];
216 faceQUp [ bottomTopFaceEnumerator(up,unknown) ] += BasisFunctionValuesLeft[order-y] * cellQ[ cellEnumerator(0,{x,y},unknown) ];
217 }}}
218#endif
219}
220
221void exahype2::dg::projectVolumetricDataOntoFaces(
222 [[maybe_unused]] const double* __restrict__ cellQ,
223 [[maybe_unused]] int order,
224 [[maybe_unused]] int unknowns,
225 [[maybe_unused]] int auxiliaryVariables,
226 [[maybe_unused]] const double* const __restrict__ BasisFunctionValuesLeft,
227 [[maybe_unused]] double* const __restrict__ faceQLeft,
228 [[maybe_unused]] double* const __restrict__ faceQRight,
229 [[maybe_unused]] double* const __restrict__ faceQBottom,
230 [[maybe_unused]] double* const __restrict__ faceQUp,
231 [[maybe_unused]] double* const __restrict__ faceQFront,
232 [[maybe_unused]] double* const __restrict__ faceQBack
233) {
234#if Dimensions==3
235 exahype2::enumerator::AoSLexicographicEnumerator cellEnumerator(1,order+1,0,unknowns,auxiliaryVariables);
236 exahype2::enumerator::FaceAoSLexicographicEnumerator leftRightFaceEnumerator(0,order+1,1,unknowns,auxiliaryVariables);
237 exahype2::enumerator::FaceAoSLexicographicEnumerator bottomTopFaceEnumerator(1,order+1,1,unknowns,auxiliaryVariables);
238 exahype2::enumerator::FaceAoSLexicographicEnumerator frontBackFaceEnumerator(2,order+1,1,unknowns,auxiliaryVariables);
239
240 // Project left/right
241 for (int y=0; y<order+1; y++) {
242 for (int z=0; z<order+1; z++) {
243 const auto left = volumeIndex(1, y, z);
244 const auto right = volumeIndex(0, y, z);
245 for (int unknown=0; unknown<unknowns+auxiliaryVariables; unknown++) {
246 faceQLeft [ leftRightFaceEnumerator(left,unknown) ] = 0.0;
247 faceQRight [ leftRightFaceEnumerator(right,unknown) ] = 0.0;
248 }}}
249
250 for (int y=0; y<order+1; y++) {
251 for (int z=0; z<order+1; z++) {
252 const auto left = volumeIndex(1, y, z);
253 const auto right = volumeIndex(0, y, z);
254 for (int x=0; x<order+1; x++) {
255 for (int unknown=0; unknown<unknowns+auxiliaryVariables; unknown++) {
256 faceQLeft [ leftRightFaceEnumerator(left,unknown) ] += BasisFunctionValuesLeft[x] * cellQ[ cellEnumerator(0,{x,y,z},unknown) ];
257 faceQRight [ leftRightFaceEnumerator(right,unknown) ] += BasisFunctionValuesLeft[order-x] * cellQ[ cellEnumerator(0,{x,y,z},unknown) ];
258 }}}}
259
260 // Project down/up
261 for (int x=0; x<order+1; x++) {
262 for (int z=0; z<order+1; z++) {
263 const auto bottom = volumeIndex(x, 1, z);
264 const auto up = volumeIndex(x, 0, z);
265 for (int unknown=0; unknown<unknowns+auxiliaryVariables; unknown++) {
266 faceQBottom [ bottomTopFaceEnumerator(bottom,unknown) ] = 0.0;
267 faceQUp [ bottomTopFaceEnumerator(up,unknown) ] = 0.0;
268 }}}
269
270 for (int x=0; x<order+1; x++) {
271 for (int z=0; z<order+1; z++) {
272 const auto bottom = volumeIndex(x, 1, z);
273 const auto up = volumeIndex(x, 0, z);
274 for (int y=0; y<order+1; y++) {
275 for (int unknown=0; unknown<unknowns+auxiliaryVariables; unknown++) {
276 faceQBottom [ bottomTopFaceEnumerator(bottom,unknown) ] += BasisFunctionValuesLeft[y] * cellQ[ cellEnumerator(0,{x,y,z},unknown) ];
277 faceQUp [ bottomTopFaceEnumerator(up,unknown) ] += BasisFunctionValuesLeft[order-y] * cellQ[ cellEnumerator(0,{x,y,z},unknown) ];
278 }}}}
279
280 // Project front/back
281 for (int x=0; x<order+1; x++) {
282 for (int y=0; y<order+1; y++) {
283 const auto front = volumeIndex(x, y, 1);
284 const auto back = volumeIndex(x, y, 0);
285 for (int unknown=0; unknown<unknowns+auxiliaryVariables; unknown++) {
286 faceQFront [ frontBackFaceEnumerator(front,unknown) ] = 0.0;
287 faceQBack [ frontBackFaceEnumerator(back,unknown) ] = 0.0;
288 }}}
289
290 for (int x=0; x<order+1; x++) {
291 for (int y=0; y<order+1; y++) {
292 const auto front = volumeIndex(x, y, 1);
293 const auto back = volumeIndex(x, y, 0);
294 for (int z=0; z<order+1; z++) {
295 for (int unknown=0; unknown<unknowns+auxiliaryVariables; unknown++) {
296 faceQFront [ frontBackFaceEnumerator(front,unknown) ] += BasisFunctionValuesLeft[z] * cellQ[ cellEnumerator(0,{x,y,z},unknown) ];
297 faceQBack [ frontBackFaceEnumerator(back,unknown) ] += BasisFunctionValuesLeft[order-z] * cellQ[ cellEnumerator(0,{x,y,z},unknown) ];
298 }}}}
299#endif
300}
301
302void exahype2::dg::projectVolumetricDataAndGradientOntoFaces(
303 [[maybe_unused]] const double* __restrict__ cellQ,
304 [[maybe_unused]] int order,
305 [[maybe_unused]] int unknowns,
306 [[maybe_unused]] int auxiliaryVariables,
307 [[maybe_unused]] const double* const __restrict__ BasisFunctionValuesLeft,
308 [[maybe_unused]] double* const __restrict__ faceQLeft,
309 [[maybe_unused]] double* const __restrict__ faceQRight,
310 [[maybe_unused]] double* const __restrict__ faceQBottom,
311 [[maybe_unused]] double* const __restrict__ faceQUp
312) {
313 assertionMsg(false, "not implemented yet");
314}
315
316void exahype2::dg::projectVolumetricDataAndGradientOntoFaces(
317 [[maybe_unused]] const double* __restrict__ cellQ,
318 [[maybe_unused]] int order,
319 [[maybe_unused]] int unknowns,
320 [[maybe_unused]] int auxiliaryVariables,
321 [[maybe_unused]] const double* const __restrict__ BasisFunctionValuesLeft,
322 [[maybe_unused]] double* const __restrict__ faceQLeft,
323 [[maybe_unused]] double* const __restrict__ faceQRight,
324 [[maybe_unused]] double* const __restrict__ faceQBottom,
325 [[maybe_unused]] double* const __restrict__ faceQUp,
326 [[maybe_unused]] double* const __restrict__ faceQFront,
327 [[maybe_unused]] double* const __restrict__ faceQBack
328) {
329 assertionMsg(false, "not implemented yet");
330}
331
332void exahype2::dg::integrateOverRiemannSolutionsAndAddToVolume_GaussLegendre(
333 [[maybe_unused]] const double* const __restrict__ faceQLeft,
334 [[maybe_unused]] const double* const __restrict__ faceQRight,
335 [[maybe_unused]] const double* const __restrict__ faceQBottom,
336 [[maybe_unused]] const double* const __restrict__ faceQUp,
337 [[maybe_unused]] int order,
338 [[maybe_unused]] int unknowns,
339 [[maybe_unused]] const int auxiliaryVariables,
340 [[maybe_unused]] const tarch::la::Vector<2,double>& cellSize,
341 [[maybe_unused]] const double* const __restrict__ BasisFunctionValuesLeft,
342 [[maybe_unused]] const double* __restrict__ MassMatrixDiagonal1d,
343 [[maybe_unused]] double* __restrict__ cellQ
344) {
345#if Dimensions==2
346 [[maybe_unused]] const int strideQ = unknowns + auxiliaryVariables;
347
348 exahype2::enumerator::AoSLexicographicEnumerator QEnumerator(
349 1, // number of cells in one memory chunk
350 order+1, // numberOfDoFsPerAxisInCell
351 0, // no halo
352 unknowns,
353 0 // no auxiliary variables here
354 );
355
356 exahype2::enumerator::FaceAoSLexicographicEnumerator leftFaceEnumerator (0, order+1, 1, unknowns, 0);
357 exahype2::enumerator::FaceAoSLexicographicEnumerator rightFaceEnumerator (0+Dimensions, order+1, 1, unknowns, 0);
358 exahype2::enumerator::FaceAoSLexicographicEnumerator bottomFaceEnumerator(1, order+1, 1, unknowns, 0);
359 exahype2::enumerator::FaceAoSLexicographicEnumerator upperFaceEnumerator (1+Dimensions, order+1, 1, unknowns, 0);
360
361 dfor( dof, order+1 ) {
362 const double coeffLeft = MassMatrixDiagonal1d[dof[1]]*BasisFunctionValuesLeft[dof[0]] * cellSize(1);
363 const double coeffRight = MassMatrixDiagonal1d[dof[1]]*BasisFunctionValuesLeft[order-dof[0]] * cellSize(1);
364
365 const double coeffBottom = MassMatrixDiagonal1d[dof[0]]*BasisFunctionValuesLeft[dof[1]] * cellSize(0);
366 const double coeffUp = MassMatrixDiagonal1d[dof[0]]*BasisFunctionValuesLeft[order-dof[1]] * cellSize(0);
367
368 tarch::la::Vector<Dimensions,int> leftDoF = dof;
369 tarch::la::Vector<Dimensions,int> rightDoF = dof;
370 tarch::la::Vector<Dimensions,int> bottomDoF = dof;
371 tarch::la::Vector<Dimensions,int> upperDoF = dof;
372
373 leftDoF(0) = 1;
374 rightDoF(0) = 0;
375 bottomDoF(1) = 1;
376 upperDoF(1) = 0;
377
378 for(int var=0; var<unknowns; var++) {
379 cellQ[ QEnumerator(0,dof,var) ] +=
380 coeffLeft * faceQLeft [leftFaceEnumerator( leftDoF, var)]
381 - coeffRight * faceQRight [rightFaceEnumerator( rightDoF, var)]
382 + coeffBottom * faceQBottom [bottomFaceEnumerator(bottomDoF,var)]
383 - coeffUp * faceQUp [upperFaceEnumerator( upperDoF, var)];
384 }
385 }
386#endif
387}
388
389void exahype2::dg::integrateOverRiemannSolutionsAndAddToVolume_GaussLegendre(
390 [[maybe_unused]] const double* const __restrict__ faceQLeft,
391 [[maybe_unused]] const double* const __restrict__ faceQRight,
392 [[maybe_unused]] const double* const __restrict__ faceQBottom,
393 [[maybe_unused]] const double* const __restrict__ faceQUp,
394 [[maybe_unused]] const double* const __restrict__ faceQFront,
395 [[maybe_unused]] const double* const __restrict__ faceQBack,
396 [[maybe_unused]] int order,
397 [[maybe_unused]] int unknowns,
398 [[maybe_unused]] const int auxiliaryVariables,
399 [[maybe_unused]] const tarch::la::Vector<3,double>& cellSize,
400 [[maybe_unused]] const double* const __restrict__ BasisFunctionValuesLeft,
401 [[maybe_unused]] const double* __restrict__ MassMatrixDiagonal1d,
402 [[maybe_unused]] double* __restrict__ cellQ
403) {
404#if Dimensions==3
405 exahype2::enumerator::AoSLexicographicEnumerator QEnumerator(
406 1, // number of cells in one memory chunk
407 order+1, // numberOfDoFsPerAxisInCell
408 0, // no halo
409 unknowns,
410 0 // no auxiliary variables here
411 );
412
413 exahype2::enumerator::FaceAoSLexicographicEnumerator leftFaceEnumerator (0, order+1, 1, unknowns, 0);
414 exahype2::enumerator::FaceAoSLexicographicEnumerator rightFaceEnumerator (0+Dimensions, order+1, 1, unknowns, 0);
415 exahype2::enumerator::FaceAoSLexicographicEnumerator bottomFaceEnumerator(1, order+1, 1, unknowns, 0);
416 exahype2::enumerator::FaceAoSLexicographicEnumerator upperFaceEnumerator( 1+Dimensions, order+1, 1, unknowns, 0);
417 exahype2::enumerator::FaceAoSLexicographicEnumerator frontFaceEnumerator( 2, order+1, 1, unknowns, 0);
418 exahype2::enumerator::FaceAoSLexicographicEnumerator backFaceEnumerator( 2+Dimensions, order+1, 1, unknowns, 0);
419
420 dfor( dof, order+1 ) {
421 //coefficients of faces in x direction
422 const double coeffLeft = BasisFunctionValuesLeft[dof[0]] * MassMatrixDiagonal1d[dof[1]] * MassMatrixDiagonal1d[dof[2]] * cellSize(1) * cellSize(2);
423 const double coeffRight = BasisFunctionValuesLeft[order-dof[0]] * MassMatrixDiagonal1d[dof[1]] * MassMatrixDiagonal1d[dof[2]] * cellSize(1) * cellSize(2);
424 //coefficients of faces in y direction
425 const double coeffBottom = MassMatrixDiagonal1d[dof[0]] * BasisFunctionValuesLeft[dof[1]] * MassMatrixDiagonal1d[dof[2]] * cellSize(0) * cellSize(2);
426 const double coeffUp = MassMatrixDiagonal1d[dof[0]] * BasisFunctionValuesLeft[order-dof[1]] * MassMatrixDiagonal1d[dof[2]] * cellSize(0) * cellSize(2);
427 //coefficients of faces in z direction
428 const double coeffFront = MassMatrixDiagonal1d[dof[0]] * MassMatrixDiagonal1d[dof[1]] * BasisFunctionValuesLeft[dof[2]] * cellSize(0) * cellSize(1);
429 const double coeffBack = MassMatrixDiagonal1d[dof[0]] * MassMatrixDiagonal1d[dof[1]] * BasisFunctionValuesLeft[order-dof[2]] * cellSize(0) * cellSize(1);
430
431 tarch::la::Vector<Dimensions,int> leftDoF = dof;
432 tarch::la::Vector<Dimensions,int> rightDoF = dof;
433 tarch::la::Vector<Dimensions,int> frontDoF = dof;
434 tarch::la::Vector<Dimensions,int> backDoF = dof;
435 tarch::la::Vector<Dimensions,int> bottomDoF = dof;
436 tarch::la::Vector<Dimensions,int> upperDoF = dof;
437
438 leftDoF(0) = 1;
439 rightDoF(0) = 0;
440 bottomDoF(1) = 1;
441 upperDoF(1) = 0;
442 frontDoF(2) = 1;
443 backDoF(2) = 0;
444
445 for(int var=0; var<unknowns; var++) {
446 cellQ[ QEnumerator(0,dof,var) ] +=
447 + coeffLeft * faceQLeft [leftFaceEnumerator (leftDoF, var)]
448 - coeffRight * faceQRight [rightFaceEnumerator (rightDoF, var)]
449 + coeffBottom * faceQBottom [bottomFaceEnumerator(bottomDoF,var)]
450 - coeffUp * faceQUp [upperFaceEnumerator (upperDoF, var)]
451 + coeffFront * faceQFront [frontFaceEnumerator (frontDoF, var)]
452 - coeffBack * faceQBack [backFaceEnumerator (backDoF, var)];
453 }
454 }
455#endif
456}
assertion2
#define assertion2(expr, param0, param1)
Definition Assertions.h:1019
assertionEquals
#define assertionEquals(lhs, rhs)
Definition Assertions.h:1034
assertionMsg
#define assertionMsg(expr, message)
Definition Assertions.h:1076
assertion
#define assertion(expr)
Definition Assertions.h:1017
dfore
#define dfore(counter, max, dim, value)
This is an exclusive d-dimensional for loop.
Definition Loop.h:942
dfor
#define dfor(counter, max)
d-dimensional Loop
Definition Loop.h:313
exahype2::dg::restrictAndAccumulateProjectedFacePolynomial
void restrictAndAccumulateProjectedFacePolynomial(const peano4::datamanagement::FaceMarker &marker, int order, int numberOfProjectedQuantities, int unknowns, int auxiliaryVariables, const double *__restrict__ RestrictionMatrix1d, const double *__restrict__ fineGridFaceQ, double *__restrict__ coarseGridFaceQ)
Counterpart of interpolateRiemannSolution().
Definition Riemann.cpp:69
exahype2::dg::clearRiemannResult
void clearRiemannResult(int order, int unknowns, double *__restrict__ faceQ)
Definition Riemann.cpp:153
exahype2::dg::clearSolutionProjection
void clearSolutionProjection(int order, int unknowns, int auxiliaryVariables, int numberOfProjectedQuantities, double *__restrict__ faceQ)
Set the whole solution projection on the face from left and right to zero.
Definition Riemann.cpp:137
exahype2::dg::projectVolumetricDataAndGradientOntoFaces
void projectVolumetricDataAndGradientOntoFaces(const double *__restrict__ cellQ, int order, int unknowns, int auxiliaryVariables, const double *const __restrict__ BasisFunctionValuesLeft, double *const __restrict__ faceQLeft, double *const __restrict__ faceQRight, double *const __restrict__ faceQBottom, double *const __restrict__ faceQUp)
Definition Riemann.cpp:302
exahype2::dg::projectVolumetricDataOntoFaces
void projectVolumetricDataOntoFaces(const double *__restrict__ cellQ, int order, int unknowns, int auxiliaryVariables, const double *const __restrict__ BasisFunctionValuesLeft, double *const __restrict__ faceQLeft, double *const __restrict__ faceQRight, double *const __restrict__ faceQBottom, double *const __restrict__ faceQUp)
Take polynomial within cell and project it onto the faces.
Definition Riemann.cpp:167
exahype2::dg::interpolateRiemannSolution
void interpolateRiemannSolution(const peano4::datamanagement::FaceMarker &marker, int order, int unknowns, const double *__restrict__ InterpolationMatrix1d, const double *__restrict__ coarseGridFaceQ, double *__restrict__ fineGridFaceQ)
Definition Riemann.cpp:10
exahype2::dg::integrateOverRiemannSolutionsAndAddToVolume_GaussLegendre
void integrateOverRiemannSolutionsAndAddToVolume_GaussLegendre(const double *const __restrict__ faceQLeft, const double *const __restrict__ faceQRight, const double *const __restrict__ faceQBottom, const double *const __restrict__ faceQUp, int order, int unknowns, const int auxiliaryVariables, const tarch::la::Vector< 2, double > &cellSize, const double *const __restrict__ BasisFunctionValuesLeft, const double *__restrict__ MassMatrixDiagonal1d, double *__restrict__ cellQ)
Given a numerical flux at the various faces, this computes and adds the Riemann integral of this flux...
Definition Riemann.cpp:332
exahype2::volumeIndex
auto volumeIndex(Args... args)
Definition VolumeIndex.h:54
peano4::datamanagement::FaceMarker
Provide information about selected face.
Definition FaceMarker.h:35
tarch::la::Vector
Simple vector class.
Definition Vector.h:134