Peano 4
Loading...
Searching...
No Matches
Interpolation.cpp
Go to the documentation of this file.
1#include "Interpolation.h"
2#include "Enumeration.h"
3
4#include "peano4/utils/Loop.h"
6
10
11namespace {
12 [[maybe_unused]] tarch::logging::Log _log( "toolbox::blockstructured" );
13
14 tarch::multicore::BooleanSemaphore interpolationMapSemaphore;
15}
16
18 #if Dimensions==2
19 const int patchSize = numberOfDoFsPerAxisInPatch * numberOfDoFsPerAxisInPatch;
20 const int numberOfPatches = 3*3;
21 #else
22 const int patchSize = numberOfDoFsPerAxisInPatch * numberOfDoFsPerAxisInPatch * numberOfDoFsPerAxisInPatch;
23 const int numberOfPatches = 3*3*3;
24 #endif
25
26 tarch::la::DynamicMatrix* result = new tarch::la::DynamicMatrix(patchSize*numberOfPatches,patchSize);
27
28 int currentRow = 0;
29 dfor3(patchIndex)
30 dfor(volumeIndex,numberOfDoFsPerAxisInPatch) {
31 tarch::la::Vector<Dimensions,int> sourceCell = (volumeIndex + patchIndex*numberOfDoFsPerAxisInPatch) / 3;
32 int sourceCellLinearised = peano4::utils::dLinearised(sourceCell,numberOfDoFsPerAxisInPatch);
33 (*result)(currentRow,sourceCellLinearised) = 1.0;
34 currentRow++;
35 }
37
38 logDebug( "createPiecewiseConstantInterpolationMatrix(int)", "created new volumetric matrix for " << numberOfDoFsPerAxisInPatch << ": " << result->toString() );
39 return result;
40}
41
43 [[maybe_unused]] int numberOfDoFsPerAxisInPatch,
44 [[maybe_unused]] int normal,
45 [[maybe_unused]] int overlap
46) {
47 tarch::la::DynamicMatrix P1d(3,1,{ {1.0},{1.0},{1.0} });
48 P1d.replicateRows( 3, numberOfDoFsPerAxisInPatch, 1, true );
49
51 P1d, numberOfDoFsPerAxisInPatch, normal
52 );
53}
54
56 const tarch::la::DynamicMatrix& P1d,
57 int numberOfDoFsPerAxisInPatch,
58 int normal
59) {
60 assertionEquals1( P1d.cols(), numberOfDoFsPerAxisInPatch, P1d.toString() );
61 assertionEquals1( P1d.rows(), 3*numberOfDoFsPerAxisInPatch, P1d.toString() );
62
63 #if Dimensions==3
64 tarch::la::DynamicMatrix* P = new tarch::la::DynamicMatrix( P1d, P1d, false);
65 #else
67 #endif
68
69 int pattern = 0;
70 switch (normal) {
71 case 0:
72 pattern = 1;
73 break;
74 case 1:
75 pattern = numberOfDoFsPerAxisInPatch;
76 break;
77 case 2:
78 pattern = numberOfDoFsPerAxisInPatch*numberOfDoFsPerAxisInPatch;
79 break;
80 }
81
82 P->insertEmptyColumns(pattern,pattern,pattern);
83 P->replicateRows( pattern, 2, pattern, false );
84 logDebug( "createInterpolationMatrixFrom1dTemplateByInsertingZeroColsAndRows(...)", "matrix for normal=" << normal << ": " << P->toString() );
85 return P;
86}
87
89 const tarch::la::DynamicMatrix& P1d,
90 int numberOfDoFsPerAxisInPatch,
91 int normal
92) {
93 assertionEquals1( P1d.cols(), numberOfDoFsPerAxisInPatch, P1d.toString() );
94 assertionEquals1( P1d.rows(), 3*numberOfDoFsPerAxisInPatch, P1d.toString() );
95
96 #if Dimensions==3
97 tarch::la::DynamicMatrix* P = new tarch::la::DynamicMatrix( P1d, P1d, false);
98 #else
100 #endif
101
102 int pattern = 0;
103 switch (normal) {
104 case 0:
105 pattern = 1;
106 break;
107 case 1:
108 pattern = numberOfDoFsPerAxisInPatch;
109 break;
110 case 2:
111 pattern = numberOfDoFsPerAxisInPatch*numberOfDoFsPerAxisInPatch;
112 break;
113 }
114
115 P->scale(0.5);
116
117 // Split matrix into blocks of pattern columns and copy each block once.
118 // Matrix grows by factor of two which reflects the fact that the overlap
119 // is of size 1.
120 P->replicateCols( pattern, 2, 0, false );
121 P->replicateRows( pattern, 2, 0, false );
122 logDebug( "createInterpolationMatrixFrom1dTemplateByLinearInterpolationAlongNormal(...)", "matrix for normal=" << normal << ": " << P->toString() );
123 return P;
124}
125
127 int numberOfDoFsPerAxisInPatch,
128 int normal,
129 bool extrapolateLinearly,
130 bool interpolateLinearlyBetweenRealAndRestrictedVolumes
131) {
133 {1.0/3.0, 2.0/3.0, 0.0},
134 { 0.0, 3.0/3.0, 0.0},
135 { 0.0, 2.0/3.0, 1.0/3.0}
136 });
137 P1d.replicateRows( 3, numberOfDoFsPerAxisInPatch, 1, true );
138 P1d.removeColumn(0);
139 P1d.removeColumn(numberOfDoFsPerAxisInPatch);
140
141 if (extrapolateLinearly) {
142 P1d(0,0) = 2.0/3.0 + 1.0/3.0 * 2.0;
143 P1d(0,1) = -1.0/3.0;
144 P1d(numberOfDoFsPerAxisInPatch*3-1,numberOfDoFsPerAxisInPatch-1) = 2.0/3.0 + 1.0/3.0 * 2.0;
145 P1d(numberOfDoFsPerAxisInPatch*3-1,numberOfDoFsPerAxisInPatch-2) = -1.0/3.0;
146 }
147 else {
148 P1d(0,0) = 1.0;
149 P1d(0,1) = 0.0;
150 P1d(numberOfDoFsPerAxisInPatch*3-1,numberOfDoFsPerAxisInPatch-1) = 1.0;
151 P1d(numberOfDoFsPerAxisInPatch*3-1,numberOfDoFsPerAxisInPatch-2) = 0.0;
152 }
153
154 logDebug( "createLinearInterpolationMatrix(...)", "1d matrix: " << P1d.toString() );
155
156 if (interpolateLinearlyBetweenRealAndRestrictedVolumes) {
158 P1d, numberOfDoFsPerAxisInPatch, normal
159 );
160 }
161 else {
163 P1d, numberOfDoFsPerAxisInPatch, normal
164 );
165 }
166}
167
169 logTraceInWith1Argument( "createLinearInterpolationMatrix(...)", numberOfDoFsPerAxisInPatch );
171 {1.0/3.0, 2.0/3.0, 0.0},
172 { 0.0, 3.0/3.0, 0.0},
173 { 0.0, 2.0/3.0, 1.0/3.0}
174 });
175 P1d.replicateRows( 3, numberOfDoFsPerAxisInPatch, 1, true );
176 P1d.removeColumn(0);
177 P1d.removeColumn(numberOfDoFsPerAxisInPatch);
178
179 if (extrapolateLinearly) {
180 P1d(0,0) = 2.0/3.0 + 1.0/3.0 * 2.0;
181 P1d(0,1) = -1.0/3.0;
182 P1d(numberOfDoFsPerAxisInPatch*3-1,numberOfDoFsPerAxisInPatch-1) = 2.0/3.0 + 1.0/3.0 * 2.0;
183 P1d(numberOfDoFsPerAxisInPatch*3-1,numberOfDoFsPerAxisInPatch-2) = -1.0/3.0;
184 }
185 else {
186 P1d(0,0) = 1.0;
187 P1d(0,1) = 0.0;
188 P1d(numberOfDoFsPerAxisInPatch*3-1,numberOfDoFsPerAxisInPatch-1) = 1.0;
189 P1d(numberOfDoFsPerAxisInPatch*3-1,numberOfDoFsPerAxisInPatch-2) = 0.0;
190 }
191
192 logDebug( "createLinearInterpolationMatrix(...)", "1d interpolation matrix: " << P1d.toString(true) );
193
194 tarch::la::DynamicMatrix P2d( P1d, P1d, false);
195 logDebug( "createLinearInterpolationMatrix(...)", "2d interpolation matrix: " << P2d.toString(true) );
196
197 logTraceOut( "createLinearInterpolationMatrix(...)" );
198 #if Dimensions==3
199 return new tarch::la::DynamicMatrix( P2d, P1d, false);
200 #else
201 return new tarch::la::DynamicMatrix( P2d );
202 #endif
203}
204
207 int numberOfDoFsPerAxisInPatch,
208 int overlap,
209 int unknowns,
210 const double* __restrict__ fineGridCellValuesLeft,
211 const double* __restrict__ fineGridCellValuesRight,
212 double* fineGridFaceValues
213) {
214 logTraceInWith3Arguments( "projectInterpolatedFineCellsOnHaloLayer_AoS(...)", marker.toString(), numberOfDoFsPerAxisInPatch, overlap );
215
216 const int normal = marker.getSelectedFaceNumber() % Dimensions;
217
218 dfore(kFine,numberOfDoFsPerAxisInPatch,normal,0) {
219 for (int iFine=0; iFine<overlap; iFine++) {
220 tarch::la::Vector<Dimensions,int> fineVolumeLeftDest = kFine;
221 tarch::la::Vector<Dimensions,int> fineVolumeLeftSrc = kFine;
222
223 fineVolumeLeftDest(normal) = iFine;
224 fineVolumeLeftSrc(normal) = numberOfDoFsPerAxisInPatch - overlap + iFine;
225
226 tarch::la::Vector<Dimensions,int> fineVolumeRightDest = kFine;
227 tarch::la::Vector<Dimensions,int> fineVolumeRightSrc = kFine;
228
229 fineVolumeRightDest(normal) = iFine + overlap;
230 fineVolumeRightSrc(normal) = iFine;
231
232 int fineVolumeLeftDestLinearised = serialiseVoxelIndexInOverlap(
233 fineVolumeLeftDest,
234 numberOfDoFsPerAxisInPatch, overlap, normal
235 );
236 int fineVolumeRightDestLinearised = serialiseVoxelIndexInOverlap(
237 fineVolumeRightDest,
238 numberOfDoFsPerAxisInPatch, overlap, normal
239 );
240 int fineVolumeLeftSrcLinearised = peano4::utils::dLinearised(
241 fineVolumeLeftSrc,
242 numberOfDoFsPerAxisInPatch
243 );
244 int fineVolumeRightSrcLinearised = peano4::utils::dLinearised(
245 fineVolumeRightSrc,
246 numberOfDoFsPerAxisInPatch
247 );
248
249 for (int j=0; j<unknowns; j++) {
250 fineGridFaceValues[fineVolumeLeftDestLinearised*unknowns+j] = fineGridCellValuesLeft[fineVolumeLeftSrcLinearised*unknowns+j];
251 fineGridFaceValues[fineVolumeRightDestLinearised*unknowns+j] = fineGridCellValuesRight[fineVolumeRightSrcLinearised*unknowns+j];
252 }
253 }
254 }
255
256 logTraceOut( "projectInterpolatedFineCellsOnHaloLayer_AoS(...)" );
257}
258
259//
260// ==========================================
261// Piece-wise constant interpolation routines
262// ==========================================
263//
266 int numberOfDoFsPerAxisInPatch,
267 int overlap,
268 int unknowns,
269 const double* __restrict__ coarseGridFaceValues,
270 double* __restrict__ fineGridFaceValues
271) {
272 logTraceInWith4Arguments( "interpolateHaloLayer_AoS_piecewise_constant(...)", marker.toString(), numberOfDoFsPerAxisInPatch, overlap, unknowns );
273
274 const int normal = marker.getSelectedFaceNumber() % Dimensions;
275
276 static internal::FaceInterpolationMap faceInterpolationMap;
277 internal::FaceInterpolationOperatorKey key(numberOfDoFsPerAxisInPatch,overlap,normal);
278
279 tarch::multicore::Lock lock(interpolationMapSemaphore);
280 if ( faceInterpolationMap.count(key)==0 ) {
281 faceInterpolationMap[key] = internal::createPiecewiseConstantInterpolationMatrix(numberOfDoFsPerAxisInPatch,normal,overlap);
282 }
283 tarch::la::DynamicMatrix* P = faceInterpolationMap[key];
284 lock.free();
285
286 #if Dimensions==2
287 const int patchSize = numberOfDoFsPerAxisInPatch*2*overlap;
288 #else
289 const int patchSize = numberOfDoFsPerAxisInPatch*numberOfDoFsPerAxisInPatch*2*overlap;
290 #endif
291
292 logDebug( "interpolateHaloLayer_AoS_piecewise_constant(...)", marker.toString() << ": " << P->toString() );
293
294 P->batchedMultiplyAoS(
295 fineGridFaceValues, // image
296 coarseGridFaceValues, // preimage
297 unknowns, // batch size, i.e. how often to apply it in one AoS rush
298 patchSize, // result size, i.e. size of image
300 marker.getRelativePositionWithinFatherCell(),
301 normal
302 ) * patchSize
303 );
304
305 logTraceOut( "interpolateHaloLayer_AoS_piecewise_constant(...)" );
306}
307
310 int numberOfDoFsPerAxisInPatch,
311 int overlap,
312 int unknowns,
313 const double* __restrict__ coarseGridCellValues,
314 const double* __restrict__ coarseGridFaceValues,
315 double* __restrict__ fineGridFaceValues
316) {
317 logTraceInWith4Arguments( "interpolateHaloLayer_AoS_piecewise_constant(...)", marker.toString(), numberOfDoFsPerAxisInPatch, overlap, unknowns );
318
319 if ( marker.isInteriorFaceWithinPatch() ) {
320 #if Dimensions==2
321 const int patchSize = numberOfDoFsPerAxisInPatch*numberOfDoFsPerAxisInPatch;
322 #else
323 const int patchSize = numberOfDoFsPerAxisInPatch*numberOfDoFsPerAxisInPatch*numberOfDoFsPerAxisInPatch;
324 #endif
325
326 double* leftAdjacentPatchData = tarch::allocateMemory<double>(patchSize * unknowns, tarch::MemoryLocation::Heap);
327 double* rightAdjacentPatchData = tarch::allocateMemory<double>(patchSize * unknowns, tarch::MemoryLocation::Heap);
328
329 tarch::la::Vector<Dimensions,int> leftAdjacentPatchIndex = marker.getRelativePositionWithinFatherCell();
330 tarch::la::Vector<Dimensions,int> rightAdjacentPatchIndex = marker.getRelativePositionWithinFatherCell();
331 leftAdjacentPatchIndex( marker.getSelectedFaceNumber()%Dimensions )--;
332
333 internal::interpolateCell_AoS_piecewise_constant(
334 leftAdjacentPatchIndex,
335 numberOfDoFsPerAxisInPatch,
336 unknowns,
337 coarseGridCellValues,
338 leftAdjacentPatchData
339 );
340 internal::interpolateCell_AoS_piecewise_constant(
341 rightAdjacentPatchIndex,
342 numberOfDoFsPerAxisInPatch,
343 unknowns,
344 coarseGridCellValues,
345 rightAdjacentPatchData
346 );
347
348 internal::projectInterpolatedFineCellsOnHaloLayer_AoS(
349 marker, numberOfDoFsPerAxisInPatch, overlap, unknowns,
350 leftAdjacentPatchData, rightAdjacentPatchData,
351 fineGridFaceValues
352 );
353
354 tarch::freeMemory(leftAdjacentPatchData, tarch::MemoryLocation::Heap );
355 tarch::freeMemory(rightAdjacentPatchData, tarch::MemoryLocation::Heap );
356 }
357 else {
359 marker,
360 numberOfDoFsPerAxisInPatch,
361 overlap,
362 unknowns,
363 coarseGridFaceValues,
364 fineGridFaceValues
365 );
366 }
367
368 logTraceOut( "interpolateHaloLayer_AoS_piecewise_constant(...)" );
369}
370
373 int numberOfDoFsPerAxisInPatch,
374 int unknowns,
375 const double* __restrict__ coarseGridCellValues,
376 double* __restrict__ fineGridCellValues
377) {
378 logTraceInWith3Arguments("interpolateCell_AoS_piecewise_constant(...)", marker.toString(), numberOfDoFsPerAxisInPatch, unknowns );
379
380 internal::interpolateCell_AoS_piecewise_constant(
381 marker.getRelativePositionWithinFatherCell(),
382 numberOfDoFsPerAxisInPatch,
383 unknowns,
384 coarseGridCellValues,
385 fineGridCellValues
386 );
387
388 logTraceOut("interpolateCell_AoS_piecewise_constant(...)");
389}
390
392 const tarch::la::Vector<Dimensions,int>& relativePositionWithinFatherCell,
393 int numberOfDoFsPerAxisInPatch,
394 int unknowns,
395 const double* __restrict__ coarseGridCellValues,
396 double* __restrict__ fineGridCellValues
397) {
398 static internal::CellInterpolationMap cellInterpolationMap;
399 internal::CellInterpolationOperatorKey key(numberOfDoFsPerAxisInPatch);
400
401 tarch::multicore::Lock lock(interpolationMapSemaphore);
402 if ( cellInterpolationMap.count(key)==0 ) {
403 cellInterpolationMap[key] = internal::createPiecewiseConstantInterpolationMatrix(numberOfDoFsPerAxisInPatch);
404 }
405 tarch::la::DynamicMatrix* P = cellInterpolationMap[key];
406 lock.free();
407
408 #if Dimensions==2
409 const int patchSize = numberOfDoFsPerAxisInPatch*numberOfDoFsPerAxisInPatch;
410 #else
411 const int patchSize = numberOfDoFsPerAxisInPatch*numberOfDoFsPerAxisInPatch*numberOfDoFsPerAxisInPatch;
412 #endif
413
414 P->batchedMultiplyAoS(
415 fineGridCellValues, // image
416 coarseGridCellValues, // preimage
417 unknowns, // batch size, i.e. how often to apply it in one AoS rush
418 patchSize, // result size, i.e. size of image
419 serialiseMarkerIn3x3PatchAssembly(relativePositionWithinFatherCell, numberOfDoFsPerAxisInPatch)
420 );
421}
422
423//
424// =============================
425// Linear interpolation routines
426// =============================
427//
430 int numberOfDoFsPerAxisInPatch,
431 int overlap,
432 int unknowns,
433 const double* __restrict__ coarseGridFaceValues,
434 double* __restrict__ fineGridFaceValues
435) {
436 logTraceInWith4Arguments( "interpolateHaloLayer_AoS_linear_with_constant_extrapolation(...)", marker.toString(), numberOfDoFsPerAxisInPatch, overlap, unknowns );
437
438 const int normal = marker.getSelectedFaceNumber() % Dimensions;
439
440 static internal::FaceInterpolationMap faceInterpolationMap;
441 internal::FaceInterpolationOperatorKey key(numberOfDoFsPerAxisInPatch,overlap,normal);
442
443 assertionEquals(overlap,1);
444
445 tarch::multicore::Lock lock(interpolationMapSemaphore);
446 if ( faceInterpolationMap.count(key)==0 ) {
447 faceInterpolationMap[key] = internal::createLinearInterpolationMatrix(numberOfDoFsPerAxisInPatch,normal,false,false);
448 }
449 tarch::la::DynamicMatrix* P = faceInterpolationMap[key];
450 lock.free();
451
452 #if Dimensions==2
453 const int patchSize = numberOfDoFsPerAxisInPatch*2*overlap;
454 #else
455 const int patchSize = numberOfDoFsPerAxisInPatch*numberOfDoFsPerAxisInPatch*2*overlap;
456 #endif
457
458 logDebug( "interpolateHaloLayer_AoS_linear_with_constant_extrapolation(...)", marker.toString() << ": " << P->toString() );
459
460 P->batchedMultiplyAoS(
461 fineGridFaceValues, // image
462 coarseGridFaceValues, // preimage
463 unknowns, // batch size, i.e. how often to apply it in one AoS rush
464 patchSize, // result size, i.e. size of image
466 marker.getRelativePositionWithinFatherCell(),
467 normal
468 ) * patchSize
469 );
470
471 logTraceOut( "interpolateHaloLayer_AoS_linear_with_constant_extrapolation(...)" );
472}
473
476 int numberOfDoFsPerAxisInPatch,
477 int overlap,
478 int unknowns,
479 const double* __restrict__ coarseGridCellValues,
480 const double* __restrict__ coarseGridFaceValues,
481 double* __restrict__ fineGridFaceValues
482) {
483 logTraceInWith4Arguments( "interpolateHaloLayer_AoS_linear_with_constant_extrapolation(...)", marker.toString(), numberOfDoFsPerAxisInPatch, overlap, unknowns );
484
485 if ( marker.isInteriorFaceWithinPatch() ) {
486 #if Dimensions==2
487 const int patchSize = numberOfDoFsPerAxisInPatch*numberOfDoFsPerAxisInPatch;
488 #else
489 const int patchSize = numberOfDoFsPerAxisInPatch*numberOfDoFsPerAxisInPatch*numberOfDoFsPerAxisInPatch;
490 #endif
491
492 double* leftAdjacentPatchData = tarch::allocateMemory<double>(patchSize * unknowns, tarch::MemoryLocation::Heap);
493 double* rightAdjacentPatchData = tarch::allocateMemory<double>(patchSize * unknowns, tarch::MemoryLocation::Heap);
494
495 tarch::la::Vector<Dimensions,int> leftAdjacentPatchIndex = marker.getRelativePositionWithinFatherCell();
496 tarch::la::Vector<Dimensions,int> rightAdjacentPatchIndex = marker.getRelativePositionWithinFatherCell();
497 leftAdjacentPatchIndex( marker.getSelectedFaceNumber()%Dimensions )--;
498
499 internal::interpolateCell_AoS_linear(
500 leftAdjacentPatchIndex,
501 numberOfDoFsPerAxisInPatch,
502 unknowns,
503 coarseGridCellValues,
504 leftAdjacentPatchData,
505 false
506 );
507 internal::interpolateCell_AoS_linear(
508 rightAdjacentPatchIndex,
509 numberOfDoFsPerAxisInPatch,
510 unknowns,
511 coarseGridCellValues,
512 rightAdjacentPatchData,
513 false
514 );
515
516 internal::projectInterpolatedFineCellsOnHaloLayer_AoS(
517 marker, numberOfDoFsPerAxisInPatch, overlap, unknowns,
518 leftAdjacentPatchData, rightAdjacentPatchData,
519 fineGridFaceValues
520 );
521
522 tarch::freeMemory(leftAdjacentPatchData, tarch::MemoryLocation::Heap );
523 tarch::freeMemory(rightAdjacentPatchData, tarch::MemoryLocation::Heap );
524 }
525 else {
527 marker,
528 numberOfDoFsPerAxisInPatch,
529 overlap,
530 unknowns,
531 coarseGridFaceValues,
532 fineGridFaceValues
533 );
534 }
535
536 logTraceOut( "interpolateHaloLayer_AoS_linear_with_constant_extrapolation(...)" );
537}
538
541 int numberOfDoFsPerAxisInPatch,
542 int unknowns,
543 const double* __restrict__ coarseGridCellValues,
544 double* __restrict__ fineGridCellValues
545) {
546 logTraceInWith3Arguments("interpolateCell_AoS_linear_with_constant_extrapolation(...)", marker.toString(), numberOfDoFsPerAxisInPatch, unknowns );
547
548 internal::interpolateCell_AoS_linear(
549 marker.getRelativePositionWithinFatherCell(),
550 numberOfDoFsPerAxisInPatch,
551 unknowns,
552 coarseGridCellValues,
553 fineGridCellValues,
554 false
555 );
556
557 logTraceOut("interpolateCell_AoS_linear_with_constant_extrapolation(...)");
558}
559
562 int numberOfDoFsPerAxisInPatch,
563 int overlap,
564 int unknowns,
565 const double* __restrict__ coarseGridFaceValues,
566 double* __restrict__ fineGridFaceValues
567) {
568 logTraceInWith4Arguments( "interpolateHaloLayer_AoS_linear_with_linear_extrapolation(...)", marker.toString(), numberOfDoFsPerAxisInPatch, overlap, unknowns );
569
570 const int normal = marker.getSelectedFaceNumber() % Dimensions;
571
572 static internal::FaceInterpolationMap faceInterpolationMap;
573 internal::FaceInterpolationOperatorKey key(numberOfDoFsPerAxisInPatch,overlap,normal);
574
575 assertionEquals(overlap,1);
576
577 tarch::multicore::Lock lock(interpolationMapSemaphore);
578 if ( faceInterpolationMap.count(key)==0 ) {
579 faceInterpolationMap[key] = internal::createLinearInterpolationMatrix(numberOfDoFsPerAxisInPatch,normal,true,false);
580 }
581 tarch::la::DynamicMatrix* P = faceInterpolationMap[key];
582 lock.free();
583
584 #if Dimensions==2
585 const int patchSize = numberOfDoFsPerAxisInPatch*2*overlap;
586 #else
587 const int patchSize = numberOfDoFsPerAxisInPatch*numberOfDoFsPerAxisInPatch*2*overlap;
588 #endif
589
590 logDebug( "interpolateHaloLayer_AoS_linear_with_constant_extrapolation(...)", marker.toString() << ": " << P->toString() );
591
592 P->batchedMultiplyAoS(
593 fineGridFaceValues, // image
594 coarseGridFaceValues, // preimage
595 unknowns, // batch size, i.e. how often to apply it in one AoS rush
596 patchSize, // result size, i.e. size of image
598 marker.getRelativePositionWithinFatherCell(),
599 normal
600 ) * patchSize
601 );
602
603 logTraceOut( "interpolateHaloLayer_AoS_linear_with_linear_extrapolation(...)" );
604}
605
608 int numberOfDoFsPerAxisInPatch,
609 int overlap,
610 int unknowns,
611 const double* __restrict__ coarseGridCellValues,
612 const double* __restrict__ coarseGridFaceValues,
613 double* __restrict__ fineGridFaceValues
614) {
615 logTraceInWith4Arguments( "interpolateHaloLayer_AoS_linear_with_linear_extrapolation(...)", marker.toString(), numberOfDoFsPerAxisInPatch, overlap, unknowns );
616
617 if ( marker.isInteriorFaceWithinPatch() ) {
618 #if Dimensions==2
619 const int patchSize = numberOfDoFsPerAxisInPatch*numberOfDoFsPerAxisInPatch;
620 #else
621 const int patchSize = numberOfDoFsPerAxisInPatch*numberOfDoFsPerAxisInPatch*numberOfDoFsPerAxisInPatch;
622 #endif
623
624 double* leftAdjacentPatchData = tarch::allocateMemory<double>(patchSize * unknowns, tarch::MemoryLocation::Heap);
625 double* rightAdjacentPatchData = tarch::allocateMemory<double>(patchSize * unknowns, tarch::MemoryLocation::Heap);
626
627 tarch::la::Vector<Dimensions,int> leftAdjacentPatchIndex = marker.getRelativePositionWithinFatherCell();
628 tarch::la::Vector<Dimensions,int> rightAdjacentPatchIndex = marker.getRelativePositionWithinFatherCell();
629 leftAdjacentPatchIndex( marker.getSelectedFaceNumber()%Dimensions )--;
630
631 internal::interpolateCell_AoS_linear(
632 leftAdjacentPatchIndex,
633 numberOfDoFsPerAxisInPatch,
634 unknowns,
635 coarseGridCellValues,
636 leftAdjacentPatchData,
637 true
638 );
639 internal::interpolateCell_AoS_linear(
640 rightAdjacentPatchIndex,
641 numberOfDoFsPerAxisInPatch,
642 unknowns,
643 coarseGridCellValues,
644 rightAdjacentPatchData,
645 true
646 );
647
648 internal::projectInterpolatedFineCellsOnHaloLayer_AoS(
649 marker, numberOfDoFsPerAxisInPatch, overlap, unknowns,
650 leftAdjacentPatchData, rightAdjacentPatchData,
651 fineGridFaceValues
652 );
653
654 tarch::freeMemory(leftAdjacentPatchData, tarch::MemoryLocation::Heap );
655 tarch::freeMemory(rightAdjacentPatchData, tarch::MemoryLocation::Heap );
656 }
657 else {
659 marker,
660 numberOfDoFsPerAxisInPatch,
661 overlap,
662 unknowns,
663 coarseGridFaceValues,
664 fineGridFaceValues
665 );
666 }
667
668 logTraceOut( "interpolateHaloLayer_AoS_linear_with_linear_extrapolation(...)" );
669}
670
673 int numberOfDoFsPerAxisInPatch,
674 int unknowns,
675 const double* __restrict__ coarseGridCellValues,
676 double* __restrict__ fineGridCellValues
677) {
678 logTraceInWith3Arguments("interpolateCell_AoS_linear_with_linear_extrapolation(...)", marker.toString(), numberOfDoFsPerAxisInPatch, unknowns );
679
680 internal::interpolateCell_AoS_linear(
681 marker.getRelativePositionWithinFatherCell(),
682 numberOfDoFsPerAxisInPatch,
683 unknowns,
684 coarseGridCellValues,
685 fineGridCellValues,
686 true
687 );
688
689 logTraceOut("interpolateCell_AoS_linear_with_linear_extrapolation(...)");
690}
691
694 int numberOfDoFsPerAxisInPatch,
695 int overlap,
696 int unknowns,
697 const double* __restrict__ coarseGridFaceValues,
698 double* __restrict__ fineGridFaceValues
699) {
700 logTraceInWith4Arguments( "interpolateHaloLayer_AoS_linear_with_constant_extrapolation_and_linear_normal_interpolation(...)", marker.toString(), numberOfDoFsPerAxisInPatch, overlap, unknowns );
701
702 const int normal = marker.getSelectedFaceNumber() % Dimensions;
703
704 static internal::FaceInterpolationMap faceInterpolationMap;
705 internal::FaceInterpolationOperatorKey key(numberOfDoFsPerAxisInPatch,overlap,normal);
706
707 assertionEquals(overlap,1);
708
709 tarch::multicore::Lock lock(interpolationMapSemaphore);
710 if ( faceInterpolationMap.count(key)==0 ) {
711 faceInterpolationMap[key] = internal::createLinearInterpolationMatrix(numberOfDoFsPerAxisInPatch,normal,false,true);
712 }
713 tarch::la::DynamicMatrix* P = faceInterpolationMap[key];
714 lock.free();
715
716 #if Dimensions==2
717 const int patchSize = numberOfDoFsPerAxisInPatch*2*overlap;
718 #else
719 const int patchSize = numberOfDoFsPerAxisInPatch*numberOfDoFsPerAxisInPatch*2*overlap;
720 #endif
721
722 logDebug( "interpolateHaloLayer_AoS_linear_with_constant_extrapolation_and_linear_normal_interpolation(...)", marker.toString() << ": " << P->toString() );
723
724 P->batchedMultiplyAoS(
725 fineGridFaceValues, // image
726 coarseGridFaceValues, // preimage
727 unknowns, // batch size, i.e. how often to apply it in one AoS rush
728 patchSize, // result size, i.e. size of image
730 marker.getRelativePositionWithinFatherCell(),
731 normal
732 ) * patchSize
733 );
734
735 logTraceOut( "interpolateHaloLayer_AoS_linear_with_constant_extrapolation_and_linear_normal_interpolation(...)" );
736}
737
740 int numberOfDoFsPerAxisInPatch,
741 int overlap,
742 int unknowns,
743 const double* __restrict__ coarseGridCellValues,
744 const double* __restrict__ coarseGridFaceValues,
745 double* __restrict__ fineGridFaceValues
746) {
747 logTraceInWith4Arguments( "interpolateHaloLayer_AoS_linear_with_constant_extrapolation_and_linear_normal_interpolation(...)", marker.toString(), numberOfDoFsPerAxisInPatch, overlap, unknowns );
748
749 if ( marker.isInteriorFaceWithinPatch() ) {
750 #if Dimensions==2
751 const int patchSize = numberOfDoFsPerAxisInPatch*numberOfDoFsPerAxisInPatch;
752 #else
753 const int patchSize = numberOfDoFsPerAxisInPatch*numberOfDoFsPerAxisInPatch*numberOfDoFsPerAxisInPatch;
754 #endif
755
756 double* leftAdjacentPatchData = tarch::allocateMemory<double>(patchSize * unknowns, tarch::MemoryLocation::Heap);
757 double* rightAdjacentPatchData = tarch::allocateMemory<double>(patchSize * unknowns, tarch::MemoryLocation::Heap);
758
759 tarch::la::Vector<Dimensions,int> leftAdjacentPatchIndex = marker.getRelativePositionWithinFatherCell();
760 tarch::la::Vector<Dimensions,int> rightAdjacentPatchIndex = marker.getRelativePositionWithinFatherCell();
761 leftAdjacentPatchIndex( marker.getSelectedFaceNumber()%Dimensions )--;
762
763 internal::interpolateCell_AoS_linear(
764 leftAdjacentPatchIndex,
765 numberOfDoFsPerAxisInPatch,
766 unknowns,
767 coarseGridCellValues,
768 leftAdjacentPatchData,
769 false
770 );
771 internal::interpolateCell_AoS_linear(
772 rightAdjacentPatchIndex,
773 numberOfDoFsPerAxisInPatch,
774 unknowns,
775 coarseGridCellValues,
776 rightAdjacentPatchData,
777 false
778 );
779
780 internal::projectInterpolatedFineCellsOnHaloLayer_AoS(
781 marker, numberOfDoFsPerAxisInPatch, overlap, unknowns,
782 leftAdjacentPatchData, rightAdjacentPatchData,
783 fineGridFaceValues
784 );
785
786 tarch::freeMemory(leftAdjacentPatchData, tarch::MemoryLocation::Heap );
787 tarch::freeMemory(rightAdjacentPatchData, tarch::MemoryLocation::Heap );
788 }
789 else {
791 marker,
792 numberOfDoFsPerAxisInPatch,
793 overlap,
794 unknowns,
795 coarseGridFaceValues,
796 fineGridFaceValues
797 );
798 }
799
800 logTraceOut( "interpolateHaloLayer_AoS_linear_with_constant_extrapolation_and_linear_normal_interpolation(...)" );
801}
802
805 int numberOfDoFsPerAxisInPatch,
806 int unknowns,
807 const double* __restrict__ coarseGridCellValues,
808 double* __restrict__ fineGridCellValues
809) {
810 logTraceInWith3Arguments("interpolateCell_AoS_linear_with_constant_extrapolation_and_linear_normal_interpolation(...)", marker.toString(), numberOfDoFsPerAxisInPatch, unknowns );
811
812 internal::interpolateCell_AoS_linear(
813 marker.getRelativePositionWithinFatherCell(),
814 numberOfDoFsPerAxisInPatch,
815 unknowns,
816 coarseGridCellValues,
817 fineGridCellValues,
818 false
819 );
820
821 logTraceOut("interpolateCell_AoS_linear_with_constant_extrapolation_and_linear_normal_interpolation(...)");
822}
823
826 int numberOfDoFsPerAxisInPatch,
827 int overlap,
828 int unknowns,
829 const double* __restrict__ coarseGridFaceValues,
830 double* __restrict__ fineGridFaceValues
831) {
832 logTraceInWith4Arguments( "interpolateHaloLayer_AoS_linear_with_linear_extrapolation_and_linear_normal_interpolation(...)", marker.toString(), numberOfDoFsPerAxisInPatch, overlap, unknowns );
833
834 const int normal = marker.getSelectedFaceNumber() % Dimensions;
835
836 static internal::FaceInterpolationMap faceInterpolationMap;
837 internal::FaceInterpolationOperatorKey key(numberOfDoFsPerAxisInPatch,overlap,normal);
838
839 assertionEquals(overlap,1);
840
841 tarch::multicore::Lock lock(interpolationMapSemaphore);
842 if ( faceInterpolationMap.count(key)==0 ) {
843 faceInterpolationMap[key] = internal::createLinearInterpolationMatrix(numberOfDoFsPerAxisInPatch,normal,true,true);
844 }
845 tarch::la::DynamicMatrix* P = faceInterpolationMap[key];
846 lock.free();
847
848 #if Dimensions==2
849 const int patchSize = numberOfDoFsPerAxisInPatch*2*overlap;
850 #else
851 const int patchSize = numberOfDoFsPerAxisInPatch*numberOfDoFsPerAxisInPatch*2*overlap;
852 #endif
853
854 logDebug( "interpolateHaloLayer_AoS_linear_with_constant_extrapolation_and_linear_normal_interpolation(...)", marker.toString() << ": " << P->toString() );
855
856 P->batchedMultiplyAoS(
857 fineGridFaceValues, // image
858 coarseGridFaceValues, // preimage
859 unknowns, // batch size, i.e. how often to apply it in one AoS rush
860 patchSize, // result size, i.e. size of image
862 marker.getRelativePositionWithinFatherCell(),
863 normal
864 ) * patchSize
865 );
866
867 logTraceOut( "interpolateHaloLayer_AoS_linear_with_linear_extrapolation_and_linear_normal_interpolation(...)" );
868}
869
872 int numberOfDoFsPerAxisInPatch,
873 int overlap,
874 int unknowns,
875 const double* __restrict__ coarseGridCellValues,
876 const double* __restrict__ coarseGridFaceValues,
877 double* __restrict__ fineGridFaceValues
878) {
879 logTraceInWith4Arguments( "interpolateHaloLayer_AoS_linear_with_linear_extrapolation_and_linear_normal_interpolation(...)", marker.toString(), numberOfDoFsPerAxisInPatch, overlap, unknowns );
880
881 if ( marker.isInteriorFaceWithinPatch() ) {
882 #if Dimensions==2
883 const int patchSize = numberOfDoFsPerAxisInPatch*numberOfDoFsPerAxisInPatch;
884 #else
885 const int patchSize = numberOfDoFsPerAxisInPatch*numberOfDoFsPerAxisInPatch*numberOfDoFsPerAxisInPatch;
886 #endif
887
888 double* leftAdjacentPatchData = tarch::allocateMemory<double>(patchSize * unknowns, tarch::MemoryLocation::Heap);
889 double* rightAdjacentPatchData = tarch::allocateMemory<double>(patchSize * unknowns, tarch::MemoryLocation::Heap);
890
891 tarch::la::Vector<Dimensions,int> leftAdjacentPatchIndex = marker.getRelativePositionWithinFatherCell();
892 tarch::la::Vector<Dimensions,int> rightAdjacentPatchIndex = marker.getRelativePositionWithinFatherCell();
893 leftAdjacentPatchIndex( marker.getSelectedFaceNumber()%Dimensions )--;
894
895 internal::interpolateCell_AoS_linear(
896 leftAdjacentPatchIndex,
897 numberOfDoFsPerAxisInPatch,
898 unknowns,
899 coarseGridCellValues,
900 leftAdjacentPatchData,
901 true
902 );
903 internal::interpolateCell_AoS_linear(
904 rightAdjacentPatchIndex,
905 numberOfDoFsPerAxisInPatch,
906 unknowns,
907 coarseGridCellValues,
908 rightAdjacentPatchData,
909 true
910 );
911
912 internal::projectInterpolatedFineCellsOnHaloLayer_AoS(
913 marker, numberOfDoFsPerAxisInPatch, overlap, unknowns,
914 leftAdjacentPatchData, rightAdjacentPatchData,
915 fineGridFaceValues
916 );
917
918 tarch::freeMemory(leftAdjacentPatchData, tarch::MemoryLocation::Heap );
919 tarch::freeMemory(rightAdjacentPatchData, tarch::MemoryLocation::Heap );
920 }
921 else {
923 marker,
924 numberOfDoFsPerAxisInPatch,
925 overlap,
926 unknowns,
927 coarseGridFaceValues,
928 fineGridFaceValues
929 );
930 }
931
932 logTraceOut( "interpolateHaloLayer_AoS_linear_with_linear_extrapolation_and_linear_normal_interpolation(...)" );
933}
934
935
937 int numberOfDoFsPerAxisInSourcePatch,
938 int numberOfDoFsPerAxisInDestinationPatch,
939 int haloSourcePatch,
940 int haloDestinationPatch,
941 int unknowns,
942 const double* __restrict__ sourceValues,
943 double* __restrict__ destinationValues,
944 ::peano4::utils::LoopPlacement parallelisation
945) {
946 assertion( numberOfDoFsPerAxisInDestinationPatch>numberOfDoFsPerAxisInSourcePatch );
947
948 simtDforWithSchedulerInstructions( destinationVolume, numberOfDoFsPerAxisInDestinationPatch+2*haloDestinationPatch, parallelisation )
949 const int baseIndexDestination = peano4::utils::dLinearised(destinationVolume,numberOfDoFsPerAxisInDestinationPatch+2*haloDestinationPatch) * unknowns;
950 for (int unknown=0; unknown<unknowns; unknown++) {
951 destinationValues[baseIndexDestination+unknown] = 0.0;
952 }
953
954 dfor( sourceVolume, numberOfDoFsPerAxisInSourcePatch+2*haloSourcePatch ) {
955 tarch::la::Vector<Dimensions, double> normalisedSourcePosition =
957 tarch::la::Vector<Dimensions, double>(1.0/numberOfDoFsPerAxisInSourcePatch),
958 tarch::la::convertScalar<double>(sourceVolume) + tarch::la::Vector<Dimensions, double>(0.5 - haloSourcePatch)
959 );
960 tarch::la::Vector<Dimensions, double> normalisedDestinationPosition =
962 tarch::la::Vector<Dimensions, double>(1.0/numberOfDoFsPerAxisInDestinationPatch),
963 tarch::la::convertScalar<double>(destinationVolume) + tarch::la::Vector<Dimensions, double>(0.5 - haloDestinationPatch)
964 );
965
966 int outsidePatchAlongCoordinateAxis = 0;
967 for (int d=0; d<Dimensions; d++) {
968// outsidePatchAlongCoordinateAxis += ( sourceVolume(d) < haloSourcePatch ) ? 1 : 0;
969// outsidePatchAlongCoordinateAxis += ( sourceVolume(d) >= haloSourcePatch+numberOfDoFsPerAxisInSourcePatch ) ? 1 : 0;
970 if ( sourceVolume(d) < haloSourcePatch ) outsidePatchAlongCoordinateAxis++;
971 if ( sourceVolume(d) >= haloSourcePatch+numberOfDoFsPerAxisInSourcePatch ) outsidePatchAlongCoordinateAxis++;
972 }
973 const bool isFromUnitialisedDiagonal = outsidePatchAlongCoordinateAxis>1;
974
975 double weight = 1.0;
976 for (int d=0; d<Dimensions; d++) {
977 double distanceAlongCoordinateAxis = std::abs( normalisedDestinationPosition(d) - normalisedSourcePosition(d) );
978 double h = 1.0 / static_cast<double>(numberOfDoFsPerAxisInSourcePatch);
979 distanceAlongCoordinateAxis = 1.0 - distanceAlongCoordinateAxis / h ;
980 weight *= std::max(0.0, distanceAlongCoordinateAxis);
981 }
982 #if !defined(SharedSYCL)
983 assertion( tarch::la::greaterEquals(weight, 0.0) );
984 assertion( tarch::la::smallerEquals(weight, 1.0) );
985 #endif
986
987 tarch::la::Vector<Dimensions, int> amendedSourceVolume = sourceVolume;
988 if (isFromUnitialisedDiagonal) {
989 #pragma unroll
990 for (int d=0; d<Dimensions; d++) {
991 amendedSourceVolume(d) = std::max( amendedSourceVolume(d), haloSourcePatch );
992 amendedSourceVolume(d) = std::min( amendedSourceVolume(d), haloSourcePatch + numberOfDoFsPerAxisInSourcePatch - 1 );
993 }
994 }
995
996 int baseIndexSource = peano4::utils::dLinearised(amendedSourceVolume, numberOfDoFsPerAxisInSourcePatch+2*haloSourcePatch) * unknowns;
997 for (int unknown=0; unknown<unknowns; unknown++) {
998 destinationValues[baseIndexDestination+unknown] += weight * sourceValues[baseIndexSource+unknown];
999 }
1000 }
1002}
1003
1005 int numberOfDoFsPerAxisInSourcePatch,
1006 int numberOfDoFsPerAxisInDestinationPatch,
1007 int haloSourcePatch,
1008 int haloDestinationPatch,
1009 int unknowns,
1010 const double* __restrict__ sourceValues,
1011 double* __restrict__ destinationValues,
1012 ::peano4::utils::LoopPlacement parallelisation
1013) {
1014 assertion( numberOfDoFsPerAxisInDestinationPatch>numberOfDoFsPerAxisInSourcePatch );
1015
1016 simtDforWithSchedulerInstructions( destinationVolume, numberOfDoFsPerAxisInDestinationPatch+2*haloDestinationPatch, parallelisation )
1017 const int baseIndexDestination = peano4::utils::dLinearised(destinationVolume,numberOfDoFsPerAxisInDestinationPatch+2*haloDestinationPatch) * unknowns;
1018 for (int unknown=0; unknown<unknowns; unknown++) {
1019 destinationValues[baseIndexDestination+unknown] = 0.0;
1020 }
1022}
1023
1024
1025
1028 int numberOfDoFsPerAxisInPatch,
1029 int unknowns,
1030 const double* __restrict__ coarseGridCellValues,
1031 double* __restrict__ fineGridCellValues
1032) {
1033 logTraceInWith3Arguments("interpolateCell_AoS_linear_with_linear_extrapolation_and_linear_normal_interpolation(...)", marker.toString(), numberOfDoFsPerAxisInPatch, unknowns );
1034
1035 internal::interpolateCell_AoS_linear(
1036 marker.getRelativePositionWithinFatherCell(),
1037 numberOfDoFsPerAxisInPatch,
1038 unknowns,
1039 coarseGridCellValues,
1040 fineGridCellValues,
1041 true
1042 );
1043
1044 logTraceOut("interpolateCell_AoS_linear_with_linear_extrapolation_and_linear_normal_interpolation(...)");
1045}
1046
1048 const tarch::la::Vector<Dimensions,int>& relativePositionWithinFatherCell,
1049 int numberOfDoFsPerAxisInPatch,
1050 int unknowns,
1051 const double* __restrict__ coarseGridCellValues,
1052 double* __restrict__ fineGridCellValues,
1053 bool extrapolateLinearly
1054) {
1055 static internal::CellInterpolationMap cellInterpolationMap;
1056 internal::CellInterpolationOperatorKey key(numberOfDoFsPerAxisInPatch);
1057
1058 tarch::multicore::Lock lock(interpolationMapSemaphore);
1059 if ( cellInterpolationMap.count(key)==0 ) {
1060 cellInterpolationMap[key] = internal::createLinearInterpolationMatrix(numberOfDoFsPerAxisInPatch,extrapolateLinearly);
1061 }
1062 tarch::la::DynamicMatrix* P = cellInterpolationMap[key];
1063 lock.free();
1064
1065 #if Dimensions==2
1066 const int patchSize = numberOfDoFsPerAxisInPatch*numberOfDoFsPerAxisInPatch;
1067 #else
1068 const int patchSize = numberOfDoFsPerAxisInPatch*numberOfDoFsPerAxisInPatch*numberOfDoFsPerAxisInPatch;
1069 #endif
1070
1071 P->batchedMultiplyAoS(
1072 fineGridCellValues, // image
1073 coarseGridCellValues, // preimage
1074 unknowns, // batch size, i.e. how often to apply it in one AoS rush
1075 patchSize, // result size, i.e. size of image
1076 serialiseMarkerIn3x3PatchAssembly(relativePositionWithinFatherCell, numberOfDoFsPerAxisInPatch)
1077 );
1078}
1079
1082 int numberOfDoFsPerAxisInPatch,
1083 int overlap,
1084 int unknowns,
1085 const double* __restrict__ normalInterpolationMatrix1d,
1086 const double* __restrict__ tangentialInterpolationMatrix1d,
1087 const double* __restrict__ coarseGridFaceValues,
1088 double* __restrict__ fineGridFaceValues
1089) {
1090 logTraceInWith3Arguments( "interpolateHaloLayer_AoS_tensor_product(...)", marker.toString(), numberOfDoFsPerAxisInPatch, overlap );
1091
1092 const int normal = marker.getSelectedFaceNumber() % Dimensions;
1093
1094 dfore(kFine,numberOfDoFsPerAxisInPatch,normal,0) {
1095 for (int iFine=0; iFine<overlap; iFine++) {
1097 dest(normal) = (marker.getSelectedFaceNumber() < Dimensions) ? iFine : 2 * overlap - iFine - 1;
1098
1099 int destLinearised = serialiseVoxelIndexInOverlap(
1100 dest,
1101 numberOfDoFsPerAxisInPatch, overlap, normal
1102 );
1103
1104 logDebug( "interpolateHaloLayer_AoS_tensor_product(...)", "clear dof " << dest );
1105 for (int unknown=0; unknown<unknowns; unknown++) {
1106 fineGridFaceValues[destLinearised*unknowns+unknown] = 0.0;
1107 }
1108
1109 dfore(kCoarse,numberOfDoFsPerAxisInPatch,normal,0) {
1110 for (int iCoarse=0; iCoarse<2*overlap; iCoarse++) {
1112 src(normal) = (marker.getSelectedFaceNumber() < Dimensions) ? iCoarse : 2 * overlap - iCoarse - 1;
1113 int srcLinearised = serialiseVoxelIndexInOverlap(
1114 src,
1115 numberOfDoFsPerAxisInPatch, overlap, normal
1116 );
1117
1118 double weight = 1.0;
1119 for (int d=0; d<Dimensions; d++) {
1120 if (d==normal) {
1121 int col = iCoarse;
1122 int row = iFine;
1123
1124 assertion4(col>=0,row,col,src,dest);
1125 assertion4(row>=0,row,col,src,dest);
1126 int index = col + row * 2 * overlap;
1127 logDebug( "interpolateHaloLayer_AoS_tensor_product(...)", "(" << row << "," << col << ")=" << index << " (normal contribution)");
1128 weight *= normalInterpolationMatrix1d[ index ];
1129 assertion(weight==weight);
1130 }
1131 else {
1132 int col = src(d);
1133 int row = dest(d) + marker.getRelativePositionWithinFatherCell()(d) * numberOfDoFsPerAxisInPatch;
1134 assertion4(col>=0,row,col,src,dest);
1135 assertion4(row>=0,row,col,src,dest);
1136 int index = col + row * numberOfDoFsPerAxisInPatch;
1137 logDebug( "interpolateHaloLayer_AoS_tensor_product(...)", "(" << row << "," << col << ")=" << index << " (tangential contribution)" );
1138 weight *= tangentialInterpolationMatrix1d[ index ];
1139 assertion(weight==weight);
1140 }
1141 }
1142
1143 logDebug( "interpolateHaloLayer_AoS_tensor_product(...)", "add dof " << src << " to " << dest << " with weight " << weight );
1144
1145 for (int unknown=0; unknown<unknowns; unknown++) {
1146 fineGridFaceValues[destLinearised*unknowns+unknown] += weight * coarseGridFaceValues[srcLinearised*unknowns+unknown];
1147 }
1148 }
1149 }
1150 }
1151 }
1152
1153 logTraceOut( "interpolateHaloLayer_AoS_tensor_product(...)" );
1154}
1155
1157 [[maybe_unused]] const peano4::datamanagement::FaceMarker& marker,
1158 [[maybe_unused]] int numberOfDoFsPerAxisInPatch,
1159 [[maybe_unused]] int overlap,
1160 [[maybe_unused]] int unknowns,
1161 [[maybe_unused]] const double* __restrict__ normalInterpolationMatrix1d,
1162 [[maybe_unused]] const double* __restrict__ tangentialInterpolationMatrix1d,
1163 [[maybe_unused]] const double* __restrict__ coarseGridCellValues,
1164 [[maybe_unused]] const double* __restrict__ coarseGridFaceValues,
1165 [[maybe_unused]] double* __restrict__ fineGridFaceValues
1166) {
1167 assertion(false);
1168}
1169
1171 [[maybe_unused]] const peano4::datamanagement::CellMarker& marker,
1172 [[maybe_unused]] int numberOfDoFsPerAxisInPatch,
1173 [[maybe_unused]] int unknowns,
1174 [[maybe_unused]] const double* __restrict__ interpolationMatrix1d,
1175 [[maybe_unused]] const double* __restrict__ coarseGridCellValues,
1176 [[maybe_unused]] double* __restrict__ fineGridCellValues
1177) {
1178 assertion(false);
1179}
#define assertion4(expr, param0, param1, param2, param3)
#define assertionEquals1(lhs, rhs, a)
#define assertionEquals(lhs, rhs)
#define assertion(expr)
And from this we can write down f$ nabla phi_i nabla phi_i dx but since we are constructing matrix let s investigate the f$ j
#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:642
#define dfore(counter, max, dim, value)
This is an exclusive d-dimensional for loop.
Definition Loop.h:937
#define dfor(counter, max)
d-dimensional Loop
Definition Loop.h:308
#define enddforx
I prefer to use this macro for dforx instead of a closing bracket as many syntax parser fail otherwis...
Definition Loop.h:928
My standard matrix is a matrix where the size is fixed at compile time.
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
#define endSimtDfor
Definition Loop.h:181
LoopPlacement
Guide loop-level parallelism.
Definition Loop.h:60
CPUGPUMethod int dLinearised(const tarch::la::Vector< Dimensions, int > &counter, int max)
Map d-dimensional vector onto integer index.
Definition Loop.cpp:106
bool greaterEquals(double lhs, double rhs, double tolerance=NUMERICAL_ZERO_DIFFERENCE)
bool smallerEquals(double lhs, double rhs, double tolerance=NUMERICAL_ZERO_DIFFERENCE)
Matrix< Rows, Cols, Scalar > multiplyComponents(const Matrix< Rows, X, Scalar > &lhs, const Matrix< X, Cols, Scalar > &rhs)
void freeMemory(void *data, MemoryLocation location)
@ Heap
Create data on the heap of the local device.
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)
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.
std::map< FaceInterpolationOperatorKey, tarch::la::DynamicMatrix * > FaceInterpolationMap
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.
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)
This interpolation should be used if a cell hosts two sets of unknowns.
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 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_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_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_linear_with_constant_extrapolation_and_linear_normal_interpolation(const peano4::datamanagement::CellMarker &marker, int numberOfDoFsPerAxisInPatch, int unknowns, const double *__restrict__ coarseGridCellValues, double *__restrict__ fineGridCellValues)
tarch::logging::Log _log("examples::unittests")
Provide information about selected face.
Definition FaceMarker.h:35
Simple vector class.
Definition Vector.h:134