23 const int patchSize = numberOfDoFsPerAxisInPatch * numberOfDoFsPerAxisInPatch;
24 const int numberOfPatches = 3 * 3;
26 const int patchSize = numberOfDoFsPerAxisInPatch * numberOfDoFsPerAxisInPatch
27 * numberOfDoFsPerAxisInPatch;
28 const int numberOfPatches = 3 * 3 * 3;
32 patchSize * numberOfPatches,
37 dfor3(patchIndex)
dfor(volumeIndex, numberOfDoFsPerAxisInPatch) {
39 sourceCell = (volumeIndex + patchIndex * numberOfDoFsPerAxisInPatch) / 3;
42 numberOfDoFsPerAxisInPatch
44 (*result)(currentRow, sourceCellLinearised) = 1.0;
50 "createPiecewiseConstantInterpolationMatrix(int)",
51 "created new volumetric matrix for "
52 << numberOfDoFsPerAxisInPatch <<
": " << result->
toString()
59 [[maybe_unused]]
int numberOfDoFsPerAxisInPatch,
60 [[maybe_unused]]
int normal,
61 [[maybe_unused]]
int overlap
68 numberOfDoFsPerAxisInPatch,
76 int numberOfDoFsPerAxisInPatch,
94 pattern = numberOfDoFsPerAxisInPatch;
97 pattern = numberOfDoFsPerAxisInPatch * numberOfDoFsPerAxisInPatch;
101 P->insertEmptyColumns(pattern, pattern, pattern);
102 P->replicateRows(pattern, 2, pattern,
false);
104 "createInterpolationMatrixFrom1dTemplateByInsertingZeroColsAndRows(...)",
105 "matrix for normal=" << normal <<
": " <<
P->toString()
113 int numberOfDoFsPerAxisInPatch,
131 pattern = numberOfDoFsPerAxisInPatch;
134 pattern = numberOfDoFsPerAxisInPatch * numberOfDoFsPerAxisInPatch;
143 P->replicateCols(pattern, 2, 0,
false);
144 P->replicateRows(pattern, 2, 0,
false);
146 "createInterpolationMatrixFrom1dTemplateByLinearInterpolationAlongNormal(...)",
147 "matrix for normal=" << normal <<
": " <<
P->toString()
154 int numberOfDoFsPerAxisInPatch,
156 bool extrapolateLinearly,
157 bool interpolateLinearlyBetweenRealAndRestrictedVolumes
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}}
168 P1d.removeColumn(numberOfDoFsPerAxisInPatch);
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;
174 numberOfDoFsPerAxisInPatch * 3 - 1,
175 numberOfDoFsPerAxisInPatch - 1
176 ) = 2.0 / 3.0 + 1.0 / 3.0 * 2.0;
178 numberOfDoFsPerAxisInPatch * 3 - 1,
179 numberOfDoFsPerAxisInPatch - 2
185 numberOfDoFsPerAxisInPatch * 3 - 1,
186 numberOfDoFsPerAxisInPatch - 1
189 numberOfDoFsPerAxisInPatch * 3 - 1,
190 numberOfDoFsPerAxisInPatch - 2
195 "createLinearInterpolationMatrix(...)",
196 "1d matrix: " << P1d.toString()
199 if (interpolateLinearlyBetweenRealAndRestrictedVolumes) {
202 numberOfDoFsPerAxisInPatch,
208 numberOfDoFsPerAxisInPatch,
216 int numberOfDoFsPerAxisInPatch,
217 bool extrapolateLinearly
220 "createLinearInterpolationMatrix(...)",
221 numberOfDoFsPerAxisInPatch
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}}
232 P1d.removeColumn(numberOfDoFsPerAxisInPatch);
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;
238 numberOfDoFsPerAxisInPatch * 3 - 1,
239 numberOfDoFsPerAxisInPatch - 1
240 ) = 2.0 / 3.0 + 1.0 / 3.0 * 2.0;
242 numberOfDoFsPerAxisInPatch * 3 - 1,
243 numberOfDoFsPerAxisInPatch - 2
249 numberOfDoFsPerAxisInPatch * 3 - 1,
250 numberOfDoFsPerAxisInPatch - 1
253 numberOfDoFsPerAxisInPatch * 3 - 1,
254 numberOfDoFsPerAxisInPatch - 2
259 "createLinearInterpolationMatrix(...)",
260 "1d interpolation matrix: " << P1d.toString(
true)
265 "createLinearInterpolationMatrix(...)",
266 "2d interpolation matrix: " << P2d.
toString(
true)
269 logTraceOut(
"createLinearInterpolationMatrix(...)");
280 int numberOfDoFsPerAxisInPatch,
283 const double* __restrict__ fineGridCellValuesLeft,
284 const double* __restrict__ fineGridCellValuesRight,
285 double* fineGridFaceValues
288 "projectInterpolatedFineCellsOnHaloLayer_AoS(...)",
290 numberOfDoFsPerAxisInPatch,
294 const int normal = marker.getSelectedFaceNumber() % Dimensions;
296 dfore(kFine, numberOfDoFsPerAxisInPatch, normal, 0) {
297 for (
int iFine = 0; iFine < overlap; iFine++) {
301 fineVolumeLeftDest(normal) = iFine;
302 fineVolumeLeftSrc(normal) = numberOfDoFsPerAxisInPatch - overlap + iFine;
307 fineVolumeRightDest(normal) = iFine + overlap;
308 fineVolumeRightSrc(normal) = iFine;
312 numberOfDoFsPerAxisInPatch,
318 numberOfDoFsPerAxisInPatch,
324 numberOfDoFsPerAxisInPatch
328 numberOfDoFsPerAxisInPatch
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];
341 logTraceOut(
"projectInterpolatedFineCellsOnHaloLayer_AoS(...)");
351 int numberOfDoFsPerAxisInPatch,
354 const double* __restrict__ coarseGridFaceValues,
355 double* __restrict__ fineGridFaceValues
358 "interpolateHaloLayer_AoS_piecewise_constant(...)",
360 numberOfDoFsPerAxisInPatch,
365 const int normal = marker.getSelectedFaceNumber() % Dimensions;
369 key(numberOfDoFsPerAxisInPatch, overlap, normal);
372 if (faceInterpolationMap.count(key) == 0) {
373 faceInterpolationMap[key] = internal::
374 createPiecewiseConstantInterpolationMatrix(
375 numberOfDoFsPerAxisInPatch,
384 const int patchSize = numberOfDoFsPerAxisInPatch * 2 * overlap;
386 const int patchSize = numberOfDoFsPerAxisInPatch * numberOfDoFsPerAxisInPatch
391 "interpolateHaloLayer_AoS_piecewise_constant(...)",
392 marker.toString() <<
": " <<
P->toString()
395 P->batchedMultiplyAoS(
397 coarseGridFaceValues,
401 marker.getRelativePositionWithinFatherCell(),
406 logTraceOut(
"interpolateHaloLayer_AoS_piecewise_constant(...)");
411 int numberOfDoFsPerAxisInPatch,
414 const double* __restrict__ coarseGridCellValues,
415 const double* __restrict__ coarseGridFaceValues,
416 double* __restrict__ fineGridFaceValues
419 "interpolateHaloLayer_AoS_piecewise_constant(...)",
421 numberOfDoFsPerAxisInPatch,
426 if (marker.isInteriorFaceWithinPatch()) {
428 const int patchSize = numberOfDoFsPerAxisInPatch
429 * numberOfDoFsPerAxisInPatch;
431 const int patchSize = numberOfDoFsPerAxisInPatch * numberOfDoFsPerAxisInPatch
432 * numberOfDoFsPerAxisInPatch;
435 double* leftAdjacentPatchData = tarch::allocateMemory<double>(
436 patchSize * unknowns,
439 double* rightAdjacentPatchData = tarch::allocateMemory<double>(
440 patchSize * unknowns,
445 leftAdjacentPatchIndex = marker.getRelativePositionWithinFatherCell();
447 rightAdjacentPatchIndex = marker.getRelativePositionWithinFatherCell();
448 leftAdjacentPatchIndex(marker.getSelectedFaceNumber() % Dimensions)--;
450 internal::interpolateCell_AoS_piecewise_constant(
451 leftAdjacentPatchIndex,
452 numberOfDoFsPerAxisInPatch,
454 coarseGridCellValues,
455 leftAdjacentPatchData
457 internal::interpolateCell_AoS_piecewise_constant(
458 rightAdjacentPatchIndex,
459 numberOfDoFsPerAxisInPatch,
461 coarseGridCellValues,
462 rightAdjacentPatchData
465 internal::projectInterpolatedFineCellsOnHaloLayer_AoS(
467 numberOfDoFsPerAxisInPatch,
470 leftAdjacentPatchData,
471 rightAdjacentPatchData,
480 numberOfDoFsPerAxisInPatch,
483 coarseGridFaceValues,
488 logTraceOut(
"interpolateHaloLayer_AoS_piecewise_constant(...)");
493 int numberOfDoFsPerAxisInPatch,
495 const double* __restrict__ coarseGridCellValues,
496 double* __restrict__ fineGridCellValues
499 "interpolateCell_AoS_piecewise_constant(...)",
501 numberOfDoFsPerAxisInPatch,
505 internal::interpolateCell_AoS_piecewise_constant(
506 marker.getRelativePositionWithinFatherCell(),
507 numberOfDoFsPerAxisInPatch,
509 coarseGridCellValues,
513 logTraceOut(
"interpolateCell_AoS_piecewise_constant(...)");
518 int numberOfDoFsPerAxisInPatch,
520 const double* __restrict__ coarseGridCellValues,
521 double* __restrict__ fineGridCellValues
523 static internal::CellInterpolationMap cellInterpolationMap;
524 internal::CellInterpolationOperatorKey key(numberOfDoFsPerAxisInPatch);
527 if (cellInterpolationMap.count(key) == 0) {
528 cellInterpolationMap[key] = internal::
529 createPiecewiseConstantInterpolationMatrix(numberOfDoFsPerAxisInPatch);
535 const int patchSize = numberOfDoFsPerAxisInPatch * numberOfDoFsPerAxisInPatch;
537 const int patchSize = numberOfDoFsPerAxisInPatch * numberOfDoFsPerAxisInPatch
538 * numberOfDoFsPerAxisInPatch;
543 coarseGridCellValues,
546 serialiseMarkerIn3x3PatchAssembly(
547 relativePositionWithinFatherCell,
548 numberOfDoFsPerAxisInPatch
561 int numberOfDoFsPerAxisInPatch,
564 const double* __restrict__ coarseGridFaceValues,
565 double* __restrict__ fineGridFaceValues
568 "interpolateHaloLayer_AoS_linear_with_constant_extrapolation(...)",
570 numberOfDoFsPerAxisInPatch,
575 const int normal = marker.getSelectedFaceNumber() % Dimensions;
579 key(numberOfDoFsPerAxisInPatch, overlap, normal);
584 if (faceInterpolationMap.count(key) == 0) {
585 faceInterpolationMap[key] = internal::createLinearInterpolationMatrix(
586 numberOfDoFsPerAxisInPatch,
596 const int patchSize = numberOfDoFsPerAxisInPatch * 2 * overlap;
598 const int patchSize = numberOfDoFsPerAxisInPatch * numberOfDoFsPerAxisInPatch
603 "interpolateHaloLayer_AoS_linear_with_constant_extrapolation(...)",
604 marker.toString() <<
": " <<
P->toString()
607 P->batchedMultiplyAoS(
609 coarseGridFaceValues,
613 marker.getRelativePositionWithinFatherCell(),
618 logTraceOut(
"interpolateHaloLayer_AoS_linear_with_constant_extrapolation(...)"
625 int numberOfDoFsPerAxisInPatch,
628 const double* __restrict__ coarseGridCellValues,
629 const double* __restrict__ coarseGridFaceValues,
630 double* __restrict__ fineGridFaceValues
633 "interpolateHaloLayer_AoS_linear_with_constant_extrapolation(...)",
635 numberOfDoFsPerAxisInPatch,
640 if (marker.isInteriorFaceWithinPatch()) {
642 const int patchSize = numberOfDoFsPerAxisInPatch
643 * numberOfDoFsPerAxisInPatch;
645 const int patchSize = numberOfDoFsPerAxisInPatch * numberOfDoFsPerAxisInPatch
646 * numberOfDoFsPerAxisInPatch;
649 double* leftAdjacentPatchData = tarch::allocateMemory<double>(
650 patchSize * unknowns,
653 double* rightAdjacentPatchData = tarch::allocateMemory<double>(
654 patchSize * unknowns,
659 leftAdjacentPatchIndex = marker.getRelativePositionWithinFatherCell();
661 rightAdjacentPatchIndex = marker.getRelativePositionWithinFatherCell();
662 leftAdjacentPatchIndex(marker.getSelectedFaceNumber() % Dimensions)--;
664 internal::interpolateCell_AoS_linear(
665 leftAdjacentPatchIndex,
666 numberOfDoFsPerAxisInPatch,
668 coarseGridCellValues,
669 leftAdjacentPatchData,
672 internal::interpolateCell_AoS_linear(
673 rightAdjacentPatchIndex,
674 numberOfDoFsPerAxisInPatch,
676 coarseGridCellValues,
677 rightAdjacentPatchData,
681 internal::projectInterpolatedFineCellsOnHaloLayer_AoS(
683 numberOfDoFsPerAxisInPatch,
686 leftAdjacentPatchData,
687 rightAdjacentPatchData,
696 numberOfDoFsPerAxisInPatch,
699 coarseGridFaceValues,
704 logTraceOut(
"interpolateHaloLayer_AoS_linear_with_constant_extrapolation(...)"
711 int numberOfDoFsPerAxisInPatch,
713 const double* __restrict__ coarseGridCellValues,
714 double* __restrict__ fineGridCellValues
717 "interpolateCell_AoS_linear_with_constant_extrapolation(...)",
719 numberOfDoFsPerAxisInPatch,
723 internal::interpolateCell_AoS_linear(
724 marker.getRelativePositionWithinFatherCell(),
725 numberOfDoFsPerAxisInPatch,
727 coarseGridCellValues,
732 logTraceOut(
"interpolateCell_AoS_linear_with_constant_extrapolation(...)");
738 int numberOfDoFsPerAxisInPatch,
741 const double* __restrict__ coarseGridFaceValues,
742 double* __restrict__ fineGridFaceValues
745 "interpolateHaloLayer_AoS_linear_with_linear_extrapolation(...)",
747 numberOfDoFsPerAxisInPatch,
752 const int normal = marker.getSelectedFaceNumber() % Dimensions;
756 key(numberOfDoFsPerAxisInPatch, overlap, normal);
761 if (faceInterpolationMap.count(key) == 0) {
762 faceInterpolationMap[key] = internal::createLinearInterpolationMatrix(
763 numberOfDoFsPerAxisInPatch,
773 const int patchSize = numberOfDoFsPerAxisInPatch * 2 * overlap;
775 const int patchSize = numberOfDoFsPerAxisInPatch * numberOfDoFsPerAxisInPatch
780 "interpolateHaloLayer_AoS_linear_with_constant_extrapolation(...)",
781 marker.toString() <<
": " <<
P->toString()
784 P->batchedMultiplyAoS(
786 coarseGridFaceValues,
790 marker.getRelativePositionWithinFatherCell(),
795 logTraceOut(
"interpolateHaloLayer_AoS_linear_with_linear_extrapolation(...)");
801 int numberOfDoFsPerAxisInPatch,
804 const double* __restrict__ coarseGridCellValues,
805 const double* __restrict__ coarseGridFaceValues,
806 double* __restrict__ fineGridFaceValues
809 "interpolateHaloLayer_AoS_linear_with_linear_extrapolation(...)",
811 numberOfDoFsPerAxisInPatch,
816 if (marker.isInteriorFaceWithinPatch()) {
818 const int patchSize = numberOfDoFsPerAxisInPatch
819 * numberOfDoFsPerAxisInPatch;
821 const int patchSize = numberOfDoFsPerAxisInPatch * numberOfDoFsPerAxisInPatch
822 * numberOfDoFsPerAxisInPatch;
825 double* leftAdjacentPatchData = tarch::allocateMemory<double>(
826 patchSize * unknowns,
829 double* rightAdjacentPatchData = tarch::allocateMemory<double>(
830 patchSize * unknowns,
835 leftAdjacentPatchIndex = marker.getRelativePositionWithinFatherCell();
837 rightAdjacentPatchIndex = marker.getRelativePositionWithinFatherCell();
838 leftAdjacentPatchIndex(marker.getSelectedFaceNumber() % Dimensions)--;
840 internal::interpolateCell_AoS_linear(
841 leftAdjacentPatchIndex,
842 numberOfDoFsPerAxisInPatch,
844 coarseGridCellValues,
845 leftAdjacentPatchData,
848 internal::interpolateCell_AoS_linear(
849 rightAdjacentPatchIndex,
850 numberOfDoFsPerAxisInPatch,
852 coarseGridCellValues,
853 rightAdjacentPatchData,
857 internal::projectInterpolatedFineCellsOnHaloLayer_AoS(
859 numberOfDoFsPerAxisInPatch,
862 leftAdjacentPatchData,
863 rightAdjacentPatchData,
872 numberOfDoFsPerAxisInPatch,
875 coarseGridFaceValues,
880 logTraceOut(
"interpolateHaloLayer_AoS_linear_with_linear_extrapolation(...)");
886 int numberOfDoFsPerAxisInPatch,
888 const double* __restrict__ coarseGridCellValues,
889 double* __restrict__ fineGridCellValues
892 "interpolateCell_AoS_linear_with_linear_extrapolation(...)",
894 numberOfDoFsPerAxisInPatch,
898 internal::interpolateCell_AoS_linear(
899 marker.getRelativePositionWithinFatherCell(),
900 numberOfDoFsPerAxisInPatch,
902 coarseGridCellValues,
907 logTraceOut(
"interpolateCell_AoS_linear_with_linear_extrapolation(...)");
913 int numberOfDoFsPerAxisInPatch,
916 const double* __restrict__ coarseGridFaceValues,
917 double* __restrict__ fineGridFaceValues
920 "interpolateHaloLayer_AoS_linear_with_constant_extrapolation_and_linear_normal_interpolation(...)",
922 numberOfDoFsPerAxisInPatch,
927 const int normal = marker.getSelectedFaceNumber() % Dimensions;
931 key(numberOfDoFsPerAxisInPatch, overlap, normal);
936 if (faceInterpolationMap.count(key) == 0) {
937 faceInterpolationMap[key] = internal::createLinearInterpolationMatrix(
938 numberOfDoFsPerAxisInPatch,
948 const int patchSize = numberOfDoFsPerAxisInPatch * 2 * overlap;
950 const int patchSize = numberOfDoFsPerAxisInPatch * numberOfDoFsPerAxisInPatch
955 "interpolateHaloLayer_AoS_linear_with_constant_extrapolation_and_linear_normal_interpolation(...)",
956 marker.toString() <<
": " <<
P->toString()
959 P->batchedMultiplyAoS(
961 coarseGridFaceValues,
965 marker.getRelativePositionWithinFatherCell(),
971 "interpolateHaloLayer_AoS_linear_with_constant_extrapolation_and_linear_normal_interpolation(...)"
978 int numberOfDoFsPerAxisInPatch,
981 const double* __restrict__ coarseGridCellValues,
982 const double* __restrict__ coarseGridFaceValues,
983 double* __restrict__ fineGridFaceValues
986 "interpolateHaloLayer_AoS_linear_with_constant_extrapolation_and_linear_normal_interpolation(...)",
988 numberOfDoFsPerAxisInPatch,
993 if (marker.isInteriorFaceWithinPatch()) {
995 const int patchSize = numberOfDoFsPerAxisInPatch
996 * numberOfDoFsPerAxisInPatch;
998 const int patchSize = numberOfDoFsPerAxisInPatch * numberOfDoFsPerAxisInPatch
999 * numberOfDoFsPerAxisInPatch;
1002 double* leftAdjacentPatchData = tarch::allocateMemory<double>(
1003 patchSize * unknowns,
1006 double* rightAdjacentPatchData = tarch::allocateMemory<double>(
1007 patchSize * unknowns,
1012 leftAdjacentPatchIndex = marker.getRelativePositionWithinFatherCell();
1014 rightAdjacentPatchIndex = marker.getRelativePositionWithinFatherCell();
1015 leftAdjacentPatchIndex(marker.getSelectedFaceNumber() % Dimensions)--;
1017 internal::interpolateCell_AoS_linear(
1018 leftAdjacentPatchIndex,
1019 numberOfDoFsPerAxisInPatch,
1021 coarseGridCellValues,
1022 leftAdjacentPatchData,
1025 internal::interpolateCell_AoS_linear(
1026 rightAdjacentPatchIndex,
1027 numberOfDoFsPerAxisInPatch,
1029 coarseGridCellValues,
1030 rightAdjacentPatchData,
1034 internal::projectInterpolatedFineCellsOnHaloLayer_AoS(
1036 numberOfDoFsPerAxisInPatch,
1039 leftAdjacentPatchData,
1040 rightAdjacentPatchData,
1049 numberOfDoFsPerAxisInPatch,
1052 coarseGridFaceValues,
1058 "interpolateHaloLayer_AoS_linear_with_constant_extrapolation_and_linear_normal_interpolation(...)"
1065 int numberOfDoFsPerAxisInPatch,
1067 const double* __restrict__ coarseGridCellValues,
1068 double* __restrict__ fineGridCellValues
1071 "interpolateCell_AoS_linear_with_constant_extrapolation_and_linear_normal_interpolation(...)",
1073 numberOfDoFsPerAxisInPatch,
1077 internal::interpolateCell_AoS_linear(
1078 marker.getRelativePositionWithinFatherCell(),
1079 numberOfDoFsPerAxisInPatch,
1081 coarseGridCellValues,
1087 "interpolateCell_AoS_linear_with_constant_extrapolation_and_linear_normal_interpolation(...)"
1094 int numberOfDoFsPerAxisInPatch,
1097 const double* __restrict__ coarseGridFaceValues,
1098 double* __restrict__ fineGridFaceValues
1101 "interpolateHaloLayer_AoS_linear_with_linear_extrapolation_and_linear_normal_interpolation(...)",
1103 numberOfDoFsPerAxisInPatch,
1108 const int normal = marker.getSelectedFaceNumber() % Dimensions;
1112 key(numberOfDoFsPerAxisInPatch, overlap, normal);
1117 if (faceInterpolationMap.count(key) == 0) {
1118 faceInterpolationMap[key] = internal::createLinearInterpolationMatrix(
1119 numberOfDoFsPerAxisInPatch,
1129 const int patchSize = numberOfDoFsPerAxisInPatch * 2 * overlap;
1131 const int patchSize = numberOfDoFsPerAxisInPatch * numberOfDoFsPerAxisInPatch
1136 "interpolateHaloLayer_AoS_linear_with_constant_extrapolation_and_linear_normal_interpolation(...)",
1137 marker.toString() <<
": " <<
P->toString()
1140 P->batchedMultiplyAoS(
1142 coarseGridFaceValues,
1146 marker.getRelativePositionWithinFatherCell(),
1152 "interpolateHaloLayer_AoS_linear_with_linear_extrapolation_and_linear_normal_interpolation(...)"
1159 int numberOfDoFsPerAxisInPatch,
1162 const double* __restrict__ coarseGridCellValues,
1163 const double* __restrict__ coarseGridFaceValues,
1164 double* __restrict__ fineGridFaceValues
1167 "interpolateHaloLayer_AoS_linear_with_linear_extrapolation_and_linear_normal_interpolation(...)",
1169 numberOfDoFsPerAxisInPatch,
1174 if (marker.isInteriorFaceWithinPatch()) {
1176 const int patchSize = numberOfDoFsPerAxisInPatch
1177 * numberOfDoFsPerAxisInPatch;
1179 const int patchSize = numberOfDoFsPerAxisInPatch * numberOfDoFsPerAxisInPatch
1180 * numberOfDoFsPerAxisInPatch;
1183 double* leftAdjacentPatchData = tarch::allocateMemory<double>(
1184 patchSize * unknowns,
1187 double* rightAdjacentPatchData = tarch::allocateMemory<double>(
1188 patchSize * unknowns,
1193 leftAdjacentPatchIndex = marker.getRelativePositionWithinFatherCell();
1195 rightAdjacentPatchIndex = marker.getRelativePositionWithinFatherCell();
1196 leftAdjacentPatchIndex(marker.getSelectedFaceNumber() % Dimensions)--;
1198 internal::interpolateCell_AoS_linear(
1199 leftAdjacentPatchIndex,
1200 numberOfDoFsPerAxisInPatch,
1202 coarseGridCellValues,
1203 leftAdjacentPatchData,
1206 internal::interpolateCell_AoS_linear(
1207 rightAdjacentPatchIndex,
1208 numberOfDoFsPerAxisInPatch,
1210 coarseGridCellValues,
1211 rightAdjacentPatchData,
1215 internal::projectInterpolatedFineCellsOnHaloLayer_AoS(
1217 numberOfDoFsPerAxisInPatch,
1220 leftAdjacentPatchData,
1221 rightAdjacentPatchData,
1230 numberOfDoFsPerAxisInPatch,
1233 coarseGridFaceValues,
1239 "interpolateHaloLayer_AoS_linear_with_linear_extrapolation_and_linear_normal_interpolation(...)"
1246 int numberOfDoFsPerAxisInSourcePatch,
1247 int numberOfDoFsPerAxisInDestinationPatch,
1248 int haloSourcePatch,
1249 int haloDestinationPatch,
1251 const double* __restrict__ sourceValues,
1252 double* __restrict__ destinationValues,
1256 numberOfDoFsPerAxisInDestinationPatch > numberOfDoFsPerAxisInSourcePatch
1259 simtDforWithSchedulerInstructions(
1261 numberOfDoFsPerAxisInDestinationPatch + 2 * haloDestinationPatch,
1263 )
const int baseIndexDestination
1266 numberOfDoFsPerAxisInDestinationPatch + 2 * haloDestinationPatch
1269 for (
int unknown = 0; unknown < unknowns; unknown++) {
1270 destinationValues[baseIndexDestination + unknown] = 0.0;
1273 dfor(sourceVolume, numberOfDoFsPerAxisInSourcePatch + 2 * haloSourcePatch) {
1278 tarch::la::convertScalar<double>(sourceVolume
1284 1.0 / numberOfDoFsPerAxisInDestinationPatch
1286 tarch::la::convertScalar<double>(destinationVolume
1290 int outsidePatchAlongCoordinateAxis = 0;
1291 for (
int d = 0; d < Dimensions; d++) {
1296 if (sourceVolume(d) < haloSourcePatch)
1297 outsidePatchAlongCoordinateAxis++;
1298 if (sourceVolume(d) >= haloSourcePatch + numberOfDoFsPerAxisInSourcePatch)
1299 outsidePatchAlongCoordinateAxis++;
1301 const bool isFromUnitialisedDiagonal = outsidePatchAlongCoordinateAxis > 1;
1303 double weight = 1.0;
1304 for (
int d = 0; d < Dimensions; d++) {
1305 double distanceAlongCoordinateAxis =
std::abs(
1306 normalisedDestinationPosition(d) - normalisedSourcePosition(d)
1308 double h = 1.0 /
static_cast<double>(numberOfDoFsPerAxisInSourcePatch);
1309 distanceAlongCoordinateAxis = 1.0 - distanceAlongCoordinateAxis / h;
1310 weight *= std::max(0.0, distanceAlongCoordinateAxis);
1312#if !defined(SharedSYCL)
1318 if (isFromUnitialisedDiagonal) {
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
1331 amendedSourceVolume,
1332 numberOfDoFsPerAxisInSourcePatch + 2 * haloSourcePatch
1335 for (
int unknown = 0; unknown < unknowns; unknown++) {
1336 destinationValues[baseIndexDestination + unknown]
1337 += weight * sourceValues[baseIndexSource + unknown];
1346 int numberOfDoFsPerAxisInSourcePatch,
1347 int numberOfDoFsPerAxisInDestinationPatch,
1348 int haloSourcePatch,
1349 int haloDestinationPatch,
1351 const double* __restrict__ interpolationData,
1352 const int* __restrict__ columnIndices,
1353 const int* __restrict__ rowIndices,
1354 const double* __restrict__ sourceValues,
1355 double* __restrict__ destinationValues,
1359 = (numberOfDoFsPerAxisInDestinationPatch + 2 * haloDestinationPatch)
1360 * (numberOfDoFsPerAxisInDestinationPatch + 2 * haloDestinationPatch)
1361 * (numberOfDoFsPerAxisInDestinationPatch + 2 * haloDestinationPatch);
1363 for (
int row = 0; row < destinationSize; row++) {
1364 for (
int unknown = 0; unknown < unknowns; unknown++) {
1365 destinationValues[row * unknowns + unknown] = 0.0;
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];
1379 int numberOfDoFsPerAxisInSourcePatch,
1380 int numberOfDoFsPerAxisInDestinationPatch,
1381 int haloSourcePatch,
1382 int haloDestinationPatch,
1384 const double* __restrict__ sourceValues,
1385 double* __restrict__ destinationValues,
1391 double* centerValues =
new double[unknowns];
1392 int numberOfDerivatives = (Dimensions == 2) ? 5 : 9;
1394 dfor(volume, numberOfDoFsPerAxisInSourcePatch + 2) {
1396 sourceVolume = volume
1398 int outsidePatchAlongCoordinateAxis = 0;
1399 for (
int d = 0; d < Dimensions; d++) {
1400 if (sourceVolume(d) < haloSourcePatch || sourceVolume(d) >= numberOfDoFsPerAxisInSourcePatch + haloSourcePatch) {
1401 outsidePatchAlongCoordinateAxis++;
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)
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;
1418 if (sourceVolume(d1) == haloSourcePatch and sourceVolume(d2) == haloSourcePatch + numberOfDoFsPerAxisInSourcePatch - 1) {
1419 stencilOffset(d1) = 1;
1420 stencilOffset(d2) = -1;
1422 if (sourceVolume(d1) == haloSourcePatch + numberOfDoFsPerAxisInSourcePatch - 1 and sourceVolume(d2) == haloSourcePatch) {
1423 stencilOffset(d1) = -1;
1424 stencilOffset(d2) = 1;
1426 if (sourceVolume(d1) == haloSourcePatch + numberOfDoFsPerAxisInSourcePatch - 1 and sourceVolume(d2) == haloSourcePatch + numberOfDoFsPerAxisInSourcePatch - 1) {
1427 stencilOffset(d1) = -1;
1428 stencilOffset(d2) = -1;
1433 if (outsidePatchAlongCoordinateAxis <= 1) {
1434 int stencilSize = 0;
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)
1447 stencil[stencilSize] = stencilOffset + offset
1450 sourceVolume + offset + stencilOffset
1452 numberOfDoFsPerAxisInSourcePatch + 2 * haloSourcePatch
1460 numberOfDoFsPerAxisInSourcePatch + 2 * haloSourcePatch
1462 for (
int unknown = 0; unknown < unknowns; unknown++) {
1463 centerValues[unknown] = sourceValues[centerIndex * unknowns + unknown];
1467 for (
int k = 0; k < stencilSize; k++) {
1468 for (
int unknown = 0; unknown < unknowns; unknown++) {
1472 ) = sourceValues[indices[k] * unknowns + unknown]
1473 - centerValues[unknown];
1477 for (
int row = 0; row < stencilSize; row++) {
1480 for (
int a = 0;
a < Dimensions;
a++) {
1481 A(row, column) = cell(
a);
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);
1501 inverseR.
rows() * c.cols(),
1509 ) * numberOfDoFsPerAxisInDestinationPatch
1510 / numberOfDoFsPerAxisInSourcePatch
1515 ) * numberOfDoFsPerAxisInDestinationPatch
1516 / numberOfDoFsPerAxisInSourcePatch
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)
1527 destinationSpan = destinationCorner2 - destinationCorner1;
1528 int destinationDimensionsProduct = 1;
1529 for (
int d = 0; d < Dimensions; d++)
1530 destinationDimensionsProduct *= destinationSpan(d);
1534 = (1.0 / numberOfDoFsPerAxisInDestinationPatch)
1537 = (1.0 / numberOfDoFsPerAxisInSourcePatch)
1540 cornerDisplacement = (normalisedCornerPosition
1541 - normalisedSourcePosition)
1542 * (
double)numberOfDoFsPerAxisInSourcePatch;
1544 double* matrixData =
new double
1545 [destinationDimensionsProduct * numberOfDerivatives];
1548 dfor(localDestinationVolume, destinationSpan) {
1550 displacement = cornerDisplacement
1551 + ((
double)numberOfDoFsPerAxisInSourcePatch
1552 / numberOfDoFsPerAxisInDestinationPatch)
1554 double>(localDestinationVolume);
1556 for (
int a = 0;
a < Dimensions;
a++) {
1557 matrixData[writeIndex] = displacement(
a);
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);
1569 destinationDimensionsProduct,
1570 numberOfDerivatives,
1575 int destinationVolumeNumber = 0;
1576 dfor(localDestinationVolume, destinationSpan) {
1578 destinationVolume = destinationCorner1 + localDestinationVolume;
1581 numberOfDoFsPerAxisInDestinationPatch + 2 * haloDestinationPatch
1583 for (
int unknown = 0; unknown < unknowns; unknown++) {
1584 destinationValues[destinationIndex * unknowns + unknown]
1585 = centerValues[unknown]
1586 + destinationValuesDelta(destinationVolumeNumber, unknown);
1588 destinationVolumeNumber++;
1592 delete[] centerValues;
1597 int numberOfDoFsPerAxisInSourcePatch,
1598 int numberOfDoFsPerAxisInDestinationPatch,
1599 int haloSourcePatch,
1600 int haloDestinationPatch,
1602 const double* __restrict__ sourceValues,
1603 double* __restrict__ destinationValues,
1607 numberOfDoFsPerAxisInDestinationPatch > numberOfDoFsPerAxisInSourcePatch
1610 simtDforWithSchedulerInstructions(
1612 numberOfDoFsPerAxisInDestinationPatch + 2 * haloDestinationPatch,
1614 )
const int baseIndexDestination
1617 numberOfDoFsPerAxisInDestinationPatch + 2 * haloDestinationPatch
1620 for (
int unknown = 0; unknown < unknowns; unknown++) {
1621 destinationValues[baseIndexDestination + unknown] = 0.0;
1630 int numberOfDoFsPerAxisInPatch,
1632 const double* __restrict__ coarseGridCellValues,
1633 double* __restrict__ fineGridCellValues
1636 "interpolateCell_AoS_linear_with_linear_extrapolation_and_linear_normal_interpolation(...)",
1638 numberOfDoFsPerAxisInPatch,
1642 internal::interpolateCell_AoS_linear(
1643 marker.getRelativePositionWithinFatherCell(),
1644 numberOfDoFsPerAxisInPatch,
1646 coarseGridCellValues,
1652 "interpolateCell_AoS_linear_with_linear_extrapolation_and_linear_normal_interpolation(...)"
1658 int numberOfDoFsPerAxisInPatch,
1660 const double* __restrict__ coarseGridCellValues,
1661 double* __restrict__ fineGridCellValues,
1662 bool extrapolateLinearly
1664 static internal::CellInterpolationMap cellInterpolationMap;
1665 internal::CellInterpolationOperatorKey key(numberOfDoFsPerAxisInPatch);
1668 if (cellInterpolationMap.count(key) == 0) {
1669 cellInterpolationMap[key] = internal::createLinearInterpolationMatrix(
1670 numberOfDoFsPerAxisInPatch,
1678 const int patchSize = numberOfDoFsPerAxisInPatch * numberOfDoFsPerAxisInPatch;
1680 const int patchSize = numberOfDoFsPerAxisInPatch * numberOfDoFsPerAxisInPatch
1681 * numberOfDoFsPerAxisInPatch;
1686 coarseGridCellValues,
1689 serialiseMarkerIn3x3PatchAssembly(
1690 relativePositionWithinFatherCell,
1691 numberOfDoFsPerAxisInPatch
1698 int numberOfDoFsPerAxisInPatch,
1701 const double* __restrict__ coarseGridFaceValues,
1702 double* __restrict__ fineGridFaceValues
1704 const int normal =
marker.getSelectedFaceNumber() % Dimensions;
1705 const int i = (
normal + 1) % Dimensions;
1707 const int j = (
normal + 2) % Dimensions;
1712 ) =
marker.getRelativePositionWithinFatherCell()(
i)
1713 * numberOfDoFsPerAxisInPatch / 3;
1716 ) =
marker.getRelativePositionWithinFatherCell()(
j)
1717 * numberOfDoFsPerAxisInPatch / 3;
1722 ) = ((
marker.getRelativePositionWithinFatherCell()(
i) + 1
1723 ) * numberOfDoFsPerAxisInPatch
1728 ) = ((
marker.getRelativePositionWithinFatherCell()(
j) + 1
1729 ) * numberOfDoFsPerAxisInPatch
1736 ) = (
marker.getRelativePositionWithinFatherCell()(
i)
1737 * numberOfDoFsPerAxisInPatch)
1741 ) = (
marker.getRelativePositionWithinFatherCell()(
j)
1742 * numberOfDoFsPerAxisInPatch)
1747 faceDimensions(normal) = (
overlap - 1) / 3 + 1;
1751 double* centerValues =
new double[
unknowns];
1753 dfor(kCoarse, faceDimensions) {
1756 ) = (
marker.getSelectedFaceNumber() < Dimensions)
1757 ? overlap - kCoarse(normal) - 1
1762 stencilOffset(i) = 1;
1763 if (
center(i) == numberOfDoFsPerAxisInPatch - 1)
1764 stencilOffset(i) = -1;
1767 stencilOffset(j) = 1;
1768 if (
center(j) == numberOfDoFsPerAxisInPatch - 1)
1769 stencilOffset(j) = -1;
1772 int stencilSize = 0;
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)
1788 center + offset + stencilOffset
1790 numberOfDoFsPerAxisInPatch,
1800 numberOfDoFsPerAxisInPatch,
1804 for (
int unknown = 0; unknown <
unknowns; unknown++) {
1805 centerValues[unknown] = coarseGridFaceValues[
index *
unknowns + unknown];
1809 for (
int k = 0;
k < stencilSize;
k++) {
1810 for (
int unknown = 0; unknown <
unknowns; unknown++) {
1814 ) = coarseGridFaceValues[indices[
k] *
unknowns + unknown]
1815 - centerValues[unknown];
1819 int numberOfDerivatives = (Dimensions == 2) ? 5 : 9;
1821 for (
int row = 0;
row < stencilSize;
row++) {
1824 for (
int a = 0;
a < Dimensions;
a++) {
1825 A(row, column) = cell(
a);
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);
1850 if (kCoarse(i) == 0)
1852 ) = (
marker.getRelativePositionWithinFatherCell()(
i)
1853 * numberOfDoFsPerAxisInPatch)
1856 if (kCoarse(j) == 0)
1858 ) = (
marker.getRelativePositionWithinFatherCell()(
j)
1859 * numberOfDoFsPerAxisInPatch)
1863 int modulus = (
marker.getRelativePositionWithinFatherCell()(
i) + 1)
1864 * numberOfDoFsPerAxisInPatch;
1865 if (modulus % 3 != 0)
1866 modulus = modulus % 3;
1869 if (kCoarse(i) == faceDimensions(i) - 1)
1870 fineEnd(i) = modulus;
1873 modulus = (
marker.getRelativePositionWithinFatherCell()(
j) + 1)
1874 * numberOfDoFsPerAxisInPatch;
1875 if (modulus % 3 != 0)
1876 modulus = modulus % 3;
1879 if (kCoarse(j) == faceDimensions(j) - 1)
1880 fineEnd(j) = modulus;
1884 fineDimensions(normal) = std::min(overlap - kCoarse(normal) * 3, 3);
1885 if (fineDimensions(normal) % 3 != 0)
1886 fineDimensions(normal) = fineDimensions(normal) % 3;
1888 fineDimensions(normal) = 3;
1890 int fineDimensionsProduct = fineDimensions(normal) * fineDimensions(i)
1891 * fineDimensions(j);
1892#elif Dimensions == 2
1893 int fineDimensionsProduct = fineDimensions(normal) * fineDimensions(i);
1896 fineMatrix(fineDimensionsProduct, numberOfDerivatives);
1899 dfor(kFine, fineDimensions) {
1903 for (
int a = 0;
a < Dimensions;
a++) {
1904 fineMatrix(row, column) = (1.0 / 3.0) * cell(
a);
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);
1917 int fineCellNumber = 0;
1919 fineCorner = kCoarse * 3 + fineStart - fineDataStart;
1921 ) = (
marker.getSelectedFaceNumber() < Dimensions)
1922 ? (kCoarse(normal) * 3)
1925 dfor(kFine, fineDimensions) {
1929 numberOfDoFsPerAxisInPatch,
1934 for (
int unknown = 0; unknown <
unknowns; unknown++) {
1936 = centerValues[unknown] + fineValues(fineCellNumber, unknown);
1942 delete[] centerValues;
1948 int numberOfDoFsPerAxisInPatch,
1951 const double* __restrict__ coarseGridFaceValues,
1952 double* __restrict__ fineGridFaceValues
1954 const int normal =
marker.getSelectedFaceNumber() % Dimensions;
1955 const int i = (
normal + 1) % Dimensions;
1957 const int j = (
normal + 2) % Dimensions;
1962 ) =
marker.getRelativePositionWithinFatherCell()(
i)
1963 * numberOfDoFsPerAxisInPatch / 3;
1966 ) =
marker.getRelativePositionWithinFatherCell()(
j)
1967 * numberOfDoFsPerAxisInPatch / 3;
1972 ) = ((
marker.getRelativePositionWithinFatherCell()(
i) + 1
1973 ) * numberOfDoFsPerAxisInPatch
1978 ) = ((
marker.getRelativePositionWithinFatherCell()(
j) + 1
1979 ) * numberOfDoFsPerAxisInPatch
1986 ) = (
marker.getRelativePositionWithinFatherCell()(
i)
1987 * numberOfDoFsPerAxisInPatch)
1991 ) = (
marker.getRelativePositionWithinFatherCell()(
j)
1992 * numberOfDoFsPerAxisInPatch)
1997 faceDimensions(normal) = (
overlap - 1) / 3 + 1;
1999 double* centerValues =
new double[
unknowns];
2004 dfor(kCoarse, faceDimensions) {
2007 ) = (
marker.getSelectedFaceNumber() < Dimensions)
2008 ? overlap - kCoarse(normal) - 1
2013 stencilOffset(i) = 1;
2014 if (
center(i) == numberOfDoFsPerAxisInPatch - 1)
2015 stencilOffset(i) = -1;
2018 stencilOffset(j) = 1;
2019 if (
center(j) == numberOfDoFsPerAxisInPatch - 1)
2020 stencilOffset(j) = -1;
2023 int stencilSize = 0;
2029 center + offset + stencilOffset
2031 numberOfDoFsPerAxisInPatch,
2039 for (
int axis = 0;
axis < Dimensions;
axis++) {
2045 center + stencilOffset + offset,
2046 numberOfDoFsPerAxisInPatch,
2052 if (
center(axis) < numberOfDoFsPerAxisInPatch - 2) {
2057 center + stencilOffset + offset,
2058 numberOfDoFsPerAxisInPatch,
2068 numberOfDoFsPerAxisInPatch,
2072 for (
int unknown = 0; unknown <
unknowns; unknown++) {
2073 centerValues[unknown] = coarseGridFaceValues[
index *
unknowns + unknown];
2076 for (
int k = 0;
k < stencilSize;
k++) {
2077 for (
int unknown = 0; unknown <
unknowns; unknown++) {
2081 ) = coarseGridFaceValues[indices[
k] *
unknowns + unknown]
2082 - centerValues[unknown];
2086 int numberOfDerivatives = (Dimensions == 2) ? 9 : 19;
2089 for (
int row = 0;
row < stencilSize;
row++) {
2092 for (
int a = 0;
a < Dimensions;
a++) {
2093 A(row, column) = cell(
a);
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);
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);
2124 if (kCoarse(i) == 0)
2126 ) = (
marker.getRelativePositionWithinFatherCell()(
i)
2127 * numberOfDoFsPerAxisInPatch)
2130 if (kCoarse(j) == 0)
2132 ) = (
marker.getRelativePositionWithinFatherCell()(
j)
2133 * numberOfDoFsPerAxisInPatch)
2137 int modulus = (
marker.getRelativePositionWithinFatherCell()(
i) + 1)
2138 * numberOfDoFsPerAxisInPatch;
2139 if (modulus % 3 != 0)
2140 modulus = modulus % 3;
2143 if (kCoarse(i) == faceDimensions(i) - 1)
2144 fineEnd(i) = modulus;
2147 modulus = (
marker.getRelativePositionWithinFatherCell()(
j) + 1)
2148 * numberOfDoFsPerAxisInPatch;
2149 if (modulus % 3 != 0)
2150 modulus = modulus % 3;
2153 if (kCoarse(j) == faceDimensions(j) - 1)
2154 fineEnd(j) = modulus;
2158 fineDimensions(normal) = std::min(overlap - kCoarse(normal) * 3, 3);
2159 if (fineDimensions(normal) % 3 != 0)
2160 fineDimensions(normal) = fineDimensions(normal) % 3;
2162 fineDimensions(normal) = 3;
2165 int fineDimensionsProduct = fineDimensions(normal) * fineDimensions(i)
2166 * fineDimensions(j);
2167#elif Dimensions == 2
2168 int fineDimensionsProduct = fineDimensions(normal) * fineDimensions(i);
2172 fineMatrix(fineDimensionsProduct, numberOfDerivatives);
2175 dfor(kFine, fineDimensions) {
2179 for (
int a = 0;
a < Dimensions;
a++) {
2180 fineMatrix(row, column) = (1.0 / 3.0) * cell(
a);
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);
2189 for (
int a = 0;
a < Dimensions;
a++) {
2190 for (
int b =
a;
b < Dimensions;
b++) {
2191 for (
int c = b;
c < Dimensions;
c++) {
2195 ) = (1.0 / 162.0) * cell(
a) * cell(b) * cell(c);
2204 int fineCellNumber = 0;
2206 fineCorner = kCoarse * 3 + fineStart - fineDataStart;
2208 ) = (
marker.getSelectedFaceNumber() < Dimensions)
2209 ? (kCoarse(normal) * 3)
2212 dfor(kFine, fineDimensions) {
2216 numberOfDoFsPerAxisInPatch,
2221 for (
int unknown = 0; unknown <
unknowns; unknown++) {
2223 = centerValues[unknown] + fineValues(fineCellNumber, unknown);
2229 delete[] centerValues;
2235 int numberOfDoFsPerAxisInPatch,
2238 const double* __restrict__ interpolationData,
2239 const int* __restrict__ columnIndices,
2240 const int* __restrict__ rowIndices,
2242 int rowIndicesStart,
2243 const double* __restrict__ coarseGridFaceValues,
2244 double* __restrict__ fineGridFaceValues
2247 numberOfDoFsPerAxisInPatch
2250 numberOfDoFsPerAxisInPatch
2253 const int normal = marker.getSelectedFaceNumber() % Dimensions;
2254 const int matrixColumns = 2 * overlap
2256 pow(numberOfDoFsPerAxisInPatch, Dimensions - 1);
2257 const int matrixRows = overlap
2258 * std::pow(numberOfDoFsPerAxisInPatch, Dimensions - 1);
2260 fineFaceDimensions(normal) = overlap;
2261 coarseFaceDimensions(normal) = 2 * overlap;
2263 double* rearrangedCoarseData =
new double[matrixColumns * unknowns];
2265 dfor(kCoarse, coarseFaceDimensions) {
2268 numberOfDoFsPerAxisInPatch,
2272 int newIndex = (marker.getSelectedFaceNumber() >= Dimensions)
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;
2280 for (
int unknown = 0; unknown < unknowns; unknown++) {
2281 rearrangedCoarseData[newIndex * unknowns + unknown] = coarseGridFaceValues
2282 [currentIndex * unknowns + unknown];
2286 dfor(kFine, fineFaceDimensions) {
2287 auto fineCoords = kFine;
2288 fineCoords(normal) = (marker.getSelectedFaceNumber() < Dimensions)
2289 ? overlap - kFine(normal) - 1
2290 : overlap + kFine(normal);
2294 numberOfDoFsPerAxisInPatch,
2299 int row = kFine(normal);
2301 for (
int d = 1; d < Dimensions; d++) {
2302 row += kFine((normal + d) % Dimensions) * base;
2303 base *= numberOfDoFsPerAxisInPatch;
2306 for (
int unknown = 0; unknown < unknowns; unknown++) {
2307 fineGridFaceValues[fineIndex * unknowns + unknown] = 0.0;
2310 int columnStart = rowIndices[rowIndicesStart + row];
2311 int columnEnd = rowIndices[rowIndicesStart + row + 1];
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];
2323 delete[] rearrangedCoarseData;
2328 int numberOfDoFsPerAxisInPatch,
2331 const double* __restrict__ normalInterpolationMatrix1d,
2332 const double* __restrict__ tangentialInterpolationMatrix1d,
2333 const double* __restrict__ coarseGridFaceValues,
2334 double* __restrict__ fineGridFaceValues
2337 "interpolateHaloLayer_AoS_tensor_product(...)",
2339 numberOfDoFsPerAxisInPatch,
2343 const int normal = marker.getSelectedFaceNumber() % Dimensions;
2345 dfore(kFine, numberOfDoFsPerAxisInPatch, normal, 0) {
2346 for (
int iFine = 0; iFine < overlap; iFine++) {
2349 ) = (marker.getSelectedFaceNumber() < Dimensions)
2351 : 2 * overlap - iFine - 1;
2355 numberOfDoFsPerAxisInPatch,
2361 "interpolateHaloLayer_AoS_tensor_product(...)",
2362 "clear dof " << dest
2364 for (
int unknown = 0; unknown < unknowns; unknown++) {
2365 fineGridFaceValues[destLinearised * unknowns + unknown] = 0.0;
2368 dfore(kCoarse, numberOfDoFsPerAxisInPatch, normal, 0) {
2369 for (
int iCoarse = 0; iCoarse < 2 * overlap; iCoarse++) {
2372 ) = (marker.getSelectedFaceNumber() < Dimensions)
2374 : 2 * overlap - iCoarse - 1;
2377 numberOfDoFsPerAxisInPatch,
2382 double weight = 1.0;
2383 for (
int d = 0; d < Dimensions; d++) {
2390 int index = col + row * 2 * overlap;
2392 "interpolateHaloLayer_AoS_tensor_product(...)",
2393 "(" << row <<
"," << col <<
")=" << index
2394 <<
" (normal contribution)"
2396 weight *= normalInterpolationMatrix1d[index];
2401 + marker.getRelativePositionWithinFatherCell()(d
2402 ) * numberOfDoFsPerAxisInPatch;
2405 int index = col + row * numberOfDoFsPerAxisInPatch;
2407 "interpolateHaloLayer_AoS_tensor_product(...)",
2408 "(" << row <<
"," << col <<
")=" << index
2409 <<
" (tangential contribution)"
2411 weight *= tangentialInterpolationMatrix1d[index];
2417 "interpolateHaloLayer_AoS_tensor_product(...)",
2418 "add dof " << src <<
" to " << dest <<
" with weight " << weight
2421 for (
int unknown = 0; unknown < unknowns; unknown++) {
2422 fineGridFaceValues[destLinearised * unknowns + unknown]
2424 * coarseGridFaceValues[srcLinearised * unknowns + unknown];
2431 logTraceOut(
"interpolateHaloLayer_AoS_tensor_product(...)");
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
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
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
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
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
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
2513 [[maybe_unused]]
int numberOfDoFsPerAxisInPatch,
2514 [[maybe_unused]]
int unknowns,
2515 [[maybe_unused]]
const double* __restrict__ coarseGridCellValues,
2516 [[maybe_unused]]
double* __restrict__ fineGridCellValues
2523 [[maybe_unused]]
int numberOfDoFsPerAxisInPatch,
2524 [[maybe_unused]]
int unknowns,
2525 [[maybe_unused]]
const double* __restrict__ coarseGridCellValues,
2526 [[maybe_unused]]
double* __restrict__ fineGridCellValues
#define assertion4(expr, param0, param1, param2, param3)
#define assertionEquals1(lhs, rhs, a)
#define assertionEquals(lhs, rhs)
#define logDebug(methodName, logMacroMessageStream)
#define logTraceOut(methodName)
#define logTraceInWith4Arguments(methodName, argument0, argument1, argument2, argument3)
#define logTraceInWith3Arguments(methodName, argument0, argument1, argument2)
#define logTraceInWith1Argument(methodName, argument0)
#define dfor3(counter)
Shortcut For dfor(counter,3)
#define dfore(counter, max, dim, value)
This is an exclusive d-dimensional for loop.
#define dfor(counter, max)
d-dimensional Loop
#define enddforx
I prefer to use this macro for dforx instead of a closing bracket as many syntax parser fail otherwis...
tarch::logging::Log _log("::")
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
Create a lock around a boolean semaphore region.
void free()
Free the lock.
const float const float const float struct part *restrict struct part *restrict const float a
LoopPlacement
Guide loop-level parallelism.
CPUGPUMethod int dLinearised(const tarch::la::Vector< Dimensions, int > &counter, int max)
Map d-dimensional vector onto integer index.
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)
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.
Provide information about selected face.