Peano
Loading...
Searching...
No Matches
toolbox::finiteelements Namespace Reference

Derive Element-wise Matrices From Stencils. More...

Namespaces

namespace  tests
 

Data Structures

class  BSplinesStencilFactory
 Stencil Factory. More...
 

Typedefs

typedef tarch::la::Matrix< TwoPowerD, TwoPowerD, doubleElementWiseAssemblyMatrix
 Element-wise Matrix.
 
typedef tarch::la::Matrix< TwoPowerD, TwoPowerD, std::complex< double > > ComplexElementWiseAssemblyMatrix
 
typedef tarch::la::Matrix< TwoPowerD, TwoPowerD, doubleInterGridTransferMatrix
 
typedef tarch::la::Matrix< TwoPowerD/2, TwoPowerD/2, doubleElementWiseAssemblyMatrixOnSpaceTimeGrid
 
typedef tarch::la::Vector< ThreePowerD, doubleStencil
 Stencil.
 
typedef tarch::la::Vector< ThreePowerD, std::complex< double > > ComplexStencil
 
typedef tarch::la::Vector< ThreePowerD *TwoPowerD, doubleVectorOfStencils
 Vectors.
 
typedef tarch::la::Vector< ThreePowerD *TwoPowerD, std::complex< double > > VectorOfComplexStencils
 
typedef tarch::la::Vector< TwoPowerD, doubleElementWiseVector
 Vector of Element (Lexicographic Ordering of Unknowns)
 

Functions

ElementWiseAssemblyMatrix getElementWiseAssemblyMatrix (const Stencil &stencil)
 Derive element-wise matrix from stencil.
 
ComplexElementWiseAssemblyMatrix getElementWiseAssemblyMatrix (const ComplexStencil &complexStencil)
 Derive element-wise matrix from stencil.
 
Stencil reconstructUniformStencilFragments (const ElementWiseAssemblyMatrix &matrix)
 Reconstruct stencils from element matrices.
 
tarch::la::Vector< TwoPowerD *ThreePowerD, doublereconstructStencilFragments (const ElementWiseAssemblyMatrix &matrix)
 Reconstruct stencils from element matrices.
 
ElementWiseAssemblyMatrix getElementWiseAssemblyMatrix (const VectorOfStencils &vectorOfStencils)
 
ElementWiseAssemblyMatrix getElementWiseAssemblyMatrix (const VectorOfStencils &vectorOfStencils, const std::bitset< ThreePowerD > &cellIsInside)
 I originally thought it would be possible to identify how many neighbours there are from the adjacent vertices' state.
 
ComplexElementWiseAssemblyMatrix getElementWiseAssemblyMatrix (const VectorOfComplexStencils &vectorOfStencils, const std::bitset< ThreePowerD > &cellIsInside)
 See non-complex variant.
 
int mapElementMatrixEntryOntoStencilEntry (int row, int col)
 An element-wise assembly matrix is distilled from a sequence of \( 2^d \) stencils which are in turn \( 3^d \) double arrays.
 
int mapElementMatrixEntryOntoStencilEntry (const tarch::la::Vector< Dimensions, int > &row, const tarch::la::Vector< Dimensions, int > &col)
 
double getDiagonalElement (const ElementWiseAssemblyMatrix &matrix)
 
double getDiagonalElement (const Stencil &stencil)
 
template<int StencilSize>
tarch::la::Vector< StencilSize *StencilSize, doublestencilProduct (const tarch::la::Vector< StencilSize, double > &a, const tarch::la::Vector< StencilSize, double > &b)
 
template<int StencilSize>
tarch::la::Vector< StencilSize *StencilSize *StencilSize, doublestencilProduct (const tarch::la::Vector< StencilSize, double > &a, const tarch::la::Vector< StencilSize, double > &b, const tarch::la::Vector< StencilSize, double > &c)
 
ElementWiseAssemblyMatrix hierarchicalTransform (const ElementWiseAssemblyMatrix &matrix, const tarch::la::Vector< Dimensions, double > &h, double weight=1.0)
 
ElementWiseAssemblyMatrix inverseHierarchicalTransform (const ElementWiseAssemblyMatrix &matrix, const tarch::la::Vector< Dimensions, double > &h, double weight=1.0)
 
int getStencilEntryLinearisedIndex (const tarch::la::Vector< Dimensions, int > stencilEntry)
 
void preprocessBoundaryStencil (Stencil &stencil, const std::bitset< Dimensions *2 > &boundaryFaceNormals)
 Preprocess boundary stencils held by boundary vertices.
 
Stencil exchangeCoordinates (const Stencil &stencil, int coord0, int coord1)
 Exchanges the coordinates of a stencil.
 
tarch::la::Vector< 3, doubleget1DMassStencil ()
 
tarch::la::Vector< 3, doubleget1DLaplaceStencil ()
 Stencil is not scaled at all with any mesh width.
 
tarch::la::Vector< 3, doubleget1DIdentity ()
 Stencil is not scaled at all with any mesh width.
 
tarch::la::Vector< 3, doubleget1DLeftIdentity ()
 
tarch::la::Vector< 3, doubleget1DRightIdentity ()
 
tarch::la::Vector< 5, doubleget1DLinearInterpolationStencil ()
 Stencil is not scaled at all with any mesh width.
 
tarch::la::Vector< 3, doubleget1DMeanValueStencil ()
 Stencil is not scaled at all with any mesh width.
 
tarch::la::Vector< 3, doubleget1DDownwindStencil ()
 Stencil is not scaled at all with any mesh width.
 
tarch::la::Vector< 3, doubleget1DUpwindStencil ()
 
tarch::la::Vector< 3, doubleget1DCentralDifferencesStencil ()
 
tarch::la::Vector< 3 *3, doublestencilProduct (const tarch::la::Vector< 3, double > &a, const tarch::la::Vector< 3, double > &b)
 Makes stencil-product of two stencils:
 
tarch::la::Vector< 5 *5, doublestencilProduct (const tarch::la::Vector< 5, double > &a, const tarch::la::Vector< 5, double > &b)
 
tarch::la::Vector< 3 *3 *3, doublestencilProduct (const tarch::la::Vector< 3, double > &a, const tarch::la::Vector< 3, double > &b, const tarch::la::Vector< 3, double > &c)
 Equals a * (b*c)
 
tarch::la::Vector< 5 *5 *5, doublestencilProduct (const tarch::la::Vector< 5, double > &a, const tarch::la::Vector< 5, double > &b, const tarch::la::Vector< 5, double > &c)
 
tarch::la::Vector< 3 *3 *3 *3, doublestencilProduct (const tarch::la::Vector< 3, double > &a, const tarch::la::Vector< 3, double > &b, const tarch::la::Vector< 3, double > &c, const tarch::la::Vector< 3, double > &d)
 
tarch::la::Vector< 5 *5 *5 *5, doublestencilProduct (const tarch::la::Vector< 5, double > &a, const tarch::la::Vector< 5, double > &b, const tarch::la::Vector< 5, double > &c, const tarch::la::Vector< 5, double > &d)
 
tarch::la::Vector< 3 *3 *3 *3 *3, doublestencilProduct (const tarch::la::Vector< 3, double > &a, const tarch::la::Vector< 3, double > &b, const tarch::la::Vector< 3, double > &c, const tarch::la::Vector< 3, double > &d, const tarch::la::Vector< 3, double > &e)
 
tarch::la::Vector< 5 *5 *5 *5 *5, doublestencilProduct (const tarch::la::Vector< 5, double > &a, const tarch::la::Vector< 5, double > &b, const tarch::la::Vector< 5, double > &c, const tarch::la::Vector< 5, double > &d, const tarch::la::Vector< 5, double > &e)
 
tarch::la::Vector< ThreePowerD, doublegetLaplacian (const tarch::la::Vector< Dimensions, double > &h, const tarch::la::Vector< Dimensions, double > &scaling=1.0)
 Computes the Laplacian.
 
tarch::la::Vector< ThreePowerD, std::complex< double > > getLaplacian (const tarch::la::Vector< Dimensions, std::complex< double > > &h, const tarch::la::Vector< Dimensions, std::complex< double > > &scaling=std::complex< double >(1.0))
 
tarch::la::Vector< ThreePowerD, std::complex< double > > getLaplacian (const tarch::la::Vector< Dimensions, double > &h, const tarch::la::Vector< Dimensions, std::complex< double > > &scaling)
 
tarch::la::Vector< ThreePowerD, doublegetUpwindDiscretisedConvection (const tarch::la::Vector< Dimensions, double > &h, const tarch::la::Vector< Dimensions, double > &convCoeff)
 Computes a convection operator, defined by the convection coefficients, discretised with the first order upwind scheme.
 
void applyNeumannBC (tarch::la::Vector< ThreePowerD, double > &stencil, double &rhs, const std::bitset< Dimensions *2 > &boundaryFaceNormals, const tarch::la::Vector< Dimensions, double > &h, double derivative)
 
tarch::la::Vector< ThreePowerD, doublegetMassMatrix (const tarch::la::Vector< Dimensions, double > &h)
 Computes the mass matrix.
 
tarch::la::Vector< ThreePowerD, std::complex< double > > getMassMatrix (const tarch::la::Vector< Dimensions, std::complex< double > > &h)
 
tarch::la::Vector< ThreePowerD, std::complex< double > > getHelmholtzShiftMatrix (const tarch::la::Vector< Dimensions, double > &h, const std::complex< double > &phi)
 Computes the Helmholtz shift term matrix.
 
tarch::la::Vector< ThreePowerD, std::complex< double > > getHelmholtzShiftMatrix (const tarch::la::Vector< Dimensions, std::complex< double > > &h, const std::complex< double > &phi)
 
tarch::la::Vector< ThreePowerD, doublegetIdentity (const tarch::la::Vector< Dimensions, double > &h)
 Computes the mass matrix.
 
tarch::la::Vector< ThreePowerD, std::complex< double > > getIdentity (const tarch::la::Vector< Dimensions, std::complex< double > > &h)
 
tarch::la::Vector< FivePowerD, doublegetDLinearInterpolation ()
 Computes the stencil for a d-linear interpolation.
 
Stencil extractElementStencil (const Stencil &stencil, const tarch::la::Vector< Dimensions, int > &adjacentCell)
 The standard stencil spans 2^d cells adjacent to a vertex.
 
ElementWiseAssemblyMatrix getPoissonMatrixWithJumpingCoefficient (const tarch::la::Vector< Dimensions, double > &cellCentre, const tarch::la::Vector< Dimensions, double > &h, int integrationPointsPerAxis, std::function< double(const tarch::la::Vector< Dimensions, double > &)>)
 Integrate over element with a subvoxel scheme (cmp.
 
tarch::tests::TestCasegetUnitTests ()
 Please destroy after usage.
 

Detailed Description

Derive Element-wise Matrices From Stencils.

Stencil factory methods.

If you build something with this toolbox, please link vs.

Many components in the Peano framework switch from a stencil notation to element-wise assembly matrices frequently. As both representations (stencils and element-wise matrices) are equivalent, it is reasonable to code only one representation, and to derive the alternative representation automatically. This class provides the required operation.

  • Stencils are stored within vertices (red lines for a 5-point stencil in 2d, e.g.).
  • However, sometimes we don't have access to neighbouring cells (big blue cell, e.g., is one coarse cell where the solver automaton descends to finer grids)
  • Consequently, we have to convert the \( 2^d \) vertex stencils within each cell into a element-wise assembly matrix \( \mathbf{R}^{2^d \times 2^d}\). With this matrix, we then can accumulate the residuals.
  • PeProt automatically generates a routine that takes the \( 2^d \) vertices and writes their \( 3^d \) stencils into one big vector of length \( 2^d \cdot 3^d \).
  • This class provides a routine that transforms this big vector into an element assembly matrix.
  • There's also a shortcut operaration that takes only one stencil. In this case, we assume that the stencil is the same for all \( 2^d \) adjacent vertices of the cell.

Example

We have the stencil

\( \left[ \begin{array}{ccc} -1 & -1 & -1 \\ -1 & 8 & -1 \\ -1 & -1 & -1 \end{array} \right] \)

for (d=2). It is stored as a vector of length nine. Within a cell, we have four vertices. The values u of these vertices can be converted into one vector of length four due to a call of the operation

tarch::la::Vector<TWO_POWER_D,double> u = MyVertices::readU(...);
float u
Simple vector class.
Definition Vector.h:134

Such a operation is generated, if you add the statement

read scalar: U

to your Peano specification file. Now, if you wanna compute matrix times vector, you basically compute $Au$ with A being a \( 4 \times 4 \) matrix. A is well-defined by the stencil. In this example, it would be

\( \left[ \begin{array}{cccc} 2 & -0.5 & -0.5 & -1.0 \\ -0.5 & 2 & -1.0 & -0.5 \\ -0.5 & -1.0 & 2 & -0.5 \\ -1.0 & -0.5 & -0.5 & 2 \end{array} \right] \)

The methods of this class convert your stencil into the element-wise stiffness matrix. There's two variants of conversion operations:

  • The simple variant takes one stencil and creates the stiffness matrix. It resembles the example.
  • A more flexible variant takes \( 2^d \) stencils and creates the element-wise system matrix. This is the method of choice if you have different stencils on the different vertices.
Author
Tobias Weinzierl

ToolboxFiniteElementxdapp where x is 2 or 2 and app is either nothing, or _debug or _trace.

Author
Tobias Weinzierl

Many components in the Peano framework switch from a stencil notation to element-wise assembly matrices. This class provides some operations to construct the stencils. To convert these stencil into element-wise assembly matrices, use the class ElementMatrix.

Author
Tobias Weinzierl

Typedef Documentation

◆ ComplexElementWiseAssemblyMatrix

◆ ComplexStencil

Definition at line 30 of file Stencil.h.

◆ ElementWiseAssemblyMatrix

◆ ElementWiseAssemblyMatrixOnSpaceTimeGrid

◆ ElementWiseVector

Vector of Element (Lexicographic Ordering of Unknowns)

Definition at line 41 of file Stencil.h.

◆ InterGridTransferMatrix

◆ Stencil

Stencil.

Definition at line 29 of file Stencil.h.

◆ VectorOfComplexStencils

◆ VectorOfStencils

Function Documentation

◆ applyNeumannBC()

void toolbox::finiteelements::applyNeumannBC ( tarch::la::Vector< ThreePowerD, double > & stencil,
double & rhs,
const std::bitset< Dimensions *2 > & boundaryFaceNormals,
const tarch::la::Vector< Dimensions, double > & h,
double derivative )
    Eliminate Neumann BC by adapting the stencil

    Please note that you often have to call preprocessBoundaryStencil() on
    the modified stencil afterwards.

    <h1>Realisation</h1>

    The realisation is very simple and relies on a central finite
    difference discretisation. This is the discretisation corresponding to
    d-linear finite elements. We explain the discretisation at hands of a
    2d setup for a vertex that is at the right boundary of the domain. The
    stencil is given by
    <pre>

s6 s7 s8 -1 -1 -1 s3 s4 s5 = -1 8 -1 s0 s1 s2 -1 -1 -1

We (virtually) extend the domain by one layer and obtain an additional virtual degree of freedom right of the boundary. Let's denote all the unknowns accessed by the stencil as

u6  u7  u8
u3  u4  u5
u0  u1  u2
        

Obviously, u2,u5 and u8 are not defined. However, we now apply the Neumann conditions with the rules:

u2 = h(0)*derivative + u0 u5 = h(0)*derivative + u3 u8 = h(0)*derivative + u6

We see that the residual evaluation results from an evaluation of the u0,u1,u3,u4,u6,u7 terms plus an additional term

h(0)*derivative*(s2+s5+s8) + s2*u0 + s5*u3 + s8*u6

This is what the routine does: it modifies the stencil (here s0 \mapsto s0+s2, s3 \mapsto s3+s5 and s6 \mapsto s6+s8) and returns the modified rhs.

Parameters
hcell sizes around the vertex
derivativePass zero if you have homogeneous Neumann

Definition at line 558 of file StencilFactory.cpp.

◆ exchangeCoordinates()

toolbox::finiteelements::Stencil toolbox::finiteelements::exchangeCoordinates ( const Stencil & stencil,
int coord0,
int coord1 )

Exchanges the coordinates of a stencil.

Operations like this operation can be used to rotate a stencil, i.e. to derive several stencils from one input stencil.

Definition at line 22 of file StencilFactory.cpp.

References assertionMsg, d3for3, and enddforx.

◆ extractElementStencil()

toolbox::finiteelements::Stencil toolbox::finiteelements::extractElementStencil ( const Stencil & stencil,
const tarch::la::Vector< Dimensions, int > & adjacentCell )

The standard stencil spans 2^d cells adjacent to a vertex.

Let (0,0) be the element left or below the current vertex, (1,0) is the one below and right of the vertex, and so forth. The stencil decomposes over all adjacent cells. Such a decomposition is used in getElementWiseAssemblyMatrix() for example to facilitate element-wise matrix-vector products.

This operation extracts from the stencil only the fragment that 'belongs' to a certain adjacent cell. It is thus similar to getElementWiseAssemblyMatrix() but realises everything from a stencil point of view rather than with an element-wise language.

Definition at line 831 of file StencilFactory.cpp.

References dfor2, peano4::utils::dLinearised(), and enddforx.

Here is the call graph for this function:

◆ get1DCentralDifferencesStencil()

tarch::la::Vector< 3, double > toolbox::finiteelements::get1DCentralDifferencesStencil ( )

Definition at line 106 of file StencilFactory.cpp.

◆ get1DDownwindStencil()

tarch::la::Vector< 3, double > toolbox::finiteelements::get1DDownwindStencil ( )

Stencil is not scaled at all with any mesh width.

Definition at line 88 of file StencilFactory.cpp.

◆ get1DIdentity()

tarch::la::Vector< 3, double > toolbox::finiteelements::get1DIdentity ( )

Stencil is not scaled at all with any mesh width.

Definition at line 61 of file StencilFactory.cpp.

Referenced by getIdentity(), and getIdentity().

Here is the caller graph for this function:

◆ get1DLaplaceStencil()

tarch::la::Vector< 3, double > toolbox::finiteelements::get1DLaplaceStencil ( )

Stencil is not scaled at all with any mesh width.

Returns
\( [-1, 2, -1] \)

Definition at line 124 of file StencilFactory.cpp.

Referenced by getLaplacian(), getLaplacian(), and getLaplacian().

Here is the caller graph for this function:

◆ get1DLeftIdentity()

tarch::la::Vector< 3, double > toolbox::finiteelements::get1DLeftIdentity ( )

Definition at line 70 of file StencilFactory.cpp.

◆ get1DLinearInterpolationStencil()

tarch::la::Vector< 5, double > toolbox::finiteelements::get1DLinearInterpolationStencil ( )

Stencil is not scaled at all with any mesh width.

Definition at line 41 of file StencilFactory.cpp.

Referenced by getDLinearInterpolation().

Here is the caller graph for this function:

◆ get1DMassStencil()

tarch::la::Vector< 3, double > toolbox::finiteelements::get1DMassStencil ( )
Returns
\( [\frac{1}{6}, \frac{2}{3}, \frac{1}{6}] \)

Definition at line 52 of file StencilFactory.cpp.

Referenced by getHelmholtzShiftMatrix(), getHelmholtzShiftMatrix(), getLaplacian(), getLaplacian(), getLaplacian(), getMassMatrix(), and getMassMatrix().

Here is the caller graph for this function:

◆ get1DMeanValueStencil()

tarch::la::Vector< 3, double > toolbox::finiteelements::get1DMeanValueStencil ( )

Stencil is not scaled at all with any mesh width.

Definition at line 115 of file StencilFactory.cpp.

◆ get1DRightIdentity()

tarch::la::Vector< 3, double > toolbox::finiteelements::get1DRightIdentity ( )

Definition at line 79 of file StencilFactory.cpp.

◆ get1DUpwindStencil()

tarch::la::Vector< 3, double > toolbox::finiteelements::get1DUpwindStencil ( )

Definition at line 97 of file StencilFactory.cpp.

◆ getDiagonalElement() [1/2]

double toolbox::finiteelements::getDiagonalElement ( const ElementWiseAssemblyMatrix & matrix)

Definition at line 59 of file ElementMatrix.cpp.

References TwoPowerD.

◆ getDiagonalElement() [2/2]

double toolbox::finiteelements::getDiagonalElement ( const Stencil & stencil)

Definition at line 64 of file ElementMatrix.cpp.

References ThreePowerD.

◆ getDLinearInterpolation()

tarch::la::Vector< FivePowerD, double > toolbox::finiteelements::getDLinearInterpolation ( )

Computes the stencil for a d-linear interpolation.

Definition at line 791 of file StencilFactory.cpp.

References assertionMsg, get1DLinearInterpolationStencil(), and stencilProduct().

Here is the call graph for this function:

◆ getElementWiseAssemblyMatrix() [1/5]

ComplexElementWiseAssemblyMatrix toolbox::finiteelements::getElementWiseAssemblyMatrix ( const ComplexStencil & complexStencil)

Derive element-wise matrix from stencil.

See also
namespace description
matrixfree::stencil::preprocessBoundaryStencil() in StencilFactory.h
getElementWiseAssemblyMatrix( const VectorOfStencils& vectorOfStencils )
Todo
Die Abbildung sollte man im Konstruktor einmal bauen und dann hier nur noch anwenden. Deshalb ist das Ding ja eine Methode und kein statisches Ding.

Definition at line 35 of file ElementMatrix.cpp.

References dfor2, peano4::utils::dLinearised(), and enddforx.

Here is the call graph for this function:

◆ getElementWiseAssemblyMatrix() [2/5]

toolbox::finiteelements::ElementWiseAssemblyMatrix toolbox::finiteelements::getElementWiseAssemblyMatrix ( const Stencil & stencil)

Derive element-wise matrix from stencil.

The operation takes the stencils from the adjacent vertices and derives the parts of them that can be evaluated in one cell. See the first row in the picture below.

The nine-point stencil is split up between the four adjacent cells. We split up all entries that could be evaluated in multiple cells equally. The diagonal element, e.g., can be evaluated by all four adjacent cells. Therefore, we evaluate it in every cell but divide the contributions by four. If we sum up, we end up with the original stencil.

This whole arguing fails if we evaluate stencils along the boundary and if there are unknowns located on the boundary. See the second and third row or preprocessBoundaryStencil() for a detailed discussion.

See also
namespace description
matrixfree::stencil::preprocessBoundaryStencil() in StencilFactory.h
Todo
Die Abbildung sollte man im Konstruktor einmal bauen und dann hier nur noch anwenden. Deshalb ist das Ding ja eine Methode und kein statisches Ding.

Definition at line 10 of file ElementMatrix.cpp.

References dfor2, peano4::utils::dLinearised(), and enddforx.

Referenced by getElementWiseAssemblyMatrix(), getPoissonMatrixWithJumpingCoefficient(), and hierarchicalTransform().

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

◆ getElementWiseAssemblyMatrix() [3/5]

toolbox::finiteelements::ComplexElementWiseAssemblyMatrix toolbox::finiteelements::getElementWiseAssemblyMatrix ( const VectorOfComplexStencils & vectorOfStencils,
const std::bitset< ThreePowerD > & cellIsInside )

See non-complex variant.

Definition at line 192 of file ElementMatrix.cpp.

References assertionMsg.

◆ getElementWiseAssemblyMatrix() [4/5]

toolbox::finiteelements::ElementWiseAssemblyMatrix toolbox::finiteelements::getElementWiseAssemblyMatrix ( const VectorOfStencils & vectorOfStencils)
   Derive element-wise matrix from @f$ 2^d @f$ stencils, i.e. this operation
   takes one big @f$ 2^d \cdot 3^d @f$ vector as argument and returns an
   @f$ 2^d \times 2^d @f$ matrix.

   Typical usage if you have data on the heap:
   <pre>

VectorOfStencils stencils; dfor2(k) tarch::la::slice( stencils, // into stencils tarch::la::Vector<THREE_POWER_D,double>(fineGridVertices[ fineGridVerticesEnumerator(k) ].getStencil()), kScalar*THREE_POWER_D ); enddforx

Please note that this operation does not work properly at the domain boundary: We do divide the contributions among all adjacent cells (see the sketch below; only the upper row, the other two rows are used to illustrate something else).

Obviously, this splitting is invalid along the domain boundary, where fewer cells are adjacent and thus other scalings have to be used. There's an extended/overloaded version of this routine to handle those cases.

See also
getElementWiseAssemblyMatrix( const Stencil& stencil );

Definition at line 100 of file ElementMatrix.cpp.

References assertionNumericalEquals2, dfor2, peano4::utils::dLinearised(), enddforx, and ThreePowerD.

Here is the call graph for this function:

◆ getElementWiseAssemblyMatrix() [5/5]

toolbox::finiteelements::ElementWiseAssemblyMatrix toolbox::finiteelements::getElementWiseAssemblyMatrix ( const VectorOfStencils & vectorOfStencils,
const std::bitset< ThreePowerD > & cellIsInside )

I originally thought it would be possible to identify how many neighbours there are from the adjacent vertices' state.

However, this is not unique. If we look at the following 3x3 cells (outside and inside)

o o i o o i i i i

and on the two vertices on the right (there are four real vertices in this ASCII art that are really surrounded by known cells). The upper one has a surrounding pattern of

o i o i

and thus two inner adjacent cells. In the upper right inner cell, we know that it is boundary. In the cell below, the lower left vertex has exactly the same pattern from inside the cell. However, as its adjacent cells have the state

o i i i

this one has three inner cells and, thus, is to be scaled differently.

Transition from other operation

If you want to replace your call to getElementWiseAssemblyMatrix() without any inside/outside flags, you simply have to add

std::bitset<THREE_POWER_D>(0).flip()

as additional argument.

Parameters
cellIsOutsideIdentifies which of the 3^d-1 surrounding cells are inside.

Definition at line 151 of file ElementMatrix.cpp.

References assertion2, assertion4, dfor2, peano4::utils::dLinearised(), enddforx, tarch::la::equals(), getElementWiseAssemblyMatrix(), and ThreePowerD.

Here is the call graph for this function:

◆ getHelmholtzShiftMatrix() [1/2]

tarch::la::Vector< ThreePowerD, std::complex< double > > toolbox::finiteelements::getHelmholtzShiftMatrix ( const tarch::la::Vector< Dimensions, double > & h,
const std::complex< double > & phi )

Computes the Helmholtz shift term matrix.

The Helmholtz shift is the term that needs to be added to the Laplacian matrix in order to get the full Helmholtz matrix

Definition at line 663 of file StencilFactory.cpp.

References assertionMsg, get1DMassStencil(), and stencilProduct().

Here is the call graph for this function:

◆ getHelmholtzShiftMatrix() [2/2]

tarch::la::Vector< ThreePowerD, std::complex< double > > toolbox::finiteelements::getHelmholtzShiftMatrix ( const tarch::la::Vector< Dimensions, std::complex< double > > & h,
const std::complex< double > & phi )

Definition at line 686 of file StencilFactory.cpp.

References assertionMsg, get1DMassStencil(), and stencilProduct().

Here is the call graph for this function:

◆ getIdentity() [1/2]

tarch::la::Vector< ThreePowerD, double > toolbox::finiteelements::getIdentity ( const tarch::la::Vector< Dimensions, double > & h)

Computes the mass matrix.

The mass matrix is returned in the finite element form, i.e.~the entries result from an integral over a h. They are scaled with the volume of h. As such it corresponds to the integral over a characteristic function of one cell.

Definition at line 726 of file StencilFactory.cpp.

References assertionMsg, get1DIdentity(), and stencilProduct().

Here is the call graph for this function:

◆ getIdentity() [2/2]

tarch::la::Vector< ThreePowerD, std::complex< double > > toolbox::finiteelements::getIdentity ( const tarch::la::Vector< Dimensions, std::complex< double > > & h)

Definition at line 750 of file StencilFactory.cpp.

References assertionMsg, get1DIdentity(), and stencilProduct().

Here is the call graph for this function:

◆ getLaplacian() [1/3]

tarch::la::Vector< ThreePowerD, double > toolbox::finiteelements::getLaplacian ( const tarch::la::Vector< Dimensions, double > & h,
const tarch::la::Vector< Dimensions, double > & scaling = 1.0 )

Computes the Laplacian.

The individual derivatives are scaled with the entries of scaling. For the standard Laplacian, you can use the scaling 1.0. For a scalar material parameter, either pass 1.0 and rescale the result or pass in a vector where all entries hold the material parameter.

Definition at line 297 of file StencilFactory.cpp.

References assertionMsg, get1DLaplaceStencil(), get1DMassStencil(), and stencilProduct().

Referenced by getPoissonMatrixWithJumpingCoefficient(), and hierarchicalTransform().

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

◆ getLaplacian() [2/3]

tarch::la::Vector< ThreePowerD, std::complex< double > > toolbox::finiteelements::getLaplacian ( const tarch::la::Vector< Dimensions, double > & h,
const tarch::la::Vector< Dimensions, std::complex< double > > & scaling )

Definition at line 454 of file StencilFactory.cpp.

References assertionMsg, get1DLaplaceStencil(), get1DMassStencil(), and stencilProduct().

Here is the call graph for this function:

◆ getLaplacian() [3/3]

tarch::la::Vector< ThreePowerD, std::complex< double > > toolbox::finiteelements::getLaplacian ( const tarch::la::Vector< Dimensions, std::complex< double > > & h,
const tarch::la::Vector< Dimensions, std::complex< double > > & scaling = std::complex<double>(1.0) )

Definition at line 341 of file StencilFactory.cpp.

References assertionMsg, get1DLaplaceStencil(), get1DMassStencil(), and stencilProduct().

Here is the call graph for this function:

◆ getMassMatrix() [1/2]

tarch::la::Vector< ThreePowerD, double > toolbox::finiteelements::getMassMatrix ( const tarch::la::Vector< Dimensions, double > & h)

Computes the mass matrix.

The mass results results from an identity in a PDE that is interpreted in the finite element sense. As such, it is scaled with h^d. And it is not a nodal evaluation but a nine-point stencil (in 2d).

Definition at line 618 of file StencilFactory.cpp.

References assertionMsg, get1DMassStencil(), and stencilProduct().

Here is the call graph for this function:

◆ getMassMatrix() [2/2]

tarch::la::Vector< ThreePowerD, std::complex< double > > toolbox::finiteelements::getMassMatrix ( const tarch::la::Vector< Dimensions, std::complex< double > > & h)

Definition at line 640 of file StencilFactory.cpp.

References assertionMsg, get1DMassStencil(), and stencilProduct().

Here is the call graph for this function:

◆ getPoissonMatrixWithJumpingCoefficient()

toolbox::finiteelements::ElementWiseAssemblyMatrix toolbox::finiteelements::getPoissonMatrixWithJumpingCoefficient ( const tarch::la::Vector< Dimensions, double > & cellCentre,
const tarch::la::Vector< Dimensions, double > & h,
int integrationPointsPerAxis,
std::function< double(const tarch::la::Vector< Dimensions, double > &)> epsilon )

Integrate over element with a subvoxel scheme (cmp.

Murray) and approximate epsilon within each of the voxels as a constant. This means we basically use the reference matrix per subvoxel, multiply it with \( \left( - \frac{i^d}{h^d} + \frac{(i+1)^d}{h^d} \right)\left( - \frac{j^d}{h^d} + \frac{(j+1)^d}{h^d} \right)\dots \) and then add it to the rows. The i and j here are the logical Cartesian coordinates of the subvoxel which have to be mirrored accordingly if we map the reference matrix onto the respective vertices.

Definition at line 848 of file StencilFactory.cpp.

References dfor, dfor2, enddforx, getElementWiseAssemblyMatrix(), getLaplacian(), tarch::la::multiplyComponents(), and tarch::la::volume().

Referenced by examples::regulargridupscaling::MyObserver::enterCell().

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

◆ getStencilEntryLinearisedIndex()

int toolbox::finiteelements::getStencilEntryLinearisedIndex ( const tarch::la::Vector< Dimensions, int > stencilEntry)

Definition at line 4 of file Stencil.cpp.

Referenced by mapElementMatrixEntryOntoStencilEntry().

Here is the caller graph for this function:

◆ getUnitTests()

tarch::tests::TestCase * toolbox::finiteelements::getUnitTests ( )

Please destroy after usage.

Definition at line 11 of file UnitTests.cpp.

Referenced by main(), and runTests().

Here is the caller graph for this function:

◆ getUpwindDiscretisedConvection()

tarch::la::Vector< ThreePowerD, double > toolbox::finiteelements::getUpwindDiscretisedConvection ( const tarch::la::Vector< Dimensions, double > & h,
const tarch::la::Vector< Dimensions, double > & convCoeff )

Computes a convection operator, defined by the convection coefficients, discretised with the first order upwind scheme.

The scheme is stabilisied through a Laplacian, i.e. it is really an upwind scheme and not central differences. However, the operation returns a Finite Difference scheme scaled with the Finite Element mesh widths. For most applications, this should do, though formally the stencil should be derived from piecewise linear shape functions.

Definition at line 498 of file StencilFactory.cpp.

References assertionMsg, and assertionNumericalEquals.

◆ hierarchicalTransform()

toolbox::finiteelements::ElementWiseAssemblyMatrix toolbox::finiteelements::hierarchicalTransform ( const ElementWiseAssemblyMatrix & matrix,
const tarch::la::Vector< Dimensions, double > & h,
double weight = 1.0 )
Parameters
weightWeight of the Laplacian. Use it to inject material parameter samples, e.g.

Definition at line 302 of file ElementMatrix.cpp.

References getElementWiseAssemblyMatrix(), getLaplacian(), and TwoPowerD.

Referenced by inverseHierarchicalTransform().

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

◆ inverseHierarchicalTransform()

toolbox::finiteelements::ElementWiseAssemblyMatrix toolbox::finiteelements::inverseHierarchicalTransform ( const ElementWiseAssemblyMatrix & matrix,
const tarch::la::Vector< Dimensions, double > & h,
double weight = 1.0 )

Definition at line 322 of file ElementMatrix.cpp.

References hierarchicalTransform().

Here is the call graph for this function:

◆ mapElementMatrixEntryOntoStencilEntry() [1/2]

int toolbox::finiteelements::mapElementMatrixEntryOntoStencilEntry ( const tarch::la::Vector< Dimensions, int > & row,
const tarch::la::Vector< Dimensions, int > & col )

Definition at line 198 of file ElementMatrix.cpp.

References getStencilEntryLinearisedIndex().

Here is the call graph for this function:

◆ mapElementMatrixEntryOntoStencilEntry() [2/2]

int toolbox::finiteelements::mapElementMatrixEntryOntoStencilEntry ( int row,
int col )

An element-wise assembly matrix is distilled from a sequence of \( 2^d \) stencils which are in turn \( 3^d \) double arrays.

This operation computes the inverse index mapping. You hand in a row of the matrix plus a column. The operation tells you which entry it is in the original stencil of the rowth vertex.

For the 2d variant, I've just hard-coded all entries. This might also be an option for 3d if it turns out that the operation is too slow. For all other dimensions, we realise a generic version of the lookup: We convert both parameters row and col into integer vectors that have either a zero or a one as argument. We then call the overloaded variant of the present function.

Definition at line 220 of file ElementMatrix.cpp.

References assertion, assertionEquals2, peano4::utils::dDelinearised(), mapElementMatrixEntryOntoStencilEntry(), and TwoPowerD.

Referenced by mapElementMatrixEntryOntoStencilEntry(), reconstructStencilFragments(), and reconstructUniformStencilFragments().

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

◆ preprocessBoundaryStencil()

void toolbox::finiteelements::preprocessBoundaryStencil ( Stencil & stencil,
const std::bitset< Dimensions *2 > & boundaryFaceNormals )

Preprocess boundary stencils held by boundary vertices.

If codes use getElementWiseAssemblyMatrix() within cells, this operation decomposes stencils assuming that all \( 2^d \) adjacent cells are visited. For interior degrees of freedom this works fine. See the top row in the cartoon below.

It does not hold if we are near to the boundary. See second row left. The six-point stencil corresponds to a homogeneous Neumann condition and shows the stencil that has to be evaluated if the boundary holds a real degree of freedom.

If we use getElementWiseAssemblyMatrix(), this stencil is decomposed according to the middle right illustration. It yields the wrong result as there are only cells visited. We could provide specialised versions of getElementWiseAssemblyMatrix() that are evaluated at the boundary. However, it is often more convenient to treat all cells exactly the same way but to modify the stencil entries (red, bottom left) such that getElementWiseAssemblyMatrix() yields the right evaluation in the end.

This is exactly the job of this preprocessing routine that accepts the 'real' stencil (top left) plus a bitfield describing whether a boundary runs vertically through the degree domain, e.g. (first entry in boundaryFaceNormals is set as the normal is parallel to the first axis). The routine then upscales the stencil entries such that an evaluation in the end yields the right result.

Parameters
boundaryFaceNormalsWhere is the boundary. The first Dimensions entry refer to the faces running through the left bottom vertex.

Definition at line 8 of file StencilFactory.cpp.

References assertion2, dfore, and peano4::utils::dLinearised().

Here is the call graph for this function:

◆ reconstructStencilFragments()

tarch::la::Vector< TwoPowerD *ThreePowerD, double > toolbox::finiteelements::reconstructStencilFragments ( const ElementWiseAssemblyMatrix & matrix)

Reconstruct stencils from element matrices.

Generalised version of reconstructUniformStencilFragments(). It distinguishes the various vertices adjacent to a cell. However, it does not make any assumptions about stencil homogeneity and thus is only able to fill in some stencil entries partially.

Definition at line 84 of file ElementMatrix.cpp.

References assertionEquals1, dfor2, enddforx, mapElementMatrixEntryOntoStencilEntry(), and ThreePowerD.

Here is the call graph for this function:

◆ reconstructUniformStencilFragments()

toolbox::finiteelements::Stencil toolbox::finiteelements::reconstructUniformStencilFragments ( const ElementWiseAssemblyMatrix & matrix)

Reconstruct stencils from element matrices.

With getElementWiseAssemblyMatrix() you can construct the element-wise assembly matrix from the stencils. This operation is the inverse. You plug in the assembly matrix and you obtain the corresponding stencils. Hereby, the operation assumes that the stencils are the same for all \( 2^d \) vertices.

Without such an assumption, we were not able to reconstruct the stencils really. The element-wise matrix defines only some stencil entries, whereas the stencils also span neighbouring cells. Only because of the assumption on the stencil uniformity, we can construct all stencil entries. As a result, this operation does not work properly if you have changing stencils (varying coefficients, e.g.).

Definition at line 70 of file ElementMatrix.cpp.

References dfor2, enddforx, and mapElementMatrixEntryOntoStencilEntry().

Here is the call graph for this function:

◆ stencilProduct() [1/10]

tarch::la::Vector< 3 *3, double > toolbox::finiteelements::stencilProduct ( const tarch::la::Vector< 3, double > & a,
const tarch::la::Vector< 3, double > & b )

Makes stencil-product of two stencils:

a * b = [a_1, a_2, a_3] o [ b_1, b_2, B_3] [a_1*b_3 a_2*b_3 a_3*b_3] = [a_1*b_2 a_2*b_2 a_3*b_2] [a_1*b_1 a_2*b_1 a_3*b_1]

Definition at line 133 of file StencilFactory.cpp.

References a.

◆ stencilProduct() [2/10]

tarch::la::Vector< 3 *3 *3, double > toolbox::finiteelements::stencilProduct ( const tarch::la::Vector< 3, double > & a,
const tarch::la::Vector< 3, double > & b,
const tarch::la::Vector< 3, double > & c )

Equals a * (b*c)

Definition at line 165 of file StencilFactory.cpp.

References a.

◆ stencilProduct() [3/10]

tarch::la::Vector< 3 *3 *3 *3, double > toolbox::finiteelements::stencilProduct ( const tarch::la::Vector< 3, double > & a,
const tarch::la::Vector< 3, double > & b,
const tarch::la::Vector< 3, double > & c,
const tarch::la::Vector< 3, double > & d )

Definition at line 203 of file StencilFactory.cpp.

References a.

◆ stencilProduct() [4/10]

tarch::la::Vector< 3 *3 *3 *3 *3, double > toolbox::finiteelements::stencilProduct ( const tarch::la::Vector< 3, double > & a,
const tarch::la::Vector< 3, double > & b,
const tarch::la::Vector< 3, double > & c,
const tarch::la::Vector< 3, double > & d,
const tarch::la::Vector< 3, double > & e )

Definition at line 247 of file StencilFactory.cpp.

References a.

◆ stencilProduct() [5/10]

tarch::la::Vector< 5 *5, double > toolbox::finiteelements::stencilProduct ( const tarch::la::Vector< 5, double > & a,
const tarch::la::Vector< 5, double > & b )

Definition at line 149 of file StencilFactory.cpp.

References a.

◆ stencilProduct() [6/10]

tarch::la::Vector< 5 *5 *5, double > toolbox::finiteelements::stencilProduct ( const tarch::la::Vector< 5, double > & a,
const tarch::la::Vector< 5, double > & b,
const tarch::la::Vector< 5, double > & c )

Definition at line 184 of file StencilFactory.cpp.

References a.

◆ stencilProduct() [7/10]

tarch::la::Vector< 5 *5 *5 *5, double > toolbox::finiteelements::stencilProduct ( const tarch::la::Vector< 5, double > & a,
const tarch::la::Vector< 5, double > & b,
const tarch::la::Vector< 5, double > & c,
const tarch::la::Vector< 5, double > & d )

Definition at line 225 of file StencilFactory.cpp.

References a.

◆ stencilProduct() [8/10]

tarch::la::Vector< 5 *5 *5 *5 *5, double > toolbox::finiteelements::stencilProduct ( const tarch::la::Vector< 5, double > & a,
const tarch::la::Vector< 5, double > & b,
const tarch::la::Vector< 5, double > & c,
const tarch::la::Vector< 5, double > & d,
const tarch::la::Vector< 5, double > & e )

Definition at line 272 of file StencilFactory.cpp.

References a.

◆ stencilProduct() [9/10]

◆ stencilProduct() [10/10]

template<int StencilSize>
tarch::la::Vector< StencilSize *StencilSize *StencilSize, double > toolbox::finiteelements::stencilProduct ( const tarch::la::Vector< StencilSize, double > & a,
const tarch::la::Vector< StencilSize, double > & b,
const tarch::la::Vector< StencilSize, double > & c )

Definition at line 20 of file ElementMatrix.cpph.

References a.