Peano 4

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::DynamicMatrix *  createLinearInterpolationMatrix (int numberOfDoFsPerAxisInPatch, int normal, bool extrapolateLinearly, bool interpolateLinearlyAlongNormal) 
The realisation relies on the following observations/follows these steps:  
tarch::la::DynamicMatrix *  createLinearInterpolationMatrix (int numberOfDoFsPerAxisInPatch, bool extrapolateLinearly) 
This is a volumetric version of the interpolation.  
tarch::la::DynamicMatrix *  createLinearInterpolationMatrix (int numberOfDoFsPerAxisInInputPatch, int numberOfDoFsPerAxisInOutputPatch, bool extrapolateLinearly) 
tarch::la::DynamicMatrix *  createPiecewiseConstantInterpolationMatrix (int numberOfDoFsPerAxisInPatch, int normal, int overlap) 
Create piecewise constant interpolation matrix.  
tarch::la::DynamicMatrix *  createPiecewiseConstantInterpolationMatrix (int numberOfDoFsPerAxisInPatch) 
This is a volumetric version of the interpolation.  
tarch::la::DynamicMatrix *  createInterpolationMatrixFrom1dTemplateByInsertingZeroColsAndRows (const tarch::la::DynamicMatrix &P1d, int numberOfDoFsPerAxisInPatch, int normal) 
This routine accepts a \( 3k \times k \) matrix as input.  
tarch::la::DynamicMatrix *  createInterpolationMatrixFrom1dTemplateByLinearInterpolationAlongNormal (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 std::map< CellInterpolationOperatorKey, tarch::la::DynamicMatrix* > toolbox::blockstructured::internal::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.
typedef std::map< FaceInterpolationOperatorKey, tarch::la::DynamicMatrix* > toolbox::blockstructured::internal::FaceInterpolationMap 
Definition at line 857 of file Interpolation.h.
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.
marker  This 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. 
numberOfDoFsPerAxisInPatch  The N in the description above 
overlap  Equals 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. 
unknowns  Number of unknowns that we store per voxel. 
clearOuterPart  If 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().
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().
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().
tarch::la::DynamicMatrix * toolbox::blockstructured::internal::createLinearInterpolationMatrix  (  int  numberOfDoFsPerAxisInInputPatch, 
int  numberOfDoFsPerAxisInOutputPatch,  
bool  extrapolateLinearly ) 
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().
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:
We first construct a 1d interpolation. Starting from the hardcoded 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.
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,
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 inbetween piecewise 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 (cleanup) step. See the guidebook for an explanation.
Definition at line 126 of file Interpolation.cpp.
References createInterpolationMatrixFrom1dTemplateByInsertingZeroColsAndRows(), createInterpolationMatrixFrom1dTemplateByLinearInterpolationAlongNormal(), logDebug, and tarch::la::DynamicMatrix::replicateRows().
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 coarsefine 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().
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 ydirection in a 2D setup, assuming threepartitioning for every patch. when we refine a face, we get \( 3^{d1} \) subfaces. we enumerate the volumes within a surface first, and then move to the next surface lexicographically.
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 fourpartitioning.
Doxygen nowadays supports markup syntax. Unfortunately, this documentation would be indented more than five spaces which makes some of the writeup 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().
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 ) 
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 ) 
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 piecewise constant interpolation. It is not that useful for dlinear interpolation.
Definition at line 481 of file Restriction.cpp.
References dfor, logTraceInWith2Arguments, logTraceOut, and tarch::la::multiplyComponents().
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().