Peano
Loading...
Searching...
No Matches
Interpolation.cpp
Go to the documentation of this file.
1#include "Interpolation.h"
2
3#include <vector>
4
5#include "Enumeration.h"
6#include "peano4/utils/Loop.h"
13
14namespace {
15 [[maybe_unused]] tarch::logging::Log _log("toolbox::blockstructured");
16
17 tarch::multicore::BooleanSemaphore interpolationMapSemaphore;
18} // namespace
19
21 createPiecewiseConstantInterpolationMatrix(int numberOfDoFsPerAxisInPatch) {
22#if Dimensions == 2
23 const int patchSize = numberOfDoFsPerAxisInPatch * numberOfDoFsPerAxisInPatch;
24 const int numberOfPatches = 3 * 3;
25#else
26 const int patchSize = numberOfDoFsPerAxisInPatch * numberOfDoFsPerAxisInPatch
27 * numberOfDoFsPerAxisInPatch;
28 const int numberOfPatches = 3 * 3 * 3;
29#endif
30
32 patchSize * numberOfPatches,
33 patchSize
34 );
35
36 int currentRow = 0;
37 dfor3(patchIndex) dfor(volumeIndex, numberOfDoFsPerAxisInPatch) {
39 sourceCell = (volumeIndex + patchIndex * numberOfDoFsPerAxisInPatch) / 3;
40 int sourceCellLinearised = peano4::utils::dLinearised(
41 sourceCell,
42 numberOfDoFsPerAxisInPatch
43 );
44 (*result)(currentRow, sourceCellLinearised) = 1.0;
45 currentRow++;
46 }
48
50 "createPiecewiseConstantInterpolationMatrix(int)",
51 "created new volumetric matrix for "
52 << numberOfDoFsPerAxisInPatch << ": " << result->toString()
53 );
54 return result;
55}
56
59 [[maybe_unused]] int numberOfDoFsPerAxisInPatch,
60 [[maybe_unused]] int normal,
61 [[maybe_unused]] int overlap
62 ) {
63 tarch::la::DynamicMatrix P1d(3, 1, {{1.0}, {1.0}, {1.0}});
64 P1d.replicateRows(3, numberOfDoFsPerAxisInPatch, 1, true);
65
67 P1d,
68 numberOfDoFsPerAxisInPatch,
69 normal
70 );
71}
72
75 const tarch::la::DynamicMatrix& P1d,
76 int numberOfDoFsPerAxisInPatch,
77 int normal
78 ) {
79 assertionEquals1(P1d.cols(), numberOfDoFsPerAxisInPatch, P1d.toString());
80 assertionEquals1(P1d.rows(), 3 * numberOfDoFsPerAxisInPatch, P1d.toString());
81
82#if Dimensions == 3
84#else
86#endif
87
88 int pattern = 0;
89 switch (normal) {
90 case 0:
91 pattern = 1;
92 break;
93 case 1:
94 pattern = numberOfDoFsPerAxisInPatch;
95 break;
96 case 2:
97 pattern = numberOfDoFsPerAxisInPatch * numberOfDoFsPerAxisInPatch;
98 break;
99 }
100
101 P->insertEmptyColumns(pattern, pattern, pattern);
102 P->replicateRows(pattern, 2, pattern, false);
103 logDebug(
104 "createInterpolationMatrixFrom1dTemplateByInsertingZeroColsAndRows(...)",
105 "matrix for normal=" << normal << ": " << P->toString()
106 );
107 return P;
108}
109
112 const tarch::la::DynamicMatrix& P1d,
113 int numberOfDoFsPerAxisInPatch,
114 int normal
115 ) {
116 assertionEquals1(P1d.cols(), numberOfDoFsPerAxisInPatch, P1d.toString());
117 assertionEquals1(P1d.rows(), 3 * numberOfDoFsPerAxisInPatch, P1d.toString());
118
119#if Dimensions == 3
121#else
123#endif
124
125 int pattern = 0;
126 switch (normal) {
127 case 0:
128 pattern = 1;
129 break;
130 case 1:
131 pattern = numberOfDoFsPerAxisInPatch;
132 break;
133 case 2:
134 pattern = numberOfDoFsPerAxisInPatch * numberOfDoFsPerAxisInPatch;
135 break;
136 }
137
138 P->scale(0.5);
139
140 // Split matrix into blocks of pattern columns and copy each block once.
141 // Matrix grows by factor of two which reflects the fact that the overlap
142 // is of size 1.
143 P->replicateCols(pattern, 2, 0, false);
144 P->replicateRows(pattern, 2, 0, false);
145 logDebug(
146 "createInterpolationMatrixFrom1dTemplateByLinearInterpolationAlongNormal(...)",
147 "matrix for normal=" << normal << ": " << P->toString()
148 );
149 return P;
150}
151
154 int numberOfDoFsPerAxisInPatch,
155 int normal,
156 bool extrapolateLinearly,
157 bool interpolateLinearlyBetweenRealAndRestrictedVolumes
158 ) {
160 3,
161 3,
162 {{1.0 / 3.0, 2.0 / 3.0, 0.0},
163 {0.0, 3.0 / 3.0, 0.0},
164 {0.0, 2.0 / 3.0, 1.0 / 3.0}}
165 );
166 P1d.replicateRows(3, numberOfDoFsPerAxisInPatch, 1, true);
167 P1d.removeColumn(0);
168 P1d.removeColumn(numberOfDoFsPerAxisInPatch);
169
170 if (extrapolateLinearly) {
171 P1d(0, 0) = 2.0 / 3.0 + 1.0 / 3.0 * 2.0;
172 P1d(0, 1) = -1.0 / 3.0;
173 P1d(
174 numberOfDoFsPerAxisInPatch * 3 - 1,
175 numberOfDoFsPerAxisInPatch - 1
176 ) = 2.0 / 3.0 + 1.0 / 3.0 * 2.0;
177 P1d(
178 numberOfDoFsPerAxisInPatch * 3 - 1,
179 numberOfDoFsPerAxisInPatch - 2
180 ) = -1.0 / 3.0;
181 } else {
182 P1d(0, 0) = 1.0;
183 P1d(0, 1) = 0.0;
184 P1d(
185 numberOfDoFsPerAxisInPatch * 3 - 1,
186 numberOfDoFsPerAxisInPatch - 1
187 ) = 1.0;
188 P1d(
189 numberOfDoFsPerAxisInPatch * 3 - 1,
190 numberOfDoFsPerAxisInPatch - 2
191 ) = 0.0;
192 }
193
194 logDebug(
195 "createLinearInterpolationMatrix(...)",
196 "1d matrix: " << P1d.toString()
197 );
198
199 if (interpolateLinearlyBetweenRealAndRestrictedVolumes) {
201 P1d,
202 numberOfDoFsPerAxisInPatch,
203 normal
204 );
205 } else {
207 P1d,
208 numberOfDoFsPerAxisInPatch,
209 normal
210 );
211 }
212}
213
216 int numberOfDoFsPerAxisInPatch,
217 bool extrapolateLinearly
218 ) {
220 "createLinearInterpolationMatrix(...)",
221 numberOfDoFsPerAxisInPatch
222 );
224 3,
225 3,
226 {{1.0 / 3.0, 2.0 / 3.0, 0.0},
227 {0.0, 3.0 / 3.0, 0.0},
228 {0.0, 2.0 / 3.0, 1.0 / 3.0}}
229 );
230 P1d.replicateRows(3, numberOfDoFsPerAxisInPatch, 1, true);
231 P1d.removeColumn(0);
232 P1d.removeColumn(numberOfDoFsPerAxisInPatch);
233
234 if (extrapolateLinearly) {
235 P1d(0, 0) = 2.0 / 3.0 + 1.0 / 3.0 * 2.0;
236 P1d(0, 1) = -1.0 / 3.0;
237 P1d(
238 numberOfDoFsPerAxisInPatch * 3 - 1,
239 numberOfDoFsPerAxisInPatch - 1
240 ) = 2.0 / 3.0 + 1.0 / 3.0 * 2.0;
241 P1d(
242 numberOfDoFsPerAxisInPatch * 3 - 1,
243 numberOfDoFsPerAxisInPatch - 2
244 ) = -1.0 / 3.0;
245 } else {
246 P1d(0, 0) = 1.0;
247 P1d(0, 1) = 0.0;
248 P1d(
249 numberOfDoFsPerAxisInPatch * 3 - 1,
250 numberOfDoFsPerAxisInPatch - 1
251 ) = 1.0;
252 P1d(
253 numberOfDoFsPerAxisInPatch * 3 - 1,
254 numberOfDoFsPerAxisInPatch - 2
255 ) = 0.0;
256 }
257
258 logDebug(
259 "createLinearInterpolationMatrix(...)",
260 "1d interpolation matrix: " << P1d.toString(true)
261 );
262
263 tarch::la::DynamicMatrix P2d(P1d, P1d, false);
264 logDebug(
265 "createLinearInterpolationMatrix(...)",
266 "2d interpolation matrix: " << P2d.toString(true)
267 );
268
269 logTraceOut("createLinearInterpolationMatrix(...)");
270#if Dimensions == 3
271 return new tarch::la::DynamicMatrix(P2d, P1d, false);
272#else
273 return new tarch::la::DynamicMatrix(P2d);
274#endif
275}
276
280 int numberOfDoFsPerAxisInPatch,
281 int overlap,
282 int unknowns,
283 const double* __restrict__ fineGridCellValuesLeft,
284 const double* __restrict__ fineGridCellValuesRight,
285 double* fineGridFaceValues
286 ) {
288 "projectInterpolatedFineCellsOnHaloLayer_AoS(...)",
289 marker.toString(),
290 numberOfDoFsPerAxisInPatch,
291 overlap
292 );
293
294 const int normal = marker.getSelectedFaceNumber() % Dimensions;
295
296 dfore(kFine, numberOfDoFsPerAxisInPatch, normal, 0) {
297 for (int iFine = 0; iFine < overlap; iFine++) {
298 tarch::la::Vector<Dimensions, int> fineVolumeLeftDest = kFine;
299 tarch::la::Vector<Dimensions, int> fineVolumeLeftSrc = kFine;
300
301 fineVolumeLeftDest(normal) = iFine;
302 fineVolumeLeftSrc(normal) = numberOfDoFsPerAxisInPatch - overlap + iFine;
303
304 tarch::la::Vector<Dimensions, int> fineVolumeRightDest = kFine;
305 tarch::la::Vector<Dimensions, int> fineVolumeRightSrc = kFine;
306
307 fineVolumeRightDest(normal) = iFine + overlap;
308 fineVolumeRightSrc(normal) = iFine;
309
310 int fineVolumeLeftDestLinearised = serialiseVoxelIndexInOverlap(
311 fineVolumeLeftDest,
312 numberOfDoFsPerAxisInPatch,
313 overlap,
314 normal
315 );
316 int fineVolumeRightDestLinearised = serialiseVoxelIndexInOverlap(
317 fineVolumeRightDest,
318 numberOfDoFsPerAxisInPatch,
319 overlap,
320 normal
321 );
322 int fineVolumeLeftSrcLinearised = peano4::utils::dLinearised(
323 fineVolumeLeftSrc,
324 numberOfDoFsPerAxisInPatch
325 );
326 int fineVolumeRightSrcLinearised = peano4::utils::dLinearised(
327 fineVolumeRightSrc,
328 numberOfDoFsPerAxisInPatch
329 );
330
331 for (int j = 0; j < unknowns; j++) {
332 fineGridFaceValues[fineVolumeLeftDestLinearised * unknowns + j]
333 = fineGridCellValuesLeft[fineVolumeLeftSrcLinearised * unknowns + j];
334 fineGridFaceValues[fineVolumeRightDestLinearised * unknowns + j]
335 = fineGridCellValuesRight
336 [fineVolumeRightSrcLinearised * unknowns + j];
337 }
338 }
339 }
340
341 logTraceOut("projectInterpolatedFineCellsOnHaloLayer_AoS(...)");
342}
343
344//
345// ==========================================
346// Piece-wise constant interpolation routines
347// ==========================================
348//
351 int numberOfDoFsPerAxisInPatch,
352 int overlap,
353 int unknowns,
354 const double* __restrict__ coarseGridFaceValues,
355 double* __restrict__ fineGridFaceValues
356) {
358 "interpolateHaloLayer_AoS_piecewise_constant(...)",
359 marker.toString(),
360 numberOfDoFsPerAxisInPatch,
361 overlap,
362 unknowns
363 );
364
365 const int normal = marker.getSelectedFaceNumber() % Dimensions;
366
367 static internal::FaceInterpolationMap faceInterpolationMap;
369 key(numberOfDoFsPerAxisInPatch, overlap, normal);
370
371 tarch::multicore::Lock lock(interpolationMapSemaphore);
372 if (faceInterpolationMap.count(key) == 0) {
373 faceInterpolationMap[key] = internal::
374 createPiecewiseConstantInterpolationMatrix(
375 numberOfDoFsPerAxisInPatch,
376 normal,
377 overlap
378 );
379 }
380 tarch::la::DynamicMatrix* P = faceInterpolationMap[key];
381 lock.free();
382
383#if Dimensions == 2
384 const int patchSize = numberOfDoFsPerAxisInPatch * 2 * overlap;
385#else
386 const int patchSize = numberOfDoFsPerAxisInPatch * numberOfDoFsPerAxisInPatch
387 * 2 * overlap;
388#endif
389
390 logDebug(
391 "interpolateHaloLayer_AoS_piecewise_constant(...)",
392 marker.toString() << ": " << P->toString()
393 );
394
395 P->batchedMultiplyAoS(
396 fineGridFaceValues, // image
397 coarseGridFaceValues, // preimage
398 unknowns, // batch size, i.e. how often to apply it in one AoS rush
399 patchSize, // result size, i.e. size of image
401 marker.getRelativePositionWithinFatherCell(),
402 normal
403 ) * patchSize
404 );
405
406 logTraceOut("interpolateHaloLayer_AoS_piecewise_constant(...)");
407}
408
411 int numberOfDoFsPerAxisInPatch,
412 int overlap,
413 int unknowns,
414 const double* __restrict__ coarseGridCellValues,
415 const double* __restrict__ coarseGridFaceValues,
416 double* __restrict__ fineGridFaceValues
417) {
419 "interpolateHaloLayer_AoS_piecewise_constant(...)",
420 marker.toString(),
421 numberOfDoFsPerAxisInPatch,
422 overlap,
423 unknowns
424 );
425
426 if (marker.isInteriorFaceWithinPatch()) {
427#if Dimensions == 2
428 const int patchSize = numberOfDoFsPerAxisInPatch
429 * numberOfDoFsPerAxisInPatch;
430#else
431 const int patchSize = numberOfDoFsPerAxisInPatch * numberOfDoFsPerAxisInPatch
432 * numberOfDoFsPerAxisInPatch;
433#endif
434
435 double* leftAdjacentPatchData = tarch::allocateMemory<double>(
436 patchSize * unknowns,
438 );
439 double* rightAdjacentPatchData = tarch::allocateMemory<double>(
440 patchSize * unknowns,
442 );
443
445 leftAdjacentPatchIndex = marker.getRelativePositionWithinFatherCell();
447 rightAdjacentPatchIndex = marker.getRelativePositionWithinFatherCell();
448 leftAdjacentPatchIndex(marker.getSelectedFaceNumber() % Dimensions)--;
449
450 internal::interpolateCell_AoS_piecewise_constant(
451 leftAdjacentPatchIndex,
452 numberOfDoFsPerAxisInPatch,
453 unknowns,
454 coarseGridCellValues,
455 leftAdjacentPatchData
456 );
457 internal::interpolateCell_AoS_piecewise_constant(
458 rightAdjacentPatchIndex,
459 numberOfDoFsPerAxisInPatch,
460 unknowns,
461 coarseGridCellValues,
462 rightAdjacentPatchData
463 );
464
465 internal::projectInterpolatedFineCellsOnHaloLayer_AoS(
466 marker,
467 numberOfDoFsPerAxisInPatch,
468 overlap,
469 unknowns,
470 leftAdjacentPatchData,
471 rightAdjacentPatchData,
472 fineGridFaceValues
473 );
474
475 tarch::freeMemory(leftAdjacentPatchData, tarch::MemoryLocation::Heap);
476 tarch::freeMemory(rightAdjacentPatchData, tarch::MemoryLocation::Heap);
477 } else {
479 marker,
480 numberOfDoFsPerAxisInPatch,
481 overlap,
482 unknowns,
483 coarseGridFaceValues,
484 fineGridFaceValues
485 );
486 }
487
488 logTraceOut("interpolateHaloLayer_AoS_piecewise_constant(...)");
489}
490
493 int numberOfDoFsPerAxisInPatch,
494 int unknowns,
495 const double* __restrict__ coarseGridCellValues,
496 double* __restrict__ fineGridCellValues
497) {
499 "interpolateCell_AoS_piecewise_constant(...)",
500 marker.toString(),
501 numberOfDoFsPerAxisInPatch,
502 unknowns
503 );
504
505 internal::interpolateCell_AoS_piecewise_constant(
506 marker.getRelativePositionWithinFatherCell(),
507 numberOfDoFsPerAxisInPatch,
508 unknowns,
509 coarseGridCellValues,
510 fineGridCellValues
511 );
512
513 logTraceOut("interpolateCell_AoS_piecewise_constant(...)");
514}
515
517 const tarch::la::Vector<Dimensions, int>& relativePositionWithinFatherCell,
518 int numberOfDoFsPerAxisInPatch,
519 int unknowns,
520 const double* __restrict__ coarseGridCellValues,
521 double* __restrict__ fineGridCellValues
522) {
523 static internal::CellInterpolationMap cellInterpolationMap;
524 internal::CellInterpolationOperatorKey key(numberOfDoFsPerAxisInPatch);
525
526 tarch::multicore::Lock lock(interpolationMapSemaphore);
527 if (cellInterpolationMap.count(key) == 0) {
528 cellInterpolationMap[key] = internal::
529 createPiecewiseConstantInterpolationMatrix(numberOfDoFsPerAxisInPatch);
530 }
531 tarch::la::DynamicMatrix* P = cellInterpolationMap[key];
532 lock.free();
533
534#if Dimensions == 2
535 const int patchSize = numberOfDoFsPerAxisInPatch * numberOfDoFsPerAxisInPatch;
536#else
537 const int patchSize = numberOfDoFsPerAxisInPatch * numberOfDoFsPerAxisInPatch
538 * numberOfDoFsPerAxisInPatch;
539#endif
540
542 fineGridCellValues, // image
543 coarseGridCellValues, // preimage
544 unknowns, // batch size, i.e. how often to apply it in one AoS rush
545 patchSize, // result size, i.e. size of image
546 serialiseMarkerIn3x3PatchAssembly(
547 relativePositionWithinFatherCell,
548 numberOfDoFsPerAxisInPatch
549 )
550 );
551}
552
553//
554// =============================
555// Linear interpolation routines
556// =============================
557//
561 int numberOfDoFsPerAxisInPatch,
562 int overlap,
563 int unknowns,
564 const double* __restrict__ coarseGridFaceValues,
565 double* __restrict__ fineGridFaceValues
566 ) {
568 "interpolateHaloLayer_AoS_linear_with_constant_extrapolation(...)",
569 marker.toString(),
570 numberOfDoFsPerAxisInPatch,
571 overlap,
572 unknowns
573 );
574
575 const int normal = marker.getSelectedFaceNumber() % Dimensions;
576
577 static internal::FaceInterpolationMap faceInterpolationMap;
579 key(numberOfDoFsPerAxisInPatch, overlap, normal);
580
581 assertionEquals(overlap, 1);
582
583 tarch::multicore::Lock lock(interpolationMapSemaphore);
584 if (faceInterpolationMap.count(key) == 0) {
585 faceInterpolationMap[key] = internal::createLinearInterpolationMatrix(
586 numberOfDoFsPerAxisInPatch,
587 normal,
588 false,
589 false
590 );
591 }
592 tarch::la::DynamicMatrix* P = faceInterpolationMap[key];
593 lock.free();
594
595#if Dimensions == 2
596 const int patchSize = numberOfDoFsPerAxisInPatch * 2 * overlap;
597#else
598 const int patchSize = numberOfDoFsPerAxisInPatch * numberOfDoFsPerAxisInPatch
599 * 2 * overlap;
600#endif
601
602 logDebug(
603 "interpolateHaloLayer_AoS_linear_with_constant_extrapolation(...)",
604 marker.toString() << ": " << P->toString()
605 );
606
607 P->batchedMultiplyAoS(
608 fineGridFaceValues, // image
609 coarseGridFaceValues, // preimage
610 unknowns, // batch size, i.e. how often to apply it in one AoS rush
611 patchSize, // result size, i.e. size of image
613 marker.getRelativePositionWithinFatherCell(),
614 normal
615 ) * patchSize
616 );
617
618 logTraceOut("interpolateHaloLayer_AoS_linear_with_constant_extrapolation(...)"
619 );
620}
621
625 int numberOfDoFsPerAxisInPatch,
626 int overlap,
627 int unknowns,
628 const double* __restrict__ coarseGridCellValues,
629 const double* __restrict__ coarseGridFaceValues,
630 double* __restrict__ fineGridFaceValues
631 ) {
633 "interpolateHaloLayer_AoS_linear_with_constant_extrapolation(...)",
634 marker.toString(),
635 numberOfDoFsPerAxisInPatch,
636 overlap,
637 unknowns
638 );
639
640 if (marker.isInteriorFaceWithinPatch()) {
641#if Dimensions == 2
642 const int patchSize = numberOfDoFsPerAxisInPatch
643 * numberOfDoFsPerAxisInPatch;
644#else
645 const int patchSize = numberOfDoFsPerAxisInPatch * numberOfDoFsPerAxisInPatch
646 * numberOfDoFsPerAxisInPatch;
647#endif
648
649 double* leftAdjacentPatchData = tarch::allocateMemory<double>(
650 patchSize * unknowns,
652 );
653 double* rightAdjacentPatchData = tarch::allocateMemory<double>(
654 patchSize * unknowns,
656 );
657
659 leftAdjacentPatchIndex = marker.getRelativePositionWithinFatherCell();
661 rightAdjacentPatchIndex = marker.getRelativePositionWithinFatherCell();
662 leftAdjacentPatchIndex(marker.getSelectedFaceNumber() % Dimensions)--;
663
664 internal::interpolateCell_AoS_linear(
665 leftAdjacentPatchIndex,
666 numberOfDoFsPerAxisInPatch,
667 unknowns,
668 coarseGridCellValues,
669 leftAdjacentPatchData,
670 false
671 );
672 internal::interpolateCell_AoS_linear(
673 rightAdjacentPatchIndex,
674 numberOfDoFsPerAxisInPatch,
675 unknowns,
676 coarseGridCellValues,
677 rightAdjacentPatchData,
678 false
679 );
680
681 internal::projectInterpolatedFineCellsOnHaloLayer_AoS(
682 marker,
683 numberOfDoFsPerAxisInPatch,
684 overlap,
685 unknowns,
686 leftAdjacentPatchData,
687 rightAdjacentPatchData,
688 fineGridFaceValues
689 );
690
691 tarch::freeMemory(leftAdjacentPatchData, tarch::MemoryLocation::Heap);
692 tarch::freeMemory(rightAdjacentPatchData, tarch::MemoryLocation::Heap);
693 } else {
695 marker,
696 numberOfDoFsPerAxisInPatch,
697 overlap,
698 unknowns,
699 coarseGridFaceValues,
700 fineGridFaceValues
701 );
702 }
703
704 logTraceOut("interpolateHaloLayer_AoS_linear_with_constant_extrapolation(...)"
705 );
706}
707
711 int numberOfDoFsPerAxisInPatch,
712 int unknowns,
713 const double* __restrict__ coarseGridCellValues,
714 double* __restrict__ fineGridCellValues
715 ) {
717 "interpolateCell_AoS_linear_with_constant_extrapolation(...)",
718 marker.toString(),
719 numberOfDoFsPerAxisInPatch,
720 unknowns
721 );
722
723 internal::interpolateCell_AoS_linear(
724 marker.getRelativePositionWithinFatherCell(),
725 numberOfDoFsPerAxisInPatch,
726 unknowns,
727 coarseGridCellValues,
728 fineGridCellValues,
729 false
730 );
731
732 logTraceOut("interpolateCell_AoS_linear_with_constant_extrapolation(...)");
733}
734
738 int numberOfDoFsPerAxisInPatch,
739 int overlap,
740 int unknowns,
741 const double* __restrict__ coarseGridFaceValues,
742 double* __restrict__ fineGridFaceValues
743 ) {
745 "interpolateHaloLayer_AoS_linear_with_linear_extrapolation(...)",
746 marker.toString(),
747 numberOfDoFsPerAxisInPatch,
748 overlap,
749 unknowns
750 );
751
752 const int normal = marker.getSelectedFaceNumber() % Dimensions;
753
754 static internal::FaceInterpolationMap faceInterpolationMap;
756 key(numberOfDoFsPerAxisInPatch, overlap, normal);
757
758 assertionEquals(overlap, 1);
759
760 tarch::multicore::Lock lock(interpolationMapSemaphore);
761 if (faceInterpolationMap.count(key) == 0) {
762 faceInterpolationMap[key] = internal::createLinearInterpolationMatrix(
763 numberOfDoFsPerAxisInPatch,
764 normal,
765 true,
766 false
767 );
768 }
769 tarch::la::DynamicMatrix* P = faceInterpolationMap[key];
770 lock.free();
771
772#if Dimensions == 2
773 const int patchSize = numberOfDoFsPerAxisInPatch * 2 * overlap;
774#else
775 const int patchSize = numberOfDoFsPerAxisInPatch * numberOfDoFsPerAxisInPatch
776 * 2 * overlap;
777#endif
778
779 logDebug(
780 "interpolateHaloLayer_AoS_linear_with_constant_extrapolation(...)",
781 marker.toString() << ": " << P->toString()
782 );
783
784 P->batchedMultiplyAoS(
785 fineGridFaceValues, // image
786 coarseGridFaceValues, // preimage
787 unknowns, // batch size, i.e. how often to apply it in one AoS rush
788 patchSize, // result size, i.e. size of image
790 marker.getRelativePositionWithinFatherCell(),
791 normal
792 ) * patchSize
793 );
794
795 logTraceOut("interpolateHaloLayer_AoS_linear_with_linear_extrapolation(...)");
796}
797
801 int numberOfDoFsPerAxisInPatch,
802 int overlap,
803 int unknowns,
804 const double* __restrict__ coarseGridCellValues,
805 const double* __restrict__ coarseGridFaceValues,
806 double* __restrict__ fineGridFaceValues
807 ) {
809 "interpolateHaloLayer_AoS_linear_with_linear_extrapolation(...)",
810 marker.toString(),
811 numberOfDoFsPerAxisInPatch,
812 overlap,
813 unknowns
814 );
815
816 if (marker.isInteriorFaceWithinPatch()) {
817#if Dimensions == 2
818 const int patchSize = numberOfDoFsPerAxisInPatch
819 * numberOfDoFsPerAxisInPatch;
820#else
821 const int patchSize = numberOfDoFsPerAxisInPatch * numberOfDoFsPerAxisInPatch
822 * numberOfDoFsPerAxisInPatch;
823#endif
824
825 double* leftAdjacentPatchData = tarch::allocateMemory<double>(
826 patchSize * unknowns,
828 );
829 double* rightAdjacentPatchData = tarch::allocateMemory<double>(
830 patchSize * unknowns,
832 );
833
835 leftAdjacentPatchIndex = marker.getRelativePositionWithinFatherCell();
837 rightAdjacentPatchIndex = marker.getRelativePositionWithinFatherCell();
838 leftAdjacentPatchIndex(marker.getSelectedFaceNumber() % Dimensions)--;
839
840 internal::interpolateCell_AoS_linear(
841 leftAdjacentPatchIndex,
842 numberOfDoFsPerAxisInPatch,
843 unknowns,
844 coarseGridCellValues,
845 leftAdjacentPatchData,
846 true
847 );
848 internal::interpolateCell_AoS_linear(
849 rightAdjacentPatchIndex,
850 numberOfDoFsPerAxisInPatch,
851 unknowns,
852 coarseGridCellValues,
853 rightAdjacentPatchData,
854 true
855 );
856
857 internal::projectInterpolatedFineCellsOnHaloLayer_AoS(
858 marker,
859 numberOfDoFsPerAxisInPatch,
860 overlap,
861 unknowns,
862 leftAdjacentPatchData,
863 rightAdjacentPatchData,
864 fineGridFaceValues
865 );
866
867 tarch::freeMemory(leftAdjacentPatchData, tarch::MemoryLocation::Heap);
868 tarch::freeMemory(rightAdjacentPatchData, tarch::MemoryLocation::Heap);
869 } else {
871 marker,
872 numberOfDoFsPerAxisInPatch,
873 overlap,
874 unknowns,
875 coarseGridFaceValues,
876 fineGridFaceValues
877 );
878 }
879
880 logTraceOut("interpolateHaloLayer_AoS_linear_with_linear_extrapolation(...)");
881}
882
886 int numberOfDoFsPerAxisInPatch,
887 int unknowns,
888 const double* __restrict__ coarseGridCellValues,
889 double* __restrict__ fineGridCellValues
890 ) {
892 "interpolateCell_AoS_linear_with_linear_extrapolation(...)",
893 marker.toString(),
894 numberOfDoFsPerAxisInPatch,
895 unknowns
896 );
897
898 internal::interpolateCell_AoS_linear(
899 marker.getRelativePositionWithinFatherCell(),
900 numberOfDoFsPerAxisInPatch,
901 unknowns,
902 coarseGridCellValues,
903 fineGridCellValues,
904 true
905 );
906
907 logTraceOut("interpolateCell_AoS_linear_with_linear_extrapolation(...)");
908}
909
913 int numberOfDoFsPerAxisInPatch,
914 int overlap,
915 int unknowns,
916 const double* __restrict__ coarseGridFaceValues,
917 double* __restrict__ fineGridFaceValues
918 ) {
920 "interpolateHaloLayer_AoS_linear_with_constant_extrapolation_and_linear_normal_interpolation(...)",
921 marker.toString(),
922 numberOfDoFsPerAxisInPatch,
923 overlap,
924 unknowns
925 );
926
927 const int normal = marker.getSelectedFaceNumber() % Dimensions;
928
929 static internal::FaceInterpolationMap faceInterpolationMap;
931 key(numberOfDoFsPerAxisInPatch, overlap, normal);
932
933 assertionEquals(overlap, 1);
934
935 tarch::multicore::Lock lock(interpolationMapSemaphore);
936 if (faceInterpolationMap.count(key) == 0) {
937 faceInterpolationMap[key] = internal::createLinearInterpolationMatrix(
938 numberOfDoFsPerAxisInPatch,
939 normal,
940 false,
941 true
942 );
943 }
944 tarch::la::DynamicMatrix* P = faceInterpolationMap[key];
945 lock.free();
946
947#if Dimensions == 2
948 const int patchSize = numberOfDoFsPerAxisInPatch * 2 * overlap;
949#else
950 const int patchSize = numberOfDoFsPerAxisInPatch * numberOfDoFsPerAxisInPatch
951 * 2 * overlap;
952#endif
953
954 logDebug(
955 "interpolateHaloLayer_AoS_linear_with_constant_extrapolation_and_linear_normal_interpolation(...)",
956 marker.toString() << ": " << P->toString()
957 );
958
959 P->batchedMultiplyAoS(
960 fineGridFaceValues, // image
961 coarseGridFaceValues, // preimage
962 unknowns, // batch size, i.e. how often to apply it in one AoS rush
963 patchSize, // result size, i.e. size of image
965 marker.getRelativePositionWithinFatherCell(),
966 normal
967 ) * patchSize
968 );
969
971 "interpolateHaloLayer_AoS_linear_with_constant_extrapolation_and_linear_normal_interpolation(...)"
972 );
973}
974
978 int numberOfDoFsPerAxisInPatch,
979 int overlap,
980 int unknowns,
981 const double* __restrict__ coarseGridCellValues,
982 const double* __restrict__ coarseGridFaceValues,
983 double* __restrict__ fineGridFaceValues
984 ) {
986 "interpolateHaloLayer_AoS_linear_with_constant_extrapolation_and_linear_normal_interpolation(...)",
987 marker.toString(),
988 numberOfDoFsPerAxisInPatch,
989 overlap,
990 unknowns
991 );
992
993 if (marker.isInteriorFaceWithinPatch()) {
994#if Dimensions == 2
995 const int patchSize = numberOfDoFsPerAxisInPatch
996 * numberOfDoFsPerAxisInPatch;
997#else
998 const int patchSize = numberOfDoFsPerAxisInPatch * numberOfDoFsPerAxisInPatch
999 * numberOfDoFsPerAxisInPatch;
1000#endif
1001
1002 double* leftAdjacentPatchData = tarch::allocateMemory<double>(
1003 patchSize * unknowns,
1005 );
1006 double* rightAdjacentPatchData = tarch::allocateMemory<double>(
1007 patchSize * unknowns,
1009 );
1010
1012 leftAdjacentPatchIndex = marker.getRelativePositionWithinFatherCell();
1014 rightAdjacentPatchIndex = marker.getRelativePositionWithinFatherCell();
1015 leftAdjacentPatchIndex(marker.getSelectedFaceNumber() % Dimensions)--;
1016
1017 internal::interpolateCell_AoS_linear(
1018 leftAdjacentPatchIndex,
1019 numberOfDoFsPerAxisInPatch,
1020 unknowns,
1021 coarseGridCellValues,
1022 leftAdjacentPatchData,
1023 false
1024 );
1025 internal::interpolateCell_AoS_linear(
1026 rightAdjacentPatchIndex,
1027 numberOfDoFsPerAxisInPatch,
1028 unknowns,
1029 coarseGridCellValues,
1030 rightAdjacentPatchData,
1031 false
1032 );
1033
1034 internal::projectInterpolatedFineCellsOnHaloLayer_AoS(
1035 marker,
1036 numberOfDoFsPerAxisInPatch,
1037 overlap,
1038 unknowns,
1039 leftAdjacentPatchData,
1040 rightAdjacentPatchData,
1041 fineGridFaceValues
1042 );
1043
1044 tarch::freeMemory(leftAdjacentPatchData, tarch::MemoryLocation::Heap);
1045 tarch::freeMemory(rightAdjacentPatchData, tarch::MemoryLocation::Heap);
1046 } else {
1048 marker,
1049 numberOfDoFsPerAxisInPatch,
1050 overlap,
1051 unknowns,
1052 coarseGridFaceValues,
1053 fineGridFaceValues
1054 );
1055 }
1056
1058 "interpolateHaloLayer_AoS_linear_with_constant_extrapolation_and_linear_normal_interpolation(...)"
1059 );
1060}
1061
1065 int numberOfDoFsPerAxisInPatch,
1066 int unknowns,
1067 const double* __restrict__ coarseGridCellValues,
1068 double* __restrict__ fineGridCellValues
1069 ) {
1071 "interpolateCell_AoS_linear_with_constant_extrapolation_and_linear_normal_interpolation(...)",
1072 marker.toString(),
1073 numberOfDoFsPerAxisInPatch,
1074 unknowns
1075 );
1076
1077 internal::interpolateCell_AoS_linear(
1078 marker.getRelativePositionWithinFatherCell(),
1079 numberOfDoFsPerAxisInPatch,
1080 unknowns,
1081 coarseGridCellValues,
1082 fineGridCellValues,
1083 false
1084 );
1085
1087 "interpolateCell_AoS_linear_with_constant_extrapolation_and_linear_normal_interpolation(...)"
1088 );
1089}
1090
1094 int numberOfDoFsPerAxisInPatch,
1095 int overlap,
1096 int unknowns,
1097 const double* __restrict__ coarseGridFaceValues,
1098 double* __restrict__ fineGridFaceValues
1099 ) {
1101 "interpolateHaloLayer_AoS_linear_with_linear_extrapolation_and_linear_normal_interpolation(...)",
1102 marker.toString(),
1103 numberOfDoFsPerAxisInPatch,
1104 overlap,
1105 unknowns
1106 );
1107
1108 const int normal = marker.getSelectedFaceNumber() % Dimensions;
1109
1110 static internal::FaceInterpolationMap faceInterpolationMap;
1112 key(numberOfDoFsPerAxisInPatch, overlap, normal);
1113
1114 assertionEquals(overlap, 1);
1115
1116 tarch::multicore::Lock lock(interpolationMapSemaphore);
1117 if (faceInterpolationMap.count(key) == 0) {
1118 faceInterpolationMap[key] = internal::createLinearInterpolationMatrix(
1119 numberOfDoFsPerAxisInPatch,
1120 normal,
1121 true,
1122 true
1123 );
1124 }
1125 tarch::la::DynamicMatrix* P = faceInterpolationMap[key];
1126 lock.free();
1127
1128#if Dimensions == 2
1129 const int patchSize = numberOfDoFsPerAxisInPatch * 2 * overlap;
1130#else
1131 const int patchSize = numberOfDoFsPerAxisInPatch * numberOfDoFsPerAxisInPatch
1132 * 2 * overlap;
1133#endif
1134
1135 logDebug(
1136 "interpolateHaloLayer_AoS_linear_with_constant_extrapolation_and_linear_normal_interpolation(...)",
1137 marker.toString() << ": " << P->toString()
1138 );
1139
1140 P->batchedMultiplyAoS(
1141 fineGridFaceValues, // image
1142 coarseGridFaceValues, // preimage
1143 unknowns, // batch size, i.e. how often to apply it in one AoS rush
1144 patchSize, // result size, i.e. size of image
1146 marker.getRelativePositionWithinFatherCell(),
1147 normal
1148 ) * patchSize
1149 );
1150
1152 "interpolateHaloLayer_AoS_linear_with_linear_extrapolation_and_linear_normal_interpolation(...)"
1153 );
1154}
1155
1159 int numberOfDoFsPerAxisInPatch,
1160 int overlap,
1161 int unknowns,
1162 const double* __restrict__ coarseGridCellValues,
1163 const double* __restrict__ coarseGridFaceValues,
1164 double* __restrict__ fineGridFaceValues
1165 ) {
1167 "interpolateHaloLayer_AoS_linear_with_linear_extrapolation_and_linear_normal_interpolation(...)",
1168 marker.toString(),
1169 numberOfDoFsPerAxisInPatch,
1170 overlap,
1171 unknowns
1172 );
1173
1174 if (marker.isInteriorFaceWithinPatch()) {
1175#if Dimensions == 2
1176 const int patchSize = numberOfDoFsPerAxisInPatch
1177 * numberOfDoFsPerAxisInPatch;
1178#else
1179 const int patchSize = numberOfDoFsPerAxisInPatch * numberOfDoFsPerAxisInPatch
1180 * numberOfDoFsPerAxisInPatch;
1181#endif
1182
1183 double* leftAdjacentPatchData = tarch::allocateMemory<double>(
1184 patchSize * unknowns,
1186 );
1187 double* rightAdjacentPatchData = tarch::allocateMemory<double>(
1188 patchSize * unknowns,
1190 );
1191
1193 leftAdjacentPatchIndex = marker.getRelativePositionWithinFatherCell();
1195 rightAdjacentPatchIndex = marker.getRelativePositionWithinFatherCell();
1196 leftAdjacentPatchIndex(marker.getSelectedFaceNumber() % Dimensions)--;
1197
1198 internal::interpolateCell_AoS_linear(
1199 leftAdjacentPatchIndex,
1200 numberOfDoFsPerAxisInPatch,
1201 unknowns,
1202 coarseGridCellValues,
1203 leftAdjacentPatchData,
1204 true
1205 );
1206 internal::interpolateCell_AoS_linear(
1207 rightAdjacentPatchIndex,
1208 numberOfDoFsPerAxisInPatch,
1209 unknowns,
1210 coarseGridCellValues,
1211 rightAdjacentPatchData,
1212 true
1213 );
1214
1215 internal::projectInterpolatedFineCellsOnHaloLayer_AoS(
1216 marker,
1217 numberOfDoFsPerAxisInPatch,
1218 overlap,
1219 unknowns,
1220 leftAdjacentPatchData,
1221 rightAdjacentPatchData,
1222 fineGridFaceValues
1223 );
1224
1225 tarch::freeMemory(leftAdjacentPatchData, tarch::MemoryLocation::Heap);
1226 tarch::freeMemory(rightAdjacentPatchData, tarch::MemoryLocation::Heap);
1227 } else {
1229 marker,
1230 numberOfDoFsPerAxisInPatch,
1231 overlap,
1232 unknowns,
1233 coarseGridFaceValues,
1234 fineGridFaceValues
1235 );
1236 }
1237
1239 "interpolateHaloLayer_AoS_linear_with_linear_extrapolation_and_linear_normal_interpolation(...)"
1240 );
1241}
1242
1243
1246 int numberOfDoFsPerAxisInSourcePatch,
1247 int numberOfDoFsPerAxisInDestinationPatch,
1248 int haloSourcePatch,
1249 int haloDestinationPatch,
1250 int unknowns,
1251 const double* __restrict__ sourceValues,
1252 double* __restrict__ destinationValues,
1253 ::peano4::utils::LoopPlacement parallelisation
1254 ) {
1255 assertion(
1256 numberOfDoFsPerAxisInDestinationPatch > numberOfDoFsPerAxisInSourcePatch
1257 );
1258
1259 simtDforWithSchedulerInstructions(
1260 destinationVolume,
1261 numberOfDoFsPerAxisInDestinationPatch + 2 * haloDestinationPatch,
1262 parallelisation
1263 ) const int baseIndexDestination
1265 destinationVolume,
1266 numberOfDoFsPerAxisInDestinationPatch + 2 * haloDestinationPatch
1267 )
1268 * unknowns;
1269 for (int unknown = 0; unknown < unknowns; unknown++) {
1270 destinationValues[baseIndexDestination + unknown] = 0.0;
1271 }
1272
1273 dfor(sourceVolume, numberOfDoFsPerAxisInSourcePatch + 2 * haloSourcePatch) {
1275 normalisedSourcePosition = tarch::la::multiplyComponents(
1276 tarch::la::
1277 Vector<Dimensions, double>(1.0 / numberOfDoFsPerAxisInSourcePatch),
1278 tarch::la::convertScalar<double>(sourceVolume
1279 ) + tarch::la::Vector<Dimensions, double>(0.5 - haloSourcePatch)
1280 );
1282 normalisedDestinationPosition = tarch::la::multiplyComponents(
1284 1.0 / numberOfDoFsPerAxisInDestinationPatch
1285 ),
1286 tarch::la::convertScalar<double>(destinationVolume
1287 ) + tarch::la::Vector<Dimensions, double>(0.5 - haloDestinationPatch)
1288 );
1289
1290 int outsidePatchAlongCoordinateAxis = 0;
1291 for (int d = 0; d < Dimensions; d++) {
1292 // outsidePatchAlongCoordinateAxis += ( sourceVolume(d) <
1293 // haloSourcePatch ) ? 1 : 0; outsidePatchAlongCoordinateAxis += (
1294 // sourceVolume(d) >=
1295 // haloSourcePatch+numberOfDoFsPerAxisInSourcePatch ) ? 1 : 0;
1296 if (sourceVolume(d) < haloSourcePatch)
1297 outsidePatchAlongCoordinateAxis++;
1298 if (sourceVolume(d) >= haloSourcePatch + numberOfDoFsPerAxisInSourcePatch)
1299 outsidePatchAlongCoordinateAxis++;
1300 }
1301 const bool isFromUnitialisedDiagonal = outsidePatchAlongCoordinateAxis > 1;
1302
1303 double weight = 1.0;
1304 for (int d = 0; d < Dimensions; d++) {
1305 double distanceAlongCoordinateAxis = std::abs(
1306 normalisedDestinationPosition(d) - normalisedSourcePosition(d)
1307 );
1308 double h = 1.0 / static_cast<double>(numberOfDoFsPerAxisInSourcePatch);
1309 distanceAlongCoordinateAxis = 1.0 - distanceAlongCoordinateAxis / h;
1310 weight *= std::max(0.0, distanceAlongCoordinateAxis);
1311 }
1312#if !defined(SharedSYCL)
1313 assertion(tarch::la::greaterEquals(weight, 0.0));
1314 assertion(tarch::la::smallerEquals(weight, 1.0));
1315#endif
1316
1317 tarch::la::Vector<Dimensions, int> amendedSourceVolume = sourceVolume;
1318 if (isFromUnitialisedDiagonal) {
1319#pragma unroll
1320 for (int d = 0; d < Dimensions; d++) {
1321 amendedSourceVolume(d
1322 ) = std::max(amendedSourceVolume(d), haloSourcePatch);
1323 amendedSourceVolume(d) = std::min(
1324 amendedSourceVolume(d),
1325 haloSourcePatch + numberOfDoFsPerAxisInSourcePatch - 1
1326 );
1327 }
1328 }
1329
1330 int baseIndexSource = peano4::utils::dLinearised(
1331 amendedSourceVolume,
1332 numberOfDoFsPerAxisInSourcePatch + 2 * haloSourcePatch
1333 )
1334 * unknowns;
1335 for (int unknown = 0; unknown < unknowns; unknown++) {
1336 destinationValues[baseIndexDestination + unknown]
1337 += weight * sourceValues[baseIndexSource + unknown];
1338 }
1339 }
1341}
1342
1343
1346 int numberOfDoFsPerAxisInSourcePatch,
1347 int numberOfDoFsPerAxisInDestinationPatch,
1348 int haloSourcePatch,
1349 int haloDestinationPatch,
1350 int unknowns,
1351 const double* __restrict__ interpolationData,
1352 const int* __restrict__ columnIndices,
1353 const int* __restrict__ rowIndices,
1354 const double* __restrict__ sourceValues,
1355 double* __restrict__ destinationValues,
1356 ::peano4::utils::LoopPlacement parallelisation
1357 ) {
1358 int destinationSize
1359 = (numberOfDoFsPerAxisInDestinationPatch + 2 * haloDestinationPatch)
1360 * (numberOfDoFsPerAxisInDestinationPatch + 2 * haloDestinationPatch)
1361 * (numberOfDoFsPerAxisInDestinationPatch + 2 * haloDestinationPatch);
1362
1363 for (int row = 0; row < destinationSize; row++) {
1364 for (int unknown = 0; unknown < unknowns; unknown++) {
1365 destinationValues[row * unknowns + unknown] = 0.0;
1366 }
1367 for (int column = rowIndices[row]; column < rowIndices[row + 1]; column++) {
1368 for (int unknown = 0; unknown < unknowns; unknown++) {
1369 destinationValues[row * unknowns + unknown]
1370 += interpolationData[column]
1371 * sourceValues[columnIndices[column] * unknowns + unknown];
1372 }
1373 }
1374 }
1375}
1376
1379 int numberOfDoFsPerAxisInSourcePatch,
1380 int numberOfDoFsPerAxisInDestinationPatch,
1381 int haloSourcePatch,
1382 int haloDestinationPatch,
1383 int unknowns,
1384 const double* __restrict__ sourceValues,
1385 double* __restrict__ destinationValues,
1386 ::peano4::utils::LoopPlacement parallelisation
1387 ) {
1388
1389 int indices[18];
1391 double* centerValues = new double[unknowns];
1392 int numberOfDerivatives = (Dimensions == 2) ? 5 : 9;
1393
1394 dfor(volume, numberOfDoFsPerAxisInSourcePatch + 2) {
1396 sourceVolume = volume
1397 + tarch::la::Vector<Dimensions, int>(haloSourcePatch - 1);
1398 int outsidePatchAlongCoordinateAxis = 0;
1399 for (int d = 0; d < Dimensions; d++) {
1400 if (sourceVolume(d) < haloSourcePatch || sourceVolume(d) >= numberOfDoFsPerAxisInSourcePatch + haloSourcePatch) {
1401 outsidePatchAlongCoordinateAxis++;
1402 }
1403 }
1404
1405 tarch::la::Vector<Dimensions, int> stencilOffset(0);
1406 for (int d = 0; d < Dimensions; d++) {
1407 if (sourceVolume(d) == 0)
1408 stencilOffset(d) = 1;
1409 if (sourceVolume(d) == numberOfDoFsPerAxisInSourcePatch + 2 * haloSourcePatch - 1)
1410 stencilOffset = -1;
1411 }
1412 for (int d1 = 0; d1 < Dimensions; d1++) {
1413 for (int d2 = d1 + 1; d2 < Dimensions; d2++) {
1414 if (sourceVolume(d1) == haloSourcePatch and sourceVolume(d2) == haloSourcePatch) {
1415 stencilOffset(d1) = 1;
1416 stencilOffset(d2) = 1;
1417 }
1418 if (sourceVolume(d1) == haloSourcePatch and sourceVolume(d2) == haloSourcePatch + numberOfDoFsPerAxisInSourcePatch - 1) {
1419 stencilOffset(d1) = 1;
1420 stencilOffset(d2) = -1;
1421 }
1422 if (sourceVolume(d1) == haloSourcePatch + numberOfDoFsPerAxisInSourcePatch - 1 and sourceVolume(d2) == haloSourcePatch) {
1423 stencilOffset(d1) = -1;
1424 stencilOffset(d2) = 1;
1425 }
1426 if (sourceVolume(d1) == haloSourcePatch + numberOfDoFsPerAxisInSourcePatch - 1 and sourceVolume(d2) == haloSourcePatch + numberOfDoFsPerAxisInSourcePatch - 1) {
1427 stencilOffset(d1) = -1;
1428 stencilOffset(d2) = -1;
1429 }
1430 }
1431 }
1432
1433 if (outsidePatchAlongCoordinateAxis <= 1) {
1434 int stencilSize = 0;
1435 dfor(offset, 3) {
1436 if (
1437 offset != tarch::la::Vector<Dimensions, int>(1) - stencilOffset and
1438#if Dimensions == 3
1439 !(offset(0) + stencilOffset(0) != 1
1440 and offset(1) + stencilOffset(1) != 1
1441 and offset(2) + stencilOffset(2) != 1)
1442#elif Dimensions == 2
1443 !(offset(0) + stencilOffset(0) != 1
1444 and offset(1) + stencilOffset(1) != 1)
1445#endif
1446 ) {
1447 stencil[stencilSize] = stencilOffset + offset
1449 indices[stencilSize] = peano4::utils::dLinearised(
1450 sourceVolume + offset + stencilOffset
1452 numberOfDoFsPerAxisInSourcePatch + 2 * haloSourcePatch
1453 );
1454 stencilSize++;
1455 }
1456 }
1457
1458 int centerIndex = peano4::utils::dLinearised(
1459 sourceVolume,
1460 numberOfDoFsPerAxisInSourcePatch + 2 * haloSourcePatch
1461 );
1462 for (int unknown = 0; unknown < unknowns; unknown++) {
1463 centerValues[unknown] = sourceValues[centerIndex * unknowns + unknown];
1464 }
1465
1466 tarch::la::DynamicMatrix deltas(stencilSize, unknowns);
1467 for (int k = 0; k < stencilSize; k++) {
1468 for (int unknown = 0; unknown < unknowns; unknown++) {
1469 deltas(
1470 k,
1471 unknown
1472 ) = sourceValues[indices[k] * unknowns + unknown]
1473 - centerValues[unknown];
1474 }
1475 }
1476 tarch::la::DynamicMatrix A(stencilSize, numberOfDerivatives);
1477 for (int row = 0; row < stencilSize; row++) {
1478 int column = 0;
1479 tarch::la::Vector<Dimensions, int>& cell = stencil[row];
1480 for (int a = 0; a < Dimensions; a++) {
1481 A(row, column) = cell(a);
1482 column++;
1483 }
1484 for (int a = 0; a < Dimensions; a++) {
1485 for (int b = a; b < Dimensions; b++) {
1486 A(row, column) = 0.5 * cell(a) * cell(b);
1487 column++;
1488 }
1489 }
1490 }
1491
1495
1496 auto Q_T = tarch::la::transpose(Q);
1497 tarch::la::DynamicMatrix c = Q_T * deltas;
1498
1500 double* result = tarch::allocateMemory(
1501 inverseR.rows() * c.cols(),
1503 );
1504 inverseR.multiplyBySmallMatrix(result, c);
1505 tarch::la::DynamicMatrix x(inverseR.rows(), c.cols(), result);
1506
1507 tarch::la::Vector<Dimensions, int> destinationCorner1
1508 = (sourceVolume - tarch::la::Vector<Dimensions, int>(haloSourcePatch)
1509 ) * numberOfDoFsPerAxisInDestinationPatch
1510 / numberOfDoFsPerAxisInSourcePatch
1511 + tarch::la::Vector<Dimensions, int>(haloDestinationPatch);
1512 tarch::la::Vector<Dimensions, int> destinationCorner2
1513 = (sourceVolume
1514 - tarch::la::Vector<Dimensions, int>(haloSourcePatch - 1)
1515 ) * numberOfDoFsPerAxisInDestinationPatch
1516 / numberOfDoFsPerAxisInSourcePatch
1517 + tarch::la::Vector<Dimensions, int>(haloDestinationPatch);
1518 for (int d = 0; d < Dimensions; d++) {
1519 destinationCorner1(d) = std::max(0, destinationCorner1(d));
1520 destinationCorner2(d) = std::min(
1521 numberOfDoFsPerAxisInDestinationPatch + 2 * haloDestinationPatch,
1522 destinationCorner2(d)
1523 );
1524 }
1525
1527 destinationSpan = destinationCorner2 - destinationCorner1;
1528 int destinationDimensionsProduct = 1;
1529 for (int d = 0; d < Dimensions; d++)
1530 destinationDimensionsProduct *= destinationSpan(d);
1531
1532 int row = 0;
1533 tarch::la::Vector<Dimensions, double> normalisedCornerPosition
1534 = (1.0 / numberOfDoFsPerAxisInDestinationPatch)
1535 * (tarch::la::convertScalar<double>(destinationCorner1) + tarch::la::Vector<Dimensions, double>(0.5 - haloDestinationPatch));
1536 tarch::la::Vector<Dimensions, double> normalisedSourcePosition
1537 = (1.0 / numberOfDoFsPerAxisInSourcePatch)
1538 * (tarch::la::convertScalar<double>(sourceVolume) + tarch::la::Vector<Dimensions, double>(0.5 - haloSourcePatch));
1540 cornerDisplacement = (normalisedCornerPosition
1541 - normalisedSourcePosition)
1542 * (double)numberOfDoFsPerAxisInSourcePatch;
1543
1544 double* matrixData = new double
1545 [destinationDimensionsProduct * numberOfDerivatives];
1546 int writeIndex = 0;
1547
1548 dfor(localDestinationVolume, destinationSpan) {
1550 displacement = cornerDisplacement
1551 + ((double)numberOfDoFsPerAxisInSourcePatch
1552 / numberOfDoFsPerAxisInDestinationPatch)
1554 double>(localDestinationVolume);
1555
1556 for (int a = 0; a < Dimensions; a++) {
1557 matrixData[writeIndex] = displacement(a);
1558 writeIndex++;
1559 }
1560 for (int a = 0; a < Dimensions; a++) {
1561 for (int b = a; b < Dimensions; b++) {
1562 matrixData[writeIndex] = 0.5 * displacement(a) * displacement(b);
1563 writeIndex++;
1564 }
1565 }
1566 }
1567
1568 tarch::la::DynamicMatrix destinationMatrix(
1569 destinationDimensionsProduct,
1570 numberOfDerivatives,
1571 matrixData
1572 );
1573 tarch::la::DynamicMatrix destinationValuesDelta = destinationMatrix * x;
1574
1575 int destinationVolumeNumber = 0;
1576 dfor(localDestinationVolume, destinationSpan) {
1578 destinationVolume = destinationCorner1 + localDestinationVolume;
1579 int destinationIndex = peano4::utils::dLinearised(
1580 destinationVolume,
1581 numberOfDoFsPerAxisInDestinationPatch + 2 * haloDestinationPatch
1582 );
1583 for (int unknown = 0; unknown < unknowns; unknown++) {
1584 destinationValues[destinationIndex * unknowns + unknown]
1585 = centerValues[unknown]
1586 + destinationValuesDelta(destinationVolumeNumber, unknown);
1587 }
1588 destinationVolumeNumber++;
1589 }
1590 }
1591 }
1592 delete[] centerValues;
1593}
1594
1597 int numberOfDoFsPerAxisInSourcePatch,
1598 int numberOfDoFsPerAxisInDestinationPatch,
1599 int haloSourcePatch,
1600 int haloDestinationPatch,
1601 int unknowns,
1602 const double* __restrict__ sourceValues,
1603 double* __restrict__ destinationValues,
1604 ::peano4::utils::LoopPlacement parallelisation
1605 ) {
1606 assertion(
1607 numberOfDoFsPerAxisInDestinationPatch > numberOfDoFsPerAxisInSourcePatch
1608 );
1609
1610 simtDforWithSchedulerInstructions(
1611 destinationVolume,
1612 numberOfDoFsPerAxisInDestinationPatch + 2 * haloDestinationPatch,
1613 parallelisation
1614 ) const int baseIndexDestination
1616 destinationVolume,
1617 numberOfDoFsPerAxisInDestinationPatch + 2 * haloDestinationPatch
1618 )
1619 * unknowns;
1620 for (int unknown = 0; unknown < unknowns; unknown++) {
1621 destinationValues[baseIndexDestination + unknown] = 0.0;
1622 }
1624}
1625
1626
1630 int numberOfDoFsPerAxisInPatch,
1631 int unknowns,
1632 const double* __restrict__ coarseGridCellValues,
1633 double* __restrict__ fineGridCellValues
1634 ) {
1636 "interpolateCell_AoS_linear_with_linear_extrapolation_and_linear_normal_interpolation(...)",
1637 marker.toString(),
1638 numberOfDoFsPerAxisInPatch,
1639 unknowns
1640 );
1641
1642 internal::interpolateCell_AoS_linear(
1643 marker.getRelativePositionWithinFatherCell(),
1644 numberOfDoFsPerAxisInPatch,
1645 unknowns,
1646 coarseGridCellValues,
1647 fineGridCellValues,
1648 true
1649 );
1650
1652 "interpolateCell_AoS_linear_with_linear_extrapolation_and_linear_normal_interpolation(...)"
1653 );
1654}
1655
1657 const tarch::la::Vector<Dimensions, int>& relativePositionWithinFatherCell,
1658 int numberOfDoFsPerAxisInPatch,
1659 int unknowns,
1660 const double* __restrict__ coarseGridCellValues,
1661 double* __restrict__ fineGridCellValues,
1662 bool extrapolateLinearly
1663) {
1664 static internal::CellInterpolationMap cellInterpolationMap;
1665 internal::CellInterpolationOperatorKey key(numberOfDoFsPerAxisInPatch);
1666
1667 tarch::multicore::Lock lock(interpolationMapSemaphore);
1668 if (cellInterpolationMap.count(key) == 0) {
1669 cellInterpolationMap[key] = internal::createLinearInterpolationMatrix(
1670 numberOfDoFsPerAxisInPatch,
1671 extrapolateLinearly
1672 );
1673 }
1674 tarch::la::DynamicMatrix* P = cellInterpolationMap[key];
1675 lock.free();
1676
1677#if Dimensions == 2
1678 const int patchSize = numberOfDoFsPerAxisInPatch * numberOfDoFsPerAxisInPatch;
1679#else
1680 const int patchSize = numberOfDoFsPerAxisInPatch * numberOfDoFsPerAxisInPatch
1681 * numberOfDoFsPerAxisInPatch;
1682#endif
1683
1685 fineGridCellValues, // image
1686 coarseGridCellValues, // preimage
1687 unknowns, // batch size, i.e. how often to apply it in one AoS rush
1688 patchSize, // result size, i.e. size of image
1689 serialiseMarkerIn3x3PatchAssembly(
1690 relativePositionWithinFatherCell,
1691 numberOfDoFsPerAxisInPatch
1692 )
1693 );
1694}
1695
1698 int numberOfDoFsPerAxisInPatch,
1699 int overlap,
1700 int unknowns,
1701 const double* __restrict__ coarseGridFaceValues,
1702 double* __restrict__ fineGridFaceValues
1703) {
1704 const int normal = marker.getSelectedFaceNumber() % Dimensions;
1705 const int i = (normal + 1) % Dimensions;
1706#if Dimensions == 3
1707 const int j = (normal + 2) % Dimensions;
1708#endif
1709
1711 coarseStart(i
1712 ) = marker.getRelativePositionWithinFatherCell()(i)
1713 * numberOfDoFsPerAxisInPatch / 3;
1714#if Dimensions == 3
1715 coarseStart(j
1716 ) = marker.getRelativePositionWithinFatherCell()(j)
1717 * numberOfDoFsPerAxisInPatch / 3;
1718#endif
1719
1721 coarseEnd(i
1722 ) = ((marker.getRelativePositionWithinFatherCell()(i) + 1
1723 ) * numberOfDoFsPerAxisInPatch
1724 + 2)
1725 / 3;
1726#if Dimensions == 3
1727 coarseEnd(j
1728 ) = ((marker.getRelativePositionWithinFatherCell()(j) + 1
1729 ) * numberOfDoFsPerAxisInPatch
1730 + 2)
1731 / 3;
1732#endif
1733
1734 tarch::la::Vector<Dimensions, int> fineDataStart(0);
1735 fineDataStart(i
1736 ) = (marker.getRelativePositionWithinFatherCell()(i)
1737 * numberOfDoFsPerAxisInPatch)
1738 % 3;
1739#if Dimensions == 3
1740 fineDataStart(j
1741 ) = (marker.getRelativePositionWithinFatherCell()(j)
1742 * numberOfDoFsPerAxisInPatch)
1743 % 3;
1744#endif
1745
1746 tarch::la::Vector<Dimensions, int> faceDimensions = coarseEnd - coarseStart;
1747 faceDimensions(normal) = (overlap - 1) / 3 + 1;
1748
1749 int indices[18];
1751 double* centerValues = new double[unknowns];
1752
1753 dfor(kCoarse, faceDimensions) {
1754 tarch::la::Vector<Dimensions, int> center = kCoarse + coarseStart;
1755 center(normal
1756 ) = (marker.getSelectedFaceNumber() < Dimensions)
1757 ? overlap - kCoarse(normal) - 1
1758 : overlap + kCoarse(normal);
1759 tarch::la::Vector<Dimensions, int> stencilOffset(0);
1760
1761 if (center(i) == 0)
1762 stencilOffset(i) = 1;
1763 if (center(i) == numberOfDoFsPerAxisInPatch - 1)
1764 stencilOffset(i) = -1;
1765#if Dimensions == 3
1766 if (center(j) == 0)
1767 stencilOffset(j) = 1;
1768 if (center(j) == numberOfDoFsPerAxisInPatch - 1)
1769 stencilOffset(j) = -1;
1770#endif
1771
1772 int stencilSize = 0;
1773 dfor(offset, 3) {
1774 if (
1775 offset != tarch::la::Vector<Dimensions, int>(1) - stencilOffset and
1776#if Dimensions == 3
1777 !(offset(0) + stencilOffset(0) != 1
1778 and offset(1) + stencilOffset(1) != 1
1779 and offset(2) + stencilOffset(2) != 1)
1780#elif Dimensions == 2
1781 !(offset(0) + stencilOffset(0) != 1
1782 and offset(1) + stencilOffset(1) != 1)
1783#endif
1784 ) {
1785 stencil[stencilSize] = stencilOffset + offset
1787 indices[stencilSize] = serialiseVoxelIndexInOverlap(
1788 center + offset + stencilOffset
1790 numberOfDoFsPerAxisInPatch,
1791 overlap,
1792 normal
1793 );
1794 stencilSize++;
1795 }
1796 }
1797
1799 center,
1800 numberOfDoFsPerAxisInPatch,
1801 overlap,
1802 normal
1803 );
1804 for (int unknown = 0; unknown < unknowns; unknown++) {
1805 centerValues[unknown] = coarseGridFaceValues[index * unknowns + unknown];
1806 }
1807
1808 tarch::la::DynamicMatrix deltas(stencilSize, unknowns);
1809 for (int k = 0; k < stencilSize; k++) {
1810 for (int unknown = 0; unknown < unknowns; unknown++) {
1811 deltas(
1812 k,
1813 unknown
1814 ) = coarseGridFaceValues[indices[k] * unknowns + unknown]
1815 - centerValues[unknown];
1816 }
1817 }
1818
1819 int numberOfDerivatives = (Dimensions == 2) ? 5 : 9;
1820 tarch::la::DynamicMatrix A(stencilSize, numberOfDerivatives);
1821 for (int row = 0; row < stencilSize; row++) {
1822 int column = 0;
1824 for (int a = 0; a < Dimensions; a++) {
1825 A(row, column) = cell(a);
1826 column++;
1827 }
1828 for (int a = 0; a < Dimensions; a++) {
1829 for (int b = a; b < Dimensions; b++) {
1830 A(row, column) = 0.5 * cell(a) * cell(b);
1831 column++;
1832 }
1833 }
1834 }
1835
1836 tarch::la::DynamicMatrix Q(A.rows(), A.cols());
1837 tarch::la::DynamicMatrix R(A.cols(), A.cols());
1838
1840
1841 auto Q_T = tarch::la::transpose(Q);
1842 tarch::la::DynamicMatrix c = Q_T * deltas;
1843
1845 tarch::la::DynamicMatrix x = inverseR * c;
1846
1849
1850 if (kCoarse(i) == 0)
1851 fineStart(i
1852 ) = (marker.getRelativePositionWithinFatherCell()(i)
1853 * numberOfDoFsPerAxisInPatch)
1854 % 3;
1855#if Dimensions == 3
1856 if (kCoarse(j) == 0)
1857 fineStart(j
1858 ) = (marker.getRelativePositionWithinFatherCell()(j)
1859 * numberOfDoFsPerAxisInPatch)
1860 % 3;
1861#endif
1862
1863 int modulus = (marker.getRelativePositionWithinFatherCell()(i) + 1)
1864 * numberOfDoFsPerAxisInPatch;
1865 if (modulus % 3 != 0)
1866 modulus = modulus % 3;
1867 else
1868 modulus = 3;
1869 if (kCoarse(i) == faceDimensions(i) - 1)
1870 fineEnd(i) = modulus;
1871
1872#if Dimensions == 3
1873 modulus = (marker.getRelativePositionWithinFatherCell()(j) + 1)
1874 * numberOfDoFsPerAxisInPatch;
1875 if (modulus % 3 != 0)
1876 modulus = modulus % 3;
1877 else
1878 modulus = 3;
1879 if (kCoarse(j) == faceDimensions(j) - 1)
1880 fineEnd(j) = modulus;
1881#endif
1882
1883 tarch::la::Vector<Dimensions, int> fineDimensions = fineEnd - fineStart;
1884 fineDimensions(normal) = std::min(overlap - kCoarse(normal) * 3, 3);
1885 if (fineDimensions(normal) % 3 != 0)
1886 fineDimensions(normal) = fineDimensions(normal) % 3;
1887 else
1888 fineDimensions(normal) = 3;
1889#if Dimensions == 3
1890 int fineDimensionsProduct = fineDimensions(normal) * fineDimensions(i)
1891 * fineDimensions(j);
1892#elif Dimensions == 2
1893 int fineDimensionsProduct = fineDimensions(normal) * fineDimensions(i);
1894#endif
1896 fineMatrix(fineDimensionsProduct, numberOfDerivatives);
1897 int row = 0;
1898
1899 dfor(kFine, fineDimensions) {
1901 cell = fineStart + kFine - tarch::la::Vector<Dimensions, int>(1);
1902 int column = 0;
1903 for (int a = 0; a < Dimensions; a++) {
1904 fineMatrix(row, column) = (1.0 / 3.0) * cell(a);
1905 column++;
1906 }
1907 for (int a = 0; a < Dimensions; a++) {
1908 for (int b = a; b < Dimensions; b++) {
1909 fineMatrix(row, column) = (1.0 / 18.0) * cell(a) * cell(b);
1910 column++;
1911 }
1912 }
1913 row++;
1914 }
1915 tarch::la::DynamicMatrix fineValues = fineMatrix * x;
1916
1917 int fineCellNumber = 0;
1919 fineCorner = kCoarse * 3 + fineStart - fineDataStart;
1920 fineCorner(normal
1921 ) = (marker.getSelectedFaceNumber() < Dimensions)
1922 ? (kCoarse(normal) * 3)
1923 : overlap + (kCoarse(normal) * 3);
1924
1925 dfor(kFine, fineDimensions) {
1926 tarch::la::Vector<Dimensions, int> fineCell = kFine + fineCorner;
1928 fineCell,
1929 numberOfDoFsPerAxisInPatch,
1930 overlap,
1931 normal
1932 );
1933
1934 for (int unknown = 0; unknown < unknowns; unknown++) {
1935 fineGridFaceValues[index * unknowns + unknown]
1936 = centerValues[unknown] + fineValues(fineCellNumber, unknown);
1937 }
1938 fineCellNumber++;
1939 }
1940 }
1941
1942 delete[] centerValues;
1943
1944}
1945
1948 int numberOfDoFsPerAxisInPatch,
1949 int overlap,
1950 int unknowns,
1951 const double* __restrict__ coarseGridFaceValues,
1952 double* __restrict__ fineGridFaceValues
1953) {
1954 const int normal = marker.getSelectedFaceNumber() % Dimensions;
1955 const int i = (normal + 1) % Dimensions;
1956#if Dimensions == 3
1957 const int j = (normal + 2) % Dimensions;
1958#endif
1959
1961 coarseStart(i
1962 ) = marker.getRelativePositionWithinFatherCell()(i)
1963 * numberOfDoFsPerAxisInPatch / 3;
1964#if Dimensions == 3
1965 coarseStart(j
1966 ) = marker.getRelativePositionWithinFatherCell()(j)
1967 * numberOfDoFsPerAxisInPatch / 3;
1968#endif
1969
1971 coarseEnd(i
1972 ) = ((marker.getRelativePositionWithinFatherCell()(i) + 1
1973 ) * numberOfDoFsPerAxisInPatch
1974 + 2)
1975 / 3;
1976#if Dimensions == 3
1977 coarseEnd(j
1978 ) = ((marker.getRelativePositionWithinFatherCell()(j) + 1
1979 ) * numberOfDoFsPerAxisInPatch
1980 + 2)
1981 / 3;
1982#endif
1983
1984 tarch::la::Vector<Dimensions, int> fineDataStart(0);
1985 fineDataStart(i
1986 ) = (marker.getRelativePositionWithinFatherCell()(i)
1987 * numberOfDoFsPerAxisInPatch)
1988 % 3;
1989#if Dimensions == 3
1990 fineDataStart(j
1991 ) = (marker.getRelativePositionWithinFatherCell()(j)
1992 * numberOfDoFsPerAxisInPatch)
1993 % 3;
1994#endif
1995
1996 tarch::la::Vector<Dimensions, int> faceDimensions = coarseEnd - coarseStart;
1997 faceDimensions(normal) = (overlap - 1) / 3 + 1;
1998
1999 double* centerValues = new double[unknowns];
2000 int indices[32];
2002
2003
2004 dfor(kCoarse, faceDimensions) {
2005 tarch::la::Vector<Dimensions, int> center = kCoarse + coarseStart;
2006 center(normal
2007 ) = (marker.getSelectedFaceNumber() < Dimensions)
2008 ? overlap - kCoarse(normal) - 1
2009 : overlap + kCoarse(normal);
2010 tarch::la::Vector<Dimensions, int> stencilOffset(0);
2011
2012 if (center(i) == 0)
2013 stencilOffset(i) = 1;
2014 if (center(i) == numberOfDoFsPerAxisInPatch - 1)
2015 stencilOffset(i) = -1;
2016#if Dimensions == 3
2017 if (center(j) == 0)
2018 stencilOffset(j) = 1;
2019 if (center(j) == numberOfDoFsPerAxisInPatch - 1)
2020 stencilOffset(j) = -1;
2021#endif
2022
2023 int stencilSize = 0;
2024 dfor(offset, 3) {
2025 if (offset != tarch::la::Vector<Dimensions, int>(1) - stencilOffset) {
2026 stencil[stencilSize] = stencilOffset + offset
2028 indices[stencilSize] = serialiseVoxelIndexInOverlap(
2029 center + offset + stencilOffset
2031 numberOfDoFsPerAxisInPatch,
2032 overlap,
2033 normal
2034 );
2035 stencilSize++;
2036 }
2037 }
2038
2039 for (int axis = 0; axis < Dimensions; axis++) {
2040 if (center(axis) > 1) {
2042 offset(axis) = -2;
2043 stencil[stencilSize] = stencilOffset + offset;
2044 indices[stencilSize] = serialiseVoxelIndexInOverlap(
2045 center + stencilOffset + offset,
2046 numberOfDoFsPerAxisInPatch,
2047 overlap,
2048 normal
2049 );
2050 stencilSize++;
2051 }
2052 if (center(axis) < numberOfDoFsPerAxisInPatch - 2) {
2054 offset(axis) = 2;
2055 stencil[stencilSize] = stencilOffset + offset;
2056 indices[stencilSize] = serialiseVoxelIndexInOverlap(
2057 center + stencilOffset + offset,
2058 numberOfDoFsPerAxisInPatch,
2059 overlap,
2060 normal
2061 );
2062 stencilSize++;
2063 }
2064 }
2065
2067 center,
2068 numberOfDoFsPerAxisInPatch,
2069 overlap,
2070 normal
2071 );
2072 for (int unknown = 0; unknown < unknowns; unknown++) {
2073 centerValues[unknown] = coarseGridFaceValues[index * unknowns + unknown];
2074 }
2075 tarch::la::DynamicMatrix deltas(stencilSize, unknowns);
2076 for (int k = 0; k < stencilSize; k++) {
2077 for (int unknown = 0; unknown < unknowns; unknown++) {
2078 deltas(
2079 k,
2080 unknown
2081 ) = coarseGridFaceValues[indices[k] * unknowns + unknown]
2082 - centerValues[unknown];
2083 }
2084 }
2085
2086 int numberOfDerivatives = (Dimensions == 2) ? 9 : 19;
2087 tarch::la::DynamicMatrix A(stencilSize, numberOfDerivatives);
2088
2089 for (int row = 0; row < stencilSize; row++) {
2090 int column = 0;
2092 for (int a = 0; a < Dimensions; a++) {
2093 A(row, column) = cell(a);
2094 column++;
2095 }
2096 for (int a = 0; a < Dimensions; a++) {
2097 for (int b = a; b < Dimensions; b++) {
2098 A(row, column) = 0.5 * cell(a) * cell(b);
2099 column++;
2100 }
2101 }
2102 for (int a = 0; a < Dimensions; a++) {
2103 for (int b = a; b < Dimensions; b++) {
2104 for (int c = b; c < Dimensions; c++) {
2105 A(row, column) = (1.0 / 6.0) * cell(a) * cell(b) * cell(c);
2106 column++;
2107 }
2108 }
2109 }
2110 }
2111 tarch::la::DynamicMatrix Q(A.rows(), A.cols());
2112 tarch::la::DynamicMatrix R(A.cols(), A.cols());
2114
2115 auto Q_T = tarch::la::transpose(Q);
2116 tarch::la::DynamicMatrix c = Q_T * deltas;
2117
2119 tarch::la::DynamicMatrix x = inverseR * c;
2122
2123
2124 if (kCoarse(i) == 0)
2125 fineStart(i
2126 ) = (marker.getRelativePositionWithinFatherCell()(i)
2127 * numberOfDoFsPerAxisInPatch)
2128 % 3;
2129#if Dimensions == 3
2130 if (kCoarse(j) == 0)
2131 fineStart(j
2132 ) = (marker.getRelativePositionWithinFatherCell()(j)
2133 * numberOfDoFsPerAxisInPatch)
2134 % 3;
2135#endif
2136
2137 int modulus = (marker.getRelativePositionWithinFatherCell()(i) + 1)
2138 * numberOfDoFsPerAxisInPatch;
2139 if (modulus % 3 != 0)
2140 modulus = modulus % 3;
2141 else
2142 modulus = 3;
2143 if (kCoarse(i) == faceDimensions(i) - 1)
2144 fineEnd(i) = modulus;
2145
2146#if Dimensions == 3
2147 modulus = (marker.getRelativePositionWithinFatherCell()(j) + 1)
2148 * numberOfDoFsPerAxisInPatch;
2149 if (modulus % 3 != 0)
2150 modulus = modulus % 3;
2151 else
2152 modulus = 3;
2153 if (kCoarse(j) == faceDimensions(j) - 1)
2154 fineEnd(j) = modulus;
2155#endif
2156
2157 tarch::la::Vector<Dimensions, int> fineDimensions = fineEnd - fineStart;
2158 fineDimensions(normal) = std::min(overlap - kCoarse(normal) * 3, 3);
2159 if (fineDimensions(normal) % 3 != 0)
2160 fineDimensions(normal) = fineDimensions(normal) % 3;
2161 else
2162 fineDimensions(normal) = 3;
2163
2164#if Dimensions == 3
2165 int fineDimensionsProduct = fineDimensions(normal) * fineDimensions(i)
2166 * fineDimensions(j);
2167#elif Dimensions == 2
2168 int fineDimensionsProduct = fineDimensions(normal) * fineDimensions(i);
2169#endif
2170
2172 fineMatrix(fineDimensionsProduct, numberOfDerivatives);
2173 int row = 0;
2174
2175 dfor(kFine, fineDimensions) {
2177 cell = fineStart + kFine - tarch::la::Vector<Dimensions, int>(1);
2178 int column = 0;
2179 for (int a = 0; a < Dimensions; a++) {
2180 fineMatrix(row, column) = (1.0 / 3.0) * cell(a);
2181 column++;
2182 }
2183 for (int a = 0; a < Dimensions; a++) {
2184 for (int b = a; b < Dimensions; b++) {
2185 fineMatrix(row, column) = (1.0 / 18.0) * cell(a) * cell(b);
2186 column++;
2187 }
2188 }
2189 for (int a = 0; a < Dimensions; a++) {
2190 for (int b = a; b < Dimensions; b++) {
2191 for (int c = b; c < Dimensions; c++) {
2192 fineMatrix(
2193 row,
2194 column
2195 ) = (1.0 / 162.0) * cell(a) * cell(b) * cell(c);
2196 column++;
2197 }
2198 }
2199 }
2200 row++;
2201 }
2202 tarch::la::DynamicMatrix fineValues = fineMatrix * x;
2203
2204 int fineCellNumber = 0;
2206 fineCorner = kCoarse * 3 + fineStart - fineDataStart;
2207 fineCorner(normal
2208 ) = (marker.getSelectedFaceNumber() < Dimensions)
2209 ? (kCoarse(normal) * 3)
2210 : overlap + (kCoarse(normal) * 3);
2211
2212 dfor(kFine, fineDimensions) {
2213 tarch::la::Vector<Dimensions, int> fineCell = kFine + fineCorner;
2215 fineCell,
2216 numberOfDoFsPerAxisInPatch,
2217 overlap,
2218 normal
2219 );
2220
2221 for (int unknown = 0; unknown < unknowns; unknown++) {
2222 fineGridFaceValues[index * unknowns + unknown]
2223 = centerValues[unknown] + fineValues(fineCellNumber, unknown);
2224 }
2225 fineCellNumber++;
2226 }
2227 }
2228
2229 delete[] centerValues;
2230}
2231
2232
2235 int numberOfDoFsPerAxisInPatch,
2236 int overlap,
2237 int unknowns,
2238 const double* __restrict__ interpolationData,
2239 const int* __restrict__ columnIndices,
2240 const int* __restrict__ rowIndices,
2241 int dataStart,
2242 int rowIndicesStart,
2243 const double* __restrict__ coarseGridFaceValues,
2244 double* __restrict__ fineGridFaceValues
2245) {
2246 tarch::la::Vector<Dimensions, int> fineFaceDimensions(
2247 numberOfDoFsPerAxisInPatch
2248 );
2249 tarch::la::Vector<Dimensions, int> coarseFaceDimensions(
2250 numberOfDoFsPerAxisInPatch
2251 );
2252
2253 const int normal = marker.getSelectedFaceNumber() % Dimensions;
2254 const int matrixColumns = 2 * overlap
2255 * std::
2256 pow(numberOfDoFsPerAxisInPatch, Dimensions - 1);
2257 const int matrixRows = overlap
2258 * std::pow(numberOfDoFsPerAxisInPatch, Dimensions - 1);
2259
2260 fineFaceDimensions(normal) = overlap;
2261 coarseFaceDimensions(normal) = 2 * overlap;
2262
2263 double* rearrangedCoarseData = new double[matrixColumns * unknowns];
2264
2265 dfor(kCoarse, coarseFaceDimensions) {
2266 int currentIndex = serialiseVoxelIndexInOverlap(
2267 kCoarse,
2268 numberOfDoFsPerAxisInPatch,
2269 overlap,
2270 normal
2271 );
2272 int newIndex = (marker.getSelectedFaceNumber() >= Dimensions)
2273 ? kCoarse(normal)
2274 : 2 * overlap - kCoarse(normal) - 1;
2275 int base = overlap * 2;
2276 for (int d = 1; d < Dimensions; d++) {
2277 newIndex += kCoarse((normal + d) % Dimensions) * base;
2278 base *= numberOfDoFsPerAxisInPatch;
2279 }
2280 for (int unknown = 0; unknown < unknowns; unknown++) {
2281 rearrangedCoarseData[newIndex * unknowns + unknown] = coarseGridFaceValues
2282 [currentIndex * unknowns + unknown];
2283 }
2284 }
2285
2286 dfor(kFine, fineFaceDimensions) {
2287 auto fineCoords = kFine;
2288 fineCoords(normal) = (marker.getSelectedFaceNumber() < Dimensions)
2289 ? overlap - kFine(normal) - 1
2290 : overlap + kFine(normal);
2291
2292 int fineIndex = serialiseVoxelIndexInOverlap(
2293 fineCoords,
2294 numberOfDoFsPerAxisInPatch,
2295 overlap,
2296 normal
2297 );
2298
2299 int row = kFine(normal);
2300 int base = overlap;
2301 for (int d = 1; d < Dimensions; d++) {
2302 row += kFine((normal + d) % Dimensions) * base;
2303 base *= numberOfDoFsPerAxisInPatch;
2304 }
2305
2306 for (int unknown = 0; unknown < unknowns; unknown++) {
2307 fineGridFaceValues[fineIndex * unknowns + unknown] = 0.0;
2308 }
2309
2310 int columnStart = rowIndices[rowIndicesStart + row];
2311 int columnEnd = rowIndices[rowIndicesStart + row + 1];
2312
2313 for (int column = columnStart; column < columnEnd; column++) {
2314 for (int unknown = 0; unknown < unknowns; unknown++) {
2315 fineGridFaceValues[fineIndex * unknowns + unknown]
2316 += interpolationData[dataStart + column]
2317 * rearrangedCoarseData
2318 [(columnIndices[dataStart + column]) * unknowns + unknown];
2319 }
2320 }
2321 }
2322
2323 delete[] rearrangedCoarseData;
2324}
2325
2328 int numberOfDoFsPerAxisInPatch,
2329 int overlap,
2330 int unknowns,
2331 const double* __restrict__ normalInterpolationMatrix1d,
2332 const double* __restrict__ tangentialInterpolationMatrix1d,
2333 const double* __restrict__ coarseGridFaceValues,
2334 double* __restrict__ fineGridFaceValues
2335) {
2337 "interpolateHaloLayer_AoS_tensor_product(...)",
2338 marker.toString(),
2339 numberOfDoFsPerAxisInPatch,
2340 overlap
2341 );
2342
2343 const int normal = marker.getSelectedFaceNumber() % Dimensions;
2344
2345 dfore(kFine, numberOfDoFsPerAxisInPatch, normal, 0) {
2346 for (int iFine = 0; iFine < overlap; iFine++) {
2348 dest(normal
2349 ) = (marker.getSelectedFaceNumber() < Dimensions)
2350 ? iFine
2351 : 2 * overlap - iFine - 1;
2352
2353 int destLinearised = serialiseVoxelIndexInOverlap(
2354 dest,
2355 numberOfDoFsPerAxisInPatch,
2356 overlap,
2357 normal
2358 );
2359
2360 logDebug(
2361 "interpolateHaloLayer_AoS_tensor_product(...)",
2362 "clear dof " << dest
2363 );
2364 for (int unknown = 0; unknown < unknowns; unknown++) {
2365 fineGridFaceValues[destLinearised * unknowns + unknown] = 0.0;
2366 }
2367
2368 dfore(kCoarse, numberOfDoFsPerAxisInPatch, normal, 0) {
2369 for (int iCoarse = 0; iCoarse < 2 * overlap; iCoarse++) {
2371 src(normal
2372 ) = (marker.getSelectedFaceNumber() < Dimensions)
2373 ? iCoarse
2374 : 2 * overlap - iCoarse - 1;
2375 int srcLinearised = serialiseVoxelIndexInOverlap(
2376 src,
2377 numberOfDoFsPerAxisInPatch,
2378 overlap,
2379 normal
2380 );
2381
2382 double weight = 1.0;
2383 for (int d = 0; d < Dimensions; d++) {
2384 if (d == normal) {
2385 int col = iCoarse;
2386 int row = iFine;
2387
2388 assertion4(col >= 0, row, col, src, dest);
2389 assertion4(row >= 0, row, col, src, dest);
2390 int index = col + row * 2 * overlap;
2391 logDebug(
2392 "interpolateHaloLayer_AoS_tensor_product(...)",
2393 "(" << row << "," << col << ")=" << index
2394 << " (normal contribution)"
2395 );
2396 weight *= normalInterpolationMatrix1d[index];
2397 assertion(weight == weight);
2398 } else {
2399 int col = src(d);
2400 int row = dest(d)
2401 + marker.getRelativePositionWithinFatherCell()(d
2402 ) * numberOfDoFsPerAxisInPatch;
2403 assertion4(col >= 0, row, col, src, dest);
2404 assertion4(row >= 0, row, col, src, dest);
2405 int index = col + row * numberOfDoFsPerAxisInPatch;
2406 logDebug(
2407 "interpolateHaloLayer_AoS_tensor_product(...)",
2408 "(" << row << "," << col << ")=" << index
2409 << " (tangential contribution)"
2410 );
2411 weight *= tangentialInterpolationMatrix1d[index];
2412 assertion(weight == weight);
2413 }
2414 }
2415
2416 logDebug(
2417 "interpolateHaloLayer_AoS_tensor_product(...)",
2418 "add dof " << src << " to " << dest << " with weight " << weight
2419 );
2420
2421 for (int unknown = 0; unknown < unknowns; unknown++) {
2422 fineGridFaceValues[destLinearised * unknowns + unknown]
2423 += weight
2424 * coarseGridFaceValues[srcLinearised * unknowns + unknown];
2425 }
2426 }
2427 }
2428 }
2429 }
2430
2431 logTraceOut("interpolateHaloLayer_AoS_tensor_product(...)");
2432}
2433
2435 [[maybe_unused]] const peano4::datamanagement::FaceMarker& marker,
2436 [[maybe_unused]] int numberOfDoFsPerAxisInPatch,
2437 [[maybe_unused]] int overlap,
2438 [[maybe_unused]] int unknowns,
2439 [[maybe_unused]] const double* __restrict__ normalInterpolationMatrix1d,
2440 [[maybe_unused]] const double* __restrict__ tangentialInterpolationMatrix1d,
2441 [[maybe_unused]] const double* __restrict__ coarseGridCellValues,
2442 [[maybe_unused]] const double* __restrict__ coarseGridFaceValues,
2443 [[maybe_unused]] double* __restrict__ fineGridFaceValues
2444) {
2445 assertion(false);
2446}
2447
2449 [[maybe_unused]] const peano4::datamanagement::FaceMarker& marker,
2450 [[maybe_unused]] int numberOfDoFsPerAxisInPatch,
2451 [[maybe_unused]] int overlap,
2452 [[maybe_unused]] int unknowns,
2453 [[maybe_unused]] const double* __restrict__ interpolationData,
2454 [[maybe_unused]] const int* __restrict__ columnIndices,
2455 [[maybe_unused]] const int* __restrict__ rowIndices,
2456 [[maybe_unused]] const double* __restrict__ coarseGridCellValues,
2457 [[maybe_unused]] const double* __restrict__ coarseGridFaceValues,
2458 [[maybe_unused]] double* __restrict__ fineGridFaceValues
2459) {
2460 assertion(false);
2461}
2462
2464 [[maybe_unused]] const peano4::datamanagement::FaceMarker& marker,
2465 [[maybe_unused]] int numberOfDoFsPerAxisInPatch,
2466 [[maybe_unused]] int overlap,
2467 [[maybe_unused]] int unknowns,
2468 [[maybe_unused]] const double* __restrict__ coarseGridCellValues,
2469 [[maybe_unused]] const double* __restrict__ coarseGridFaceValues,
2470 [[maybe_unused]] double* __restrict__ fineGridFaceValues
2471) {
2472 assertion(false);
2473}
2474
2476 [[maybe_unused]] const peano4::datamanagement::FaceMarker& marker,
2477 [[maybe_unused]] int numberOfDoFsPerAxisInPatch,
2478 [[maybe_unused]] int overlap,
2479 [[maybe_unused]] int unknowns,
2480 [[maybe_unused]] const double* __restrict__ coarseGridCellValues,
2481 [[maybe_unused]] const double* __restrict__ coarseGridFaceValues,
2482 [[maybe_unused]] double* __restrict__ fineGridFaceValues
2483) {
2484 assertion(false);
2485}
2486
2488 [[maybe_unused]] const peano4::datamanagement::CellMarker& marker,
2489 [[maybe_unused]] int numberOfDoFsPerAxisInPatch,
2490 [[maybe_unused]] int unknowns,
2491 [[maybe_unused]] const double* __restrict__ interpolationMatrix1d,
2492 [[maybe_unused]] const double* __restrict__ coarseGridCellValues,
2493 [[maybe_unused]] double* __restrict__ fineGridCellValues
2494) {
2495 assertion(false);
2496}
2497
2499 [[maybe_unused]] const peano4::datamanagement::CellMarker& marker,
2500 [[maybe_unused]] int numberOfDoFsPerAxisInPatch,
2501 [[maybe_unused]] int unknowns,
2502 [[maybe_unused]] const double* __restrict__ interpolationData,
2503 [[maybe_unused]] const int* __restrict__ columnIndices,
2504 [[maybe_unused]] const int* __restrict__ rowIndices,
2505 [[maybe_unused]] const double* __restrict__ coarseGridCellValues,
2506 [[maybe_unused]] double* __restrict__ fineGridCellValues
2507) {
2508 assertion(false);
2509}
2510
2512 [[maybe_unused]] const peano4::datamanagement::CellMarker& marker,
2513 [[maybe_unused]] int numberOfDoFsPerAxisInPatch,
2514 [[maybe_unused]] int unknowns,
2515 [[maybe_unused]] const double* __restrict__ coarseGridCellValues,
2516 [[maybe_unused]] double* __restrict__ fineGridCellValues
2517) {
2518 assertion(false);
2519}
2520
2522 [[maybe_unused]] const peano4::datamanagement::CellMarker& marker,
2523 [[maybe_unused]] int numberOfDoFsPerAxisInPatch,
2524 [[maybe_unused]] int unknowns,
2525 [[maybe_unused]] const double* __restrict__ coarseGridCellValues,
2526 [[maybe_unused]] double* __restrict__ fineGridCellValues
2527) {
2528 assertion(false);
2529}
#define assertion4(expr, param0, param1, param2, param3)
#define assertionEquals1(lhs, rhs, a)
#define assertionEquals(lhs, rhs)
#define assertion(expr)
#define logDebug(methodName, logMacroMessageStream)
Definition Log.h:50
#define logTraceOut(methodName)
Definition Log.h:379
#define logTraceInWith4Arguments(methodName, argument0, argument1, argument2, argument3)
Definition Log.h:373
#define logTraceInWith3Arguments(methodName, argument0, argument1, argument2)
Definition Log.h:372
#define logTraceInWith1Argument(methodName, argument0)
Definition Log.h:370
#define dfor3(counter)
Shortcut For dfor(counter,3)
Definition Loop.h:647
#define dfore(counter, max, dim, value)
This is an exclusive d-dimensional for loop.
Definition Loop.h:942
#define dfor(counter, max)
d-dimensional Loop
Definition Loop.h:313
#define enddforx
I prefer to use this macro for dforx instead of a closing bracket as many syntax parser fail otherwis...
Definition Loop.h:933
tarch::logging::Log _log("::")
float P
Definition vec.h:7
My standard matrix is a matrix where the size is fixed at compile time.
void batchedMultiplyAoS(double *__restrict__ result, const double *__restrict__ x, int batchCount, int resultSize, int firstRow)
This operation assumes that x holds a whole batch of vectors in AoS format.
void multiplyBySmallMatrix(double *result, const DynamicMatrix &matrix) const
void replicateRows(int blockSize, int numberOfReplications, int shiftAfterEveryReplication, bool extendColumnsToAccommodateShifts)
Split the matrix into blocks of rows of size blockSize.
std::string toString(bool addLineBreaks=false) const
Log Device.
Definition Log.h:516
Create a lock around a boolean semaphore region.
Definition Lock.h:19
void free()
Free the lock.
Definition Lock.cpp:37
#define endSimtDfor
Definition Loop.h:181
const float const float const float struct part *restrict struct part *restrict const float a
Definition hydro_iact.h:55
list offset
Definition acoustic.py:42
j
Definition euler.py:99
index
Definition makeIC.py:38
LoopPlacement
Guide loop-level parallelism.
Definition Loop.h:65
CPUGPUMethod int dLinearised(const tarch::la::Vector< Dimensions, int > &counter, int max)
Map d-dimensional vector onto integer index.
Definition Loop.cpp:106
list stencil
Definition poisson.py:32
CF abs(const CF &cf)
b
Definition swe.py:86
My collection of tiny vector operations.
DynamicMatrix transpose(const DynamicMatrix &matrix)
bool greaterEquals(double lhs, double rhs, double tolerance=NUMERICAL_ZERO_DIFFERENCE)
DynamicMatrix invertUpperTriangular(const DynamicMatrix &R)
void modifiedGramSchmidt(Matrix< Rows, Cols, Scalar > A, Matrix< Rows, Cols, Scalar > &Q, Matrix< Cols, Cols, Scalar > &R)
Produces an economy-size QR decomposition of a matrix A, A is changed.
bool smallerEquals(double lhs, double rhs, double tolerance=NUMERICAL_ZERO_DIFFERENCE)
tarch::la::Vector< Size, NewScalarType > convertScalar(const tarch::la::Vector< Size, Scalar > &vector)
Definition Vector.cpph:94
Matrix< Rows, Cols, Scalar > multiplyComponents(const Matrix< Rows, X, Scalar > &lhs, const Matrix< X, Cols, Scalar > &rhs)
Vector< Cols, Scalar > row(const Matrix< Rows, Cols, Scalar > &matrix, int whichRow)
Extract row from matrix.
void freeMemory(void *data, MemoryLocation location, int device=accelerator::Device::HostDevice)
Free memory.
@ Heap
Create data on the heap of the local device.
T * allocateMemory(std::size_t count, MemoryLocation location, int device=accelerator::Device::HostDevice)
Allocates memory on the specified memory location.
Definition accelerator.h:99
tarch::la::DynamicMatrix * createInterpolationMatrixFrom1dTemplateByInsertingZeroColsAndRows(const tarch::la::DynamicMatrix &P1d, int numberOfDoFsPerAxisInPatch, int normal)
This routine accepts a matrix as input.
void interpolateCell_AoS_linear(const tarch::la::Vector< Dimensions, int > &relativePositionWithinFatherCell, int numberOfDoFsPerAxisInPatch, int unknowns, const double *__restrict__ fineGridValues, double *coarseGridValues, bool extrapolateLinearly)
std::map< FaceInterpolationOperatorKey, tarch::la::DynamicMatrix * > FaceInterpolationMap
tarch::la::DynamicMatrix * createInterpolationMatrixFrom1dTemplateByLinearInterpolationAlongNormal(const tarch::la::DynamicMatrix &P1d, int numberOfDoFsPerAxisInPatch, int normal)
This is an extension of createInterpolationMatrixFrom1dTemplateByInsertingZeroColsAndRows().
void interpolateCell_AoS_piecewise_constant(const tarch::la::Vector< Dimensions, int > &relativePositionWithinFatherCell, int numberOfDoFsPerAxisInPatch, int unknowns, const double *__restrict__ fineGridValues, double *coarseGridValues)
void projectInterpolatedFineCellsOnHaloLayer_AoS(const peano4::datamanagement::FaceMarker &marker, int numberOfDoFsPerAxisInPatch, int overlap, int unknowns, const double *__restrict__ fineGridCellValuesLeft, const double *__restrict__ fineGridCellValuesRight, double *fineGridFaceValues)
This routine assumes that you have the whole patch data of hte left and right adjacent patch at hands...
tarch::la::DynamicMatrix * createPiecewiseConstantInterpolationMatrix(int numberOfDoFsPerAxisInPatch, int normal, int overlap)
Create piecewise constant interpolation matrix.
tarch::la::DynamicMatrix * createLinearInterpolationMatrix(int numberOfDoFsPerAxisInPatch, int normal, bool extrapolateLinearly, bool interpolateLinearlyAlongNormal)
The realisation relies on the following observations/follows these steps:
void interpolateHaloLayer_AoS_tensor_product(const peano4::datamanagement::FaceMarker &marker, int numberOfDoFsPerAxisInPatch, int overlap, int unknowns, const double *__restrict__ coarseGridFaceValues, double *__restrict__ fineGridFaceValues)
This is a wrapper around the toolbox routines.
void interpolateCell_AoS_matrix(const peano4::datamanagement::CellMarker &marker, int numberOfDoFsPerAxisInPatch, int unknowns, const double *__restrict__ coarseGridCellValues, double *__restrict__ fineGridCellValues)
int serialisePatchIndexInOverlap(const tarch::la::Vector< Dimensions, int > &patchIndex, int normal)
Patches along a face are basically organised as arrays, i.e.
void interpolateCellDataAssociatedToVolumesIntoOverlappingCell_linear(int numberOfDoFsPerAxisInSourcePatch, int numberOfDoFsPerAxisInDestinationPatch, int haloSourcePatch, int haloDestinationPatch, int unknowns, const double *__restrict__ sourceValues, double *__restrict__ destinationValues, ::peano4::utils::LoopPlacement parallelisation)
void interpolateHaloLayer_AoS_matrix(const peano4::datamanagement::FaceMarker &marker, int numberOfDoFsPerAxisInPatch, int overlap, int unknowns, const double *__restrict__ coarseGridCellValues, const double *__restrict__ coarseGridFaceValues, double *__restrict__ fineGridFaceValues)
void interpolateCell_AoS_linear_with_constant_extrapolation(const peano4::datamanagement::CellMarker &marker, int numberOfDoFsPerAxisInPatch, int unknowns, const double *__restrict__ coarseGridCellValues, double *__restrict__ fineGridCellValues)
void interpolateCell_AoS_piecewise_constant(const peano4::datamanagement::CellMarker &marker, int numberOfDoFsPerAxisInPatch, int unknowns, const double *__restrict__ coarseGridCellValues, double *__restrict__ fineGridCellValues)
This routine is called if we create a new cell (dynamic AMR)
void interpolateCell_AoS_linear_with_linear_extrapolation(const peano4::datamanagement::CellMarker &marker, int numberOfDoFsPerAxisInPatch, int unknowns, const double *__restrict__ coarseGridCellValues, double *__restrict__ fineGridCellValues)
void interpolateCellDataAssociatedToVolumesIntoOverlappingCell_secondOrder(int numberOfDoFsPerAxisInSourcePatch, int numberOfDoFsPerAxisInDestinationPatch, int haloSourcePatch, int haloDestinationPatch, int unknowns, const double *__restrict__ sourceValues, double *__restrict__ destinationValues, ::peano4::utils::LoopPlacement parallelisation)
void interpolateCell_AoS_linear_with_linear_extrapolation_and_linear_normal_interpolation(const peano4::datamanagement::CellMarker &marker, int numberOfDoFsPerAxisInPatch, int unknowns, const double *__restrict__ coarseGridCellValues, double *__restrict__ fineGridCellValues)
void interpolateHaloLayer_AoS_linear_with_constant_extrapolation(const peano4::datamanagement::FaceMarker &marker, int numberOfDoFsPerAxisInPatch, int overlap, int unknowns, const double *__restrict__ coarseGridFaceValues, double *__restrict__ fineGridFaceValues)
int serialiseVoxelIndexInOverlap(const tarch::la::Vector< Dimensions, int > &overlapCell, int numberOfDoFsPerAxisInPatch, int overlap, int normal)
The volumes or elements within an overlap are always enumerated lexicographically.
void interpolateHaloLayer_AoS_piecewise_constant(const peano4::datamanagement::FaceMarker &marker, int numberOfDoFsPerAxisInPatch, int overlap, int unknowns, const double *__restrict__ coarseGridFaceValues, double *__restrict__ fineGridFaceValues)
Take the coarse grid values and interpolate them onto the fine grid.
void interpolateCell_AoS_tensor_product(const peano4::datamanagement::CellMarker &marker, int numberOfDoFsPerAxisInPatch, int unknowns, const double *__restrict__ coarseGridCellValues, double *__restrict__ fineGridCellValues)
void interpolateHaloLayer_AoS_linear_with_linear_extrapolation_and_linear_normal_interpolation(const peano4::datamanagement::FaceMarker &marker, int numberOfDoFsPerAxisInPatch, int overlap, int unknowns, const double *__restrict__ coarseGridFaceValues, double *__restrict__ fineGridFaceValues)
void interpolateHaloLayer_AoS_third_order(const peano4::datamanagement::FaceMarker &marker, int numberOfDoFsPerAxisInPatch, int overlap, int unknowns, const double *__restrict__ coarseGridCellValues, const double *__restrict__ coarseGridFaceValues, double *__restrict__ fineGridFaceValues)
void interpolateHaloLayer_AoS_linear_with_linear_extrapolation(const peano4::datamanagement::FaceMarker &marker, int numberOfDoFsPerAxisInPatch, int overlap, int unknowns, const double *__restrict__ coarseGridFaceValues, double *__restrict__ fineGridFaceValues)
void interpolateCellDataAssociatedToVolumesIntoOverlappingCell_fourthOrder(int numberOfDoFsPerAxisInSourcePatch, int numberOfDoFsPerAxisInDestinationPatch, int haloSourcePatch, int haloDestinationPatch, int unknowns, const double *__restrict__ sourceValues, double *__restrict__ destinationValues, ::peano4::utils::LoopPlacement parallelisation)
void interpolateHaloLayer_AoS_second_order(const peano4::datamanagement::FaceMarker &marker, int numberOfDoFsPerAxisInPatch, int overlap, int unknowns, const double *__restrict__ coarseGridCellValues, const double *__restrict__ coarseGridFaceValues, double *__restrict__ fineGridFaceValues)
void interpolateHaloLayer_AoS_linear_with_constant_extrapolation_and_linear_normal_interpolation(const peano4::datamanagement::FaceMarker &marker, int numberOfDoFsPerAxisInPatch, int overlap, int unknowns, const double *__restrict__ coarseGridFaceValues, double *__restrict__ fineGridFaceValues)
void interpolateCell_AoS_third_order(const peano4::datamanagement::CellMarker &marker, int numberOfDoFsPerAxisInPatch, int unknowns, const double *__restrict__ coarseGridCellValues, double *__restrict__ fineGridCellValues)
void interpolateCell_AoS_linear_with_constant_extrapolation_and_linear_normal_interpolation(const peano4::datamanagement::CellMarker &marker, int numberOfDoFsPerAxisInPatch, int unknowns, const double *__restrict__ coarseGridCellValues, double *__restrict__ fineGridCellValues)
void interpolateCellDataAssociatedToVolumesIntoOverlappingCell_matrix(int numberOfDoFsPerAxisInSourcePatch, int numberOfDoFsPerAxisInDestinationPatch, int haloSourcePatch, int haloDestinationPatch, int unknowns, const double *__restrict__ interpolationData, const int *__restrict__ columnIndices, const int *__restrict__ rowIndices, const double *__restrict__ sourceValues, double *__restrict__ destinationValues, ::peano4::utils::LoopPlacement parallelisation)
This interpolation should be used if a cell hosts two sets of unknowns.
void interpolateCell_AoS_second_order(const peano4::datamanagement::CellMarker &marker, int numberOfDoFsPerAxisInPatch, int unknowns, const double *__restrict__ coarseGridCellValues, double *__restrict__ fineGridCellValues)
Provide information about selected face.
Definition FaceMarker.h:35
Simple vector class.
Definition Vector.h:134