Commit b0e6a264 authored by Henrik Stolzmann's avatar Henrik Stolzmann
Browse files

Add comments and explanations to the construction algorithm for high-order RaviartThomas functions.

Furthermore removed the unused member variable `row1_`.
parent fdc3c686
......@@ -12,6 +12,24 @@
namespace Dune
{
/*
* `RTPreBasisFactory` provides a basis for the Raviart-Thomas function space.
* `RaviartThomasL2InterpolationFactory` provides the linear functionals.
*
* `Defaultbasisfactory::create` first builds the function space and the linear functionals.
* Then the constructor of `BasisMatrix` gets called. There the matrix
*
* \begin{equation}
* A_{i,j} := N_j(\phi_i)
* \end{equation}
*
* with linear functionals $N_j$ and basisfunctions $\phi_i$ gets assembled.
* Then the matrix gets inverted and is then used as a coefficent matrix for the standard monomial basis.
*
* For more details on the theory see the first chapter "Construction of Local Finite Element Spaces Using the Generic Reference Elements"
* of the book "Advances in Dune" by Dedner, Flemisch and Klöfkorn published in 2012.
*/
template< unsigned int dim, class SF, class CF >
struct RaviartThomasBasisFactory
: public DefaultBasisFactory< RTPreBasisFactory<dim,CF>,
......
......@@ -103,10 +103,16 @@ namespace Dune
struct RTL2InterpolationBuilder
{
static const unsigned int dimension = dim;
// for the dofs associated to the element
typedef OrthonormalBasisFactory< dimension, Field > TestBasisFactory;
typedef typename TestBasisFactory::Object TestBasis;
// for the dofs associated to the faces
typedef OrthonormalBasisFactory< dimension-1, Field > TestFaceBasisFactory;
typedef typename TestFaceBasisFactory::Object TestFaceBasis;
// the normals of the faces
typedef FieldVector< Field, dimension > Normal;
RTL2InterpolationBuilder () = default;
......@@ -127,11 +133,16 @@ namespace Dune
std::size_t order () const { return order_; }
// number of faces
unsigned int faceSize () const { return faceSize_; }
// basis associated to the element
TestBasis *testBasis () const { return testBasis_; }
// basis associated to face f
TestFaceBasis *testFaceBasis ( unsigned int f ) const { assert( f < faceSize() ); return faceStructure_[ f ].basis_; }
// normal of face f
const Normal &normal ( unsigned int f ) const { assert( f < faceSize() ); return *(faceStructure_[ f ].normal_); }
template< class Topology >
......@@ -147,6 +158,17 @@ namespace Dune
faceStructure_.reserve( faceSize_ );
for( unsigned int face = 0; face < faceSize_; ++face )
{
/* For simplices or cubes of arbitrary dimension you could just use
*
* ```
* GeometryType faceGeometry = Impl::getBase(geometry_);
* TestFaceBasis *faceBasis = TestFaceBasisFactory::template create< faceGeometry >( order );
* ```
*
* For i.e. Prisms and Pyramids in 3d this does not work because they contain squares and triangles as faces.
* And depending on the dynamic face index a different face geometry is needed.
*
*/
TestFaceBasis *faceBasis = Impl::IfTopology< CreateFaceBasis, dimension-1 >::apply( refElement.type( face, 1 ).id(), order );
faceStructure_.emplace_back( faceBasis, refElement.integrationOuterNormal( face ) );
}
......@@ -327,7 +349,15 @@ namespace Dune
}
private:
/** /brief evaluate boundary functionals **/
/** \brief evaluate boundary functionals
*
* \param startRow row of matrix to start
* \param mVal value of the testBasis at a quadrature point on a face
* \param rtVal value of the RaviartThomasBasis at a quadrature point on a face
* \param normal the normal of the face
* \param weight quadrature weight
* \param matrix result gets written into matrix starting with row: row
*/
template <class MVal, class RTVal,class Matrix>
void fillBnd (unsigned int startRow,
const MVal &mVal,
......@@ -350,6 +380,14 @@ namespace Dune
assert( miter == mVal.end() );
}
}
/** \brief evaluate interior functionals
*
* \param startRow row of matrix to start
* \param mVal value of the testBasis at a quadrature point in the interior of the ReferenceElement
* \param rtVal value of the RaviartThomasBasis at a quadrature point in the interior of the ReferenceElement
* \param weight quadrature weight
* \param matrix result gets written into matrix starting with row: row
*/
template <class MVal, class RTVal,class Matrix>
void fillInterior (unsigned int startRow,
const MVal &mVal,
......
......@@ -36,13 +36,13 @@ namespace Dune
{
RTVecMatrix<Topology,Field> vecMatrix(order);
MBasis *mbasis = MBasisFactory::template create<Topology>(order+1);
typename std::remove_const<Object>::type *tmBasis =
new typename std::remove_const<Object>::type(*mbasis);
typename std::remove_const<Object>::type *tmBasis = new typename std::remove_const<Object>::type(*mbasis);
tmBasis->fill(vecMatrix);
return tmBasis;
}
static void release( Object *object ) { delete object; }
};
template <class Topology, class Field>
struct RTVecMatrix
{
......@@ -51,28 +51,84 @@ namespace Dune
typedef MonomialBasis<Topology,MI> MIBasis;
RTVecMatrix(std::size_t order)
{
/*
* Construction of Raviart-Thomas elements in high dimensions see "Mixed Finite Elements in \R^3" by Nedelec, 1980.
*
* Let $\P_{n,k}$ be the space of polynomials in $n$ variables with degree $\leq k$.
* The space of Raviart-Thomas functions in $n$ dimensions with index $k$ is defined as
*
* \begin{equation*}
* RT_k := (\P_{k-1})^n \oplus \widetilde \P_k x
* \end{equation*}
* with $x=(x_1,x_2,\dots, x_n)$ in $n$ dimensions and $\widetilde \P_k$ the homogeneous polynomials of degree $k$.
*
* For $RT_k$ holds
* \begin{equation*}
* (\P_{k-1})^n \subset RT_k \subset (\P_k)^n.
* \end{equation*}
*
* We construct $(\P_k)^n$ and and only use the monomials contained in $RT_k$.
*
*/
MIBasis basis(order+1);
FieldVector< MI, dim > x;
/*
* Init MultiIndices
* x[0]=(1,0,0) x
* x[1]=(0,1,0) y
* x[2]=(0,0,1) z
*/
for( unsigned int i = 0; i < dim; ++i )
x[ i ].set( i, 1 );
std::vector< MI > val( basis.size() );
// val now contains all monomials in $n$ dimensions with degree $\leq order+1$
basis.evaluate( x, val );
col_ = basis.size();
// get $\dim (\P_{order-1})$
unsigned int notHomogen = 0;
if (order>0)
notHomogen = basis.sizes()[order-1];
// get $\dim \widetilde (\P_order)$
unsigned int homogen = basis.sizes()[order]-notHomogen;
/*
*
* The set $RT_k$ is defined as
*
* \begin{equation}
* RT_k := (\P_k)^dim + \widetilde \P_k x with x\in \R^n.
* \end{equation}
*
* The space $\P_k$ is split in $\P_k = \P_{k-1} + \widetilde \P_k$.
*
* \begin{align}
* RT_k &= (\P_{k-1} + \widetilde \P_k)^dim + \widetilde \P_k x with x\in \R^n
* &= (\P_{k-1})^n + (\widetilde \P_k)^n + \widetilde \P_k x with x\in \R^n
* \end{align}
*
* Thus $\dim RT_k = n * \dim \P_{k-1} + (n+1)*\dim (\widetilde \P_k)$
*/
// row_ = \dim RT_k *dim
row_ = (notHomogen*dim+homogen*(dim+1))*dim;
row1_ = basis.sizes()[order]*dim*dim;
mat_ = new Field*[row_];
int row = 0;
/* Assign the correct values for the nonhomogeneous polymonials $p\in (\P_{oder-1})^dim$
* A basis function is represented by $dim$ rows.
*/
for (unsigned int i=0; i<notHomogen+homogen; ++i)
{
for (unsigned int r=0; r<dim; ++r)
{
for (unsigned int rr=0; rr<dim; ++rr)
{
// init row to zero
mat_[row] = new Field[col_];
for (unsigned int j=0; j<col_; ++j)
{
......@@ -84,16 +140,23 @@ namespace Dune
}
}
}
/* Assign the correct values for the homogeneous polymonials $p\in RT_k \backslash (\P_{oder-1})^dim$
* A basis function is represented by $dim$ rows.
*/
for (unsigned int i=0; i<homogen; ++i)
{
for (unsigned int r=0; r<dim; ++r)
{
// init rows to zero
mat_[row] = new Field[col_];
for (unsigned int j=0; j<col_; ++j)
{
mat_[row][j] = 0.;
}
unsigned int w;
// get a monomial $ p \in \widetilde \P_{order}$
MI xval = val[notHomogen+i];
xval *= x[r];
for (w=homogen+notHomogen; w<val.size(); ++w)
......@@ -109,6 +172,7 @@ namespace Dune
}
}
}
~RTVecMatrix()
{
for (unsigned int i=0; i<rows(); ++i) {
......@@ -116,12 +180,15 @@ namespace Dune
}
delete [] mat_;
}
unsigned int cols() const {
return col_;
}
unsigned int rows() const {
return row_;
}
template <class Vector>
void row( const unsigned int row, Vector &vec ) const
{
......@@ -130,7 +197,7 @@ namespace Dune
for (unsigned int i=0; i<N; ++i)
field_cast(mat_[row][i],vec[i]);
}
unsigned int row_,col_,row1_;
unsigned int row_,col_;
Field **mat_;
};
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment