Peano 4
Loading...
Searching...
No Matches
toolbox::blockstructured::internal Namespace Reference

Data Structures

struct  CellInterpolationOperatorKey
 
struct  FaceInterpolationOperatorKey
 

Typedefs

typedef std::map< FaceInterpolationOperatorKey, tarch::la::DynamicMatrix * > FaceInterpolationMap
 
typedef std::map< CellInterpolationOperatorKey, tarch::la::DynamicMatrix * > CellInterpolationMap
 A cell interpolation is a huge matrix.
 

Functions

tarch::la::DynamicMatrixcreateLinearInterpolationMatrix (int numberOfDoFsPerAxisInPatch, int normal, bool extrapolateLinearly, bool interpolateLinearlyAlongNormal)
 The realisation relies on the following observations/follows these steps:
 
tarch::la::DynamicMatrixcreateLinearInterpolationMatrix (int numberOfDoFsPerAxisInPatch, bool extrapolateLinearly)
 This is a volumetric version of the interpolation.
 
tarch::la::DynamicMatrixcreateLinearInterpolationMatrix (int numberOfDoFsPerAxisInInputPatch, int numberOfDoFsPerAxisInOutputPatch, bool extrapolateLinearly)
 
tarch::la::DynamicMatrixcreatePiecewiseConstantInterpolationMatrix (int numberOfDoFsPerAxisInPatch, int normal, int overlap)
 Create piecewise constant interpolation matrix.
 
tarch::la::DynamicMatrixcreatePiecewiseConstantInterpolationMatrix (int numberOfDoFsPerAxisInPatch)
 This is a volumetric version of the interpolation.
 
tarch::la::DynamicMatrixcreateInterpolationMatrixFrom1dTemplateByInsertingZeroColsAndRows (const tarch::la::DynamicMatrix &P1d, int numberOfDoFsPerAxisInPatch, int normal)
 This routine accepts a \( 3k \times k \) matrix as input.
 
tarch::la::DynamicMatrixcreateInterpolationMatrixFrom1dTemplateByLinearInterpolationAlongNormal (const tarch::la::DynamicMatrix &P1d, int numberOfDoFsPerAxisInPatch, int normal)
 This is an extension of createInterpolationMatrixFrom1dTemplateByInsertingZeroColsAndRows().
 
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.
 
void interpolateCell_AoS_piecewise_constant (const tarch::la::Vector< Dimensions, int > &relativePositionWithinFatherCell, int numberOfDoFsPerAxisInPatch, int unknowns, const double *__restrict__ fineGridValues, double *coarseGridValues)
 
void interpolateCell_AoS_linear (const tarch::la::Vector< Dimensions, int > &relativePositionWithinFatherCell, int numberOfDoFsPerAxisInPatch, int unknowns, const double *__restrict__ fineGridValues, double *coarseGridValues, bool extrapolateLinearly)
 
void clearHalfOfHaloLayerAoS (const peano4::datamanagement::FaceMarker &marker, int numberOfDoFsPerAxisInPatch, int overlap, int unknowns, bool clearInnerPart, double *values)
 Clear half of a halo layer.
 
void projectCells_AoS (const peano4::datamanagement::CellMarker &fineGridCellMarker, int numberOfDoFsPerAxisInPatch, std::function< void(tarch::la::Vector< Dimensions, int > coarseVolume, tarch::la::Vector< Dimensions, int > fineVolume, tarch::la::Vector< Dimensions, double > coarseVolumeCentre, tarch::la::Vector< Dimensions, double > fineVolumeCentre, double coarseVolumeH, double fineVolumeH)> update)
 Helper function.
 

Typedef Documentation

◆ CellInterpolationMap

A cell interpolation is a huge matrix.

If a patch consists of \( p \) entries per coordinate axis, then the matrix has the dimensions \( p^d \cdot 3^d \times p^d \). Its input is always the \( p^d \) vector of the coarse mesh. The output is typically a \( p^d \) vector as well, but there are \( 3^d \) of these guys.

We store one huge interpolation matrix and then pick the right band out of it.

Definition at line 867 of file Interpolation.h.

◆ FaceInterpolationMap

Function Documentation

◆ clearHalfOfHaloLayerAoS()

void toolbox::blockstructured::internal::clearHalfOfHaloLayerAoS ( const peano4::datamanagement::FaceMarker & marker,
int numberOfDoFsPerAxisInPatch,
int overlap,
int unknowns,
bool clearInnerPart,
double * values )

Clear half of a halo layer.

Faces in blockstructured codes with patches of size NxNxN typically carry overlaps of size MxNxN. The M is an even number with M<2N. If you use overlaps/halos of size 1, then M=2: a face manages one layer from the left and one layer from the right.

In blockstructured codes, you often have to erase half of this halo layer. If you run into an adaptive grid, for examples, you want to erase hanging layers before you interpolate or restrict into them. This is what this routine does.

Parameters
markerThis marker identifies a face. See the class documentation in particular for details about the enumeration (selected faces) of the faces from a cell's point of view.
numberOfDoFsPerAxisInPatchThe N in the description above
overlapEquals M/2 in the description above, i.e. if you each patch is surrounded by one halo layer, then the total overlap equals two and this argument equals 1.
unknownsNumber of unknowns that we store per voxel.
clearOuterPartIf this flag is set, we clear the outer half of the overlap (relative to the normal identified by marker).

Definition at line 451 of file Restriction.cpp.

References assertion, dfore, j, k, and toolbox::blockstructured::serialiseVoxelIndexInOverlap().

Here is the call graph for this function:

◆ createInterpolationMatrixFrom1dTemplateByInsertingZeroColsAndRows()

tarch::la::DynamicMatrix * toolbox::blockstructured::internal::createInterpolationMatrixFrom1dTemplateByInsertingZeroColsAndRows ( const tarch::la::DynamicMatrix & P1d,
int numberOfDoFsPerAxisInPatch,
int normal )

This routine accepts a \( 3k \times k \) matrix as input.

This is the 1d matrix. It clearly maps a sequence of k finite volumes onto a new sequence of 3k volumes.

We now exploit the fact any interpolation routine is basically a rotation or mirroring of this operator, while we have to take into account that we need two times k entries as we always store left and right halo overlaps.

Definition at line 55 of file Interpolation.cpp.

References assertionEquals1, tarch::la::DynamicMatrix::cols(), logDebug, tarch::la::DynamicMatrix::rows(), and tarch::la::DynamicMatrix::toString().

Referenced by createLinearInterpolationMatrix(), and createPiecewiseConstantInterpolationMatrix().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ createInterpolationMatrixFrom1dTemplateByLinearInterpolationAlongNormal()

tarch::la::DynamicMatrix * toolbox::blockstructured::internal::createInterpolationMatrixFrom1dTemplateByLinearInterpolationAlongNormal ( const tarch::la::DynamicMatrix & P1d,
int numberOfDoFsPerAxisInPatch,
int normal )

This is an extension of createInterpolationMatrixFrom1dTemplateByInsertingZeroColsAndRows().

In the setup as sketched below

we have two options to initialise bright green values: We can either use the bright red volumes only, or we can interpolate between the bright red ones and the dark red ones. This operation does the latter.

Definition at line 88 of file Interpolation.cpp.

References assertionEquals1, tarch::la::DynamicMatrix::cols(), logDebug, tarch::la::DynamicMatrix::rows(), and tarch::la::DynamicMatrix::toString().

Referenced by createLinearInterpolationMatrix().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ createLinearInterpolationMatrix() [1/3]

tarch::la::DynamicMatrix * toolbox::blockstructured::internal::createLinearInterpolationMatrix ( int numberOfDoFsPerAxisInInputPatch,
int numberOfDoFsPerAxisInOutputPatch,
bool extrapolateLinearly )

◆ createLinearInterpolationMatrix() [2/3]

tarch::la::DynamicMatrix * toolbox::blockstructured::internal::createLinearInterpolationMatrix ( int numberOfDoFsPerAxisInPatch,
bool extrapolateLinearly )

This is a volumetric version of the interpolation.

Definition at line 168 of file Interpolation.cpp.

References logDebug, logTraceInWith1Argument, logTraceOut, tarch::la::DynamicMatrix::replicateRows(), and tarch::la::DynamicMatrix::toString().

Here is the call graph for this function:

◆ createLinearInterpolationMatrix() [3/3]

tarch::la::DynamicMatrix * toolbox::blockstructured::internal::createLinearInterpolationMatrix ( int numberOfDoFsPerAxisInPatch,
int normal,
bool extrapolateLinearly,
bool interpolateLinearlyAlongNormal )

The realisation relies on the following observations/follows these steps:

A 1d setup

We first construct a 1d interpolation. Starting from the hard-coded pattern

{1.0/3.0, 2.0/3.0,     0.0},
{    0.0, 3.0/3.0,     0.0},
{    0.0, 2.0/3.0, 1.0/3.0}
   

we repeat this whole thing numberOfDoFsPerAxisInPatch times and each time shift by one. For numberOfDoFsPerAxisInPatch=4 this yields for example

{0.333333,0.666667,       0,0,0,0},
{0,              1,       0,0,0,0},
{0,       0.666667,0.333333,0,0,0},
{0,       0.333333,0.666667,0,0,0},
{0,              0,       1,0,0,0},
...
   

Let this be the matrix P1d. If we had a 1d setup, we could compute fine = P1d * coarse and would get all fine grid values in one rush. Obviously, we do not have a 1d setup.

2d extensions

Once we know P1d, we have to "blow is up". The input does not have numberOfDoFsPerAxisInPatch entries but has numberOfDoFsPerAxisInPatch*2 entries. Also the image has to be multiplied by two. We fist determine the normal and look whether we are left or right of this normal for the fine grid face. Next, we reiterate that we work with a lexicographic enumeration. If we want to study the left face of a cell,

The issue with the missing diagonal element

We do not know the diagonal element which we'd need to interpolate the outmost cells in our data. Without knowing them, two options are on the table: We can interpolate constantly at the edges, or we can extrapolate. While extrapolation might sound to be a good choice, I found that it sometimes yields physically invalid results. The Euler equations, for example, will immediately give you negative (tiny) energies or negative pressure.

The picture above illustrates the problem: We have by default no access to diagonal voxels in the mesh. This missing piece of information is illustrated by a voxel which is crossed out. If we now map the bright red voxels onto the bright green voxels, we don't really know what to do with the outermost voxels: We can extrapolate constantly or linearly, but we will always be slightly wrong, as we lack the information from the neighbour.

Extrapolating constantly seems to be pretty robust, but actually makes our interpolation scheme deterioriate towards something in-between piece-wise constant and linear. Extrapolating linearly is better, but might yield unphysical solutions. For example, assume there is a density shock with a steep drop between the two outermost voxels. If we extrapolate linearly, we might end up with an interpolated/extrapolated fine grid voxel with negative density. This yields a wrong solution. To fix such a case, you might want to work with linear extrapolation but also use a postprocessing (clean-up) step. See the guidebook for an explanation.

Definition at line 126 of file Interpolation.cpp.

References createInterpolationMatrixFrom1dTemplateByInsertingZeroColsAndRows(), createInterpolationMatrixFrom1dTemplateByLinearInterpolationAlongNormal(), logDebug, and tarch::la::DynamicMatrix::replicateRows().

Here is the call graph for this function:

◆ createPiecewiseConstantInterpolationMatrix() [1/2]

tarch::la::DynamicMatrix * toolbox::blockstructured::internal::createPiecewiseConstantInterpolationMatrix ( int numberOfDoFsPerAxisInPatch)

This is a volumetric version of the interpolation.

It is close to trivial, as we can employ internal::projectCells_AoS() to identify pairs of coarse-fine volumes that overlap and then we simply copy over results.

Definition at line 17 of file Interpolation.cpp.

References dfor, dfor3, peano4::utils::dLinearised(), enddforx, logDebug, and tarch::la::DynamicMatrix::toString().

Here is the call graph for this function:

◆ createPiecewiseConstantInterpolationMatrix() [2/2]

tarch::la::DynamicMatrix * toolbox::blockstructured::internal::createPiecewiseConstantInterpolationMatrix ( int numberOfDoFsPerAxisInPatch,
int normal,
int overlap )

Create piecewise constant interpolation matrix.

The matrix is very simple as it holds exclusively zeroes and ones. To construct it, we run over the image and and look for each image voxel which coarser face voxel is overlaps. In this context, we exploit the fact that all voxels are enumerated lexicographically:

The resulting matrix for the example above (where the normal equals 0) is 24x8 matrix which equals

\( \begin{array}{cccccccc} 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \end{array} \)

Below is another example for a face on y-direction in a 2D setup, assuming three-partitioning for every patch. when we refine a face, we get \( 3^{d-1} \) subfaces. we enumerate the volumes within a surface first, and then move to the next surface lexicographically.

Construction algorithm

I have two variants of the actual construction in here. One of them is commented out. It is the manual flavour which I wanted to keep to allow people to understand what's going on.

The one shipped is way more elegant (in my opinion) than the manual implementation. It starts from the observation that a 1d interpolation from a mesh onto a subdivided mesh with a subdivision factor by three is very simplistic matrix.

\[ \begin{array}{cccccccc} 1 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 1 \end{array} \]

Here, I use four-partitioning.

Formatting of documentation (doxygen)

Doxygen nowadays supports markup syntax. Unfortunately, this documentation would be indented more than five spaces which makes some of the write-up preformatted in markup. So the doxygen parser completely messes up things.

For this reason, the documentation here is not properly indented.

Definition at line 42 of file Interpolation.cpp.

References createInterpolationMatrixFrom1dTemplateByInsertingZeroColsAndRows(), and tarch::la::DynamicMatrix::replicateRows().

Here is the call graph for this function:

◆ interpolateCell_AoS_linear()

void toolbox::blockstructured::internal::interpolateCell_AoS_linear ( const tarch::la::Vector< Dimensions, int > & relativePositionWithinFatherCell,
int numberOfDoFsPerAxisInPatch,
int unknowns,
const double *__restrict__ fineGridValues,
double * coarseGridValues,
bool extrapolateLinearly )

◆ interpolateCell_AoS_piecewise_constant()

void toolbox::blockstructured::internal::interpolateCell_AoS_piecewise_constant ( const tarch::la::Vector< Dimensions, int > & relativePositionWithinFatherCell,
int numberOfDoFsPerAxisInPatch,
int unknowns,
const double *__restrict__ fineGridValues,
double * coarseGridValues )
See also
Other interpolateCell_AoS_piecewise_constant() variant which delegates to this one which is also used by the face interpolation.

◆ projectCells_AoS()

void toolbox::blockstructured::internal::projectCells_AoS ( const peano4::datamanagement::CellMarker & fineGridCellMarker,
int numberOfDoFsPerAxisInPatch,
std::function< void(tarch::la::Vector< Dimensions, int > coarseVolume, tarch::la::Vector< Dimensions, int > fineVolume, tarch::la::Vector< Dimensions, double > coarseVolumeCentre, tarch::la::Vector< Dimensions, double > fineVolumeCentre, double coarseVolumeH, double fineVolumeH)> update )

Helper function.

This function runs through all fine grid/coarse grid cell combinations which are identified via fineGridCellMarker. So we call the functor once per fine grid cell, but \( 3^d \) times or not at all for any coarse grid cell. This routine is used as a helper function for piece-wise constant interpolation. It is not that useful for d-linear interpolation.

Definition at line 481 of file Restriction.cpp.

References dfor, logTraceInWith2Arguments, logTraceOut, and tarch::la::multiplyComponents().

Here is the call graph for this function:

◆ projectInterpolatedFineCellsOnHaloLayer_AoS()

void toolbox::blockstructured::internal::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.

We identify the overlap with the fine halo layer and copy values over from these two arrays onto the face data.

Definition at line 205 of file Interpolation.cpp.

References dfore, peano4::utils::dLinearised(), j, logTraceInWith3Arguments, logTraceOut, and toolbox::blockstructured::serialiseVoxelIndexInOverlap().

Here is the call graph for this function: