![]() |
Peano
|
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 1039 of file Interpolation.h.
typedef std::map<FaceInterpolationOperatorKey, tarch::la::DynamicMatrix*> toolbox::blockstructured::internal::FaceInterpolationMap |
Definition at line 1026 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 1310 of file Restriction.cpp.
References assertion, dfore, 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 73 of file Interpolation.cpp.
References assertionEquals1, tarch::la::DynamicMatrix::cols(), logDebug, P, 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 110 of file Interpolation.cpp.
References assertionEquals1, tarch::la::DynamicMatrix::cols(), logDebug, P, 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 214 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 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.
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 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 152 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 coarse-fine volumes that overlap and then we simply copy over results.
Definition at line 20 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 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.
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.
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 57 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 piece-wise constant interpolation. It is not that useful for d-linear interpolation.
Definition at line 1341 of file Restriction.cpp.
References tarch::la::convertScalar(), 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 277 of file Interpolation.cpp.
References dfore, peano4::utils::dLinearised(), logTraceInWith3Arguments, logTraceOut, and toolbox::blockstructured::serialiseVoxelIndexInOverlap().