5#ifndef DUNE_LOCALFUNCTIONS_RAVIARTTHOMAS_RAVIARTTHOMASSIMPLEX_RAVIARTTHOMASSIMPLEXINTERPOLATION_HH
6#define DUNE_LOCALFUNCTIONS_RAVIARTTHOMAS_RAVIARTTHOMASSIMPLEX_RAVIARTTHOMASSIMPLEXINTERPOLATION_HH
14#include <dune/geometry/referenceelements.hh>
18#include <dune/localfunctions/common/localkey.hh>
19#include <dune/localfunctions/utility/interpolationhelper.hh>
20#include <dune/localfunctions/utility/polynomialbasis.hh>
21#include <dune/localfunctions/orthonormal/orthonormalbasis.hh>
29 template <
unsigned int dim,
class Field >
30 struct RaviartThomasL2InterpolationFactory;
37 class LocalCoefficientsContainer
39 typedef LocalCoefficientsContainer This;
42 template <
class Setter>
43 LocalCoefficientsContainer (
const Setter &setter )
45 setter.setLocalKeys(localKey_);
48 const LocalKey &localKey (
const unsigned int i )
const
51 return localKey_[ i ];
54 std::size_t
size ()
const
56 return localKey_.size();
60 std::vector< LocalKey > localKey_;
68 template <
unsigned int dim >
69 struct RaviartThomasCoefficientsFactory
71 typedef std::size_t Key;
72 typedef const LocalCoefficientsContainer Object;
74 template< GeometryType::Id geometryId >
75 static Object *create(
const Key &key )
77 typedef RaviartThomasL2InterpolationFactory< dim, double > InterpolationFactory;
78 if( !supports< geometryId >( key ) )
80 typename InterpolationFactory::Object *interpolation = InterpolationFactory::template create< geometryId >( key );
81 Object *localKeys =
new Object( *interpolation );
82 InterpolationFactory::release( interpolation );
86 template< GeometryType::Id geometryId >
87 static bool supports (
const Key &key )
91 static void release( Object *
object ) {
delete object; }
105 template<
unsigned int dim,
class Field >
106 struct RTL2InterpolationBuilder
108 static const unsigned int dimension = dim;
111 typedef OrthonormalBasisFactory< dimension, Field > TestBasisFactory;
112 typedef typename TestBasisFactory::Object TestBasis;
115 typedef OrthonormalBasisFactory< dimension-1, Field > TestFaceBasisFactory;
116 typedef typename TestFaceBasisFactory::Object TestFaceBasis;
121 RTL2InterpolationBuilder () =
default;
123 RTL2InterpolationBuilder (
const RTL2InterpolationBuilder & ) =
delete;
124 RTL2InterpolationBuilder ( RTL2InterpolationBuilder && ) =
delete;
126 ~RTL2InterpolationBuilder ()
128 TestBasisFactory::release( testBasis_ );
129 for( FaceStructure &f : faceStructure_ )
130 TestFaceBasisFactory::release( f.basis_ );
135 std::size_t order ()
const {
return order_; }
138 unsigned int faceSize ()
const {
return faceSize_; }
141 TestBasis *testBasis ()
const {
return testBasis_; }
144 TestFaceBasis *testFaceBasis (
unsigned int f )
const { assert( f < faceSize() );
return faceStructure_[ f ].basis_; }
147 const Normal normal (
unsigned int f )
const { assert( f < faceSize() );
return faceStructure_[ f ].normal_; }
149 template< GeometryType::Id geometryId >
150 void build ( std::size_t order )
153 geometry_ = geometry;
156 testBasis_ = (order > 0 ? TestBasisFactory::template create< geometry >( order-1 ) : nullptr);
159 faceSize_ = refElement.size( 1 );
160 faceStructure_.reserve( faceSize_ );
161 for(
unsigned int face = 0; face < faceSize_; ++face )
174 TestFaceBasis *faceBasis = Impl::toGeometryTypeIdConstant<dimension-1>(refElement.type( face, 1 ), [&](
auto faceGeometryTypeId) {
175 return TestFaceBasisFactory::template create< decltype(faceGeometryTypeId)::value >( order );
177 faceStructure_.emplace_back( faceBasis, refElement.integrationOuterNormal( face ) );
179 assert( faceStructure_.size() == faceSize_ );
185 FaceStructure( TestFaceBasis *tfb,
const Normal n )
186 : basis_( tfb ), normal_( n )
189 TestFaceBasis *basis_;
193 std::vector< FaceStructure > faceStructure_;
194 TestBasis *testBasis_ =
nullptr;
196 unsigned int faceSize_;
210 template<
unsigned int dimension,
class F>
212 :
public InterpolationHelper< F ,dimension >
215 typedef InterpolationHelper<F,dimension> Base;
219 typedef RTL2InterpolationBuilder<dimension,Field> Builder;
225 template<
class Function,
class Vector,
226 decltype(std::declval<Vector>().size(),
bool{}) =
true,
227 decltype(std::declval<Vector>().resize(0u),
bool{}) =
true>
228 void interpolate (
const Function &function, Vector &coefficients )
const
230 coefficients.resize(size());
231 typename Base::template Helper<Function,Vector,true> func( function,coefficients );
235 template<
class Basis,
class Matrix,
236 decltype(std::declval<Matrix>().rows(),
bool{}) =
true,
237 decltype(std::declval<Matrix>().cols(),
bool{}) =
true,
238 decltype(std::declval<Matrix>().resize(0u,0u),
bool{}) =
true>
239 void interpolate (
const Basis &basis,
Matrix &matrix )
const
241 matrix.resize( size(), basis.size() );
242 typename Base::template Helper<Basis,Matrix,false> func( basis,matrix );
246 std::size_t order()
const
250 std::size_t size()
const
254 template <GeometryType::Id geometryId>
255 void build( std::size_t order )
259 builder_.template build<geometryId>(order_);
260 if (builder_.testBasis())
261 size_ += dimension*builder_.testBasis()->size();
262 for (
unsigned int f=0; f<builder_.faceSize(); ++f )
263 if (builder_.testFaceBasis(f))
264 size_ += builder_.testFaceBasis(f)->size();
267 void setLocalKeys(std::vector< LocalKey > &keys)
const
270 unsigned int row = 0;
271 for (
unsigned int f=0; f<builder_.faceSize(); ++f)
273 if (builder_.faceSize())
274 for (
unsigned int i=0; i<builder_.testFaceBasis(f)->size(); ++i,++row)
277 if (builder_.testBasis())
278 for (
unsigned int i=0; i<builder_.testBasis()->size()*dimension; ++i,++row)
280 assert( row == size() );
284 template<
class Func,
class Container,
bool type >
285 void interpolate (
typename Base::template Helper<Func,Container,type> &func )
const
289 std::vector< Field > testBasisVal;
291 for (
unsigned int i=0; i<size(); ++i)
292 for (
unsigned int j=0; j<func.size(); ++j)
295 unsigned int row = 0;
303 for (
unsigned int f=0; f<builder_.faceSize(); ++f)
305 if (!builder_.testFaceBasis(f))
307 testBasisVal.resize(builder_.testFaceBasis(f)->size());
309 const auto &geometry = refElement.template geometry< 1 >( f );
311 const FaceQuadrature &faceQuad = FaceQuadratureRules::rule( subGeoType, 2*order_+2 );
313 const unsigned int quadratureSize = faceQuad.size();
314 for(
unsigned int qi = 0; qi < quadratureSize; ++qi )
317 builder_.testFaceBasis(f)->template evaluate<0>(faceQuad[qi].position(),testBasisVal);
319 testBasisVal[0] = 1.;
320 fillBnd( row, testBasisVal,
321 func.evaluate( geometry.global( faceQuad[qi].position() ) ),
322 builder_.normal(f), faceQuad[qi].weight(),
326 row += builder_.testFaceBasis(f)->size();
329 if (builder_.testBasis())
331 testBasisVal.resize(builder_.testBasis()->size());
337 const unsigned int quadratureSize = elemQuad.size();
338 for(
unsigned int qi = 0; qi < quadratureSize; ++qi )
340 builder_.testBasis()->template evaluate<0>(elemQuad[qi].position(),testBasisVal);
341 fillInterior( row, testBasisVal,
342 func.evaluate(elemQuad[qi].position()),
343 elemQuad[qi].weight(),
347 row += builder_.testBasis()->size()*dimension;
362 template <
class MVal,
class RTVal,
class Matrix>
363 void fillBnd (
unsigned int startRow,
370 const unsigned int endRow = startRow+mVal.size();
371 typename RTVal::const_iterator rtiter = rtVal.begin();
372 for (
unsigned int col = 0; col < rtVal.size() ; ++rtiter,++col)
374 Field cFactor = (*rtiter)*normal;
375 typename MVal::const_iterator miter = mVal.
begin();
376 for (
unsigned int row = startRow;
377 row!=endRow; ++miter, ++row )
379 matrix.add(row,col, (weight*cFactor)*(*miter) );
381 assert( miter == mVal.end() );
392 template <
class MVal,
class RTVal,
class Matrix>
393 void fillInterior (
unsigned int startRow,
399 const unsigned int endRow = startRow+mVal.size()*dimension;
400 typename RTVal::const_iterator rtiter = rtVal.begin();
401 for (
unsigned int col = 0; col < rtVal.size() ; ++rtiter,++col)
403 typename MVal::const_iterator miter = mVal.begin();
404 for (
unsigned int row = startRow;
405 row!=endRow; ++miter,row+=dimension )
407 for (
unsigned int i=0; i<dimension; ++i)
409 matrix.add(row+i,col, (weight*(*miter))*(*rtiter)[i] );
412 assert( miter == mVal.end() );
421 template <
unsigned int dim,
class Field >
422 struct RaviartThomasL2InterpolationFactory
424 typedef RTL2InterpolationBuilder<dim,Field> Builder;
426 typedef std::size_t Key;
427 typedef typename std::remove_const<Object>::type NonConstObject;
429 template <GeometryType::Id geometryId>
430 static Object *create(
const Key &key )
432 if ( !supports<geometryId>(key) )
434 NonConstObject *interpol =
new NonConstObject();
435 interpol->template build<geometryId>(key);
438 template< GeometryType::Id geometryId >
439 static bool supports (
const Key &key )
441 return GeometryType(geometryId).isSimplex();
443 static void release( Object *
object ) {
delete object; }
Iterator begin()
begin iterator
Definition: densevector.hh:347
Unique label for each type of entities that can occur in DUNE grids.
Definition: type.hh:114
Describe position of one degree of freedom.
Definition: localkey.hh:24
A generic dynamic dense matrix.
Definition: matrix.hh:561
Abstract base class for quadrature rules.
Definition: quadraturerules.hh:214
A container for all quadrature rules of dimension dim
Definition: quadraturerules.hh:260
static const QuadratureRule & rule(const GeometryType &t, int p, QuadratureType::Enum qt=QuadratureType::GaussLegendre)
select the appropriate QuadratureRule for GeometryType t and order p
Definition: quadraturerules.hh:326
An L2-based interpolation for Raviart Thomas.
Definition: raviartthomassimplexinterpolation.hh:213
A few common exception classes.
GeometryType
Type representing VTK's entity geometry types.
Definition: common.hh:132
Dune namespace.
Definition: alignedallocator.hh:13
constexpr std::integral_constant< std::size_t, sizeof...(II)> size(std::integer_sequence< T, II... >)
Return the size of the sequence.
Definition: integersequence.hh:75
static const ReferenceElement & general(const GeometryType &type)
get general reference elements
Definition: referenceelements.hh:156
A unique label for each type of element that can occur in a grid.
Helper classes to provide indices for geometrytypes for use in a vector.