3#ifndef DUNE_LOCALFUNCTIONS_RAVIARTTHOMAS_RAVIARTTHOMASSIMPLEX_RAVIARTTHOMASSIMPLEXINTERPOLATION_HH
4#define DUNE_LOCALFUNCTIONS_RAVIARTTHOMAS_RAVIARTTHOMASSIMPLEX_RAVIARTTHOMASSIMPLEXINTERPOLATION_HH
12#include <dune/geometry/referenceelements.hh>
15#include <dune/localfunctions/common/localkey.hh>
16#include <dune/localfunctions/utility/interpolationhelper.hh>
17#include <dune/localfunctions/utility/polynomialbasis.hh>
18#include <dune/localfunctions/orthonormal/orthonormalbasis.hh>
26 template <
unsigned int dim,
class Field >
27 struct RaviartThomasL2InterpolationFactory;
34 class LocalCoefficientsContainer
36 typedef LocalCoefficientsContainer This;
39 template <
class Setter>
40 LocalCoefficientsContainer (
const Setter &setter )
42 setter.setLocalKeys(localKey_);
45 const LocalKey &localKey (
const unsigned int i )
const
48 return localKey_[ i ];
51 std::size_t size ()
const
53 return localKey_.size();
57 std::vector< LocalKey > localKey_;
65 template <
unsigned int dim >
66 struct RaviartThomasCoefficientsFactory
68 typedef std::size_t Key;
69 typedef const LocalCoefficientsContainer Object;
71 template<
class Topology >
72 static Object *create(
const Key &key )
74 typedef RaviartThomasL2InterpolationFactory< dim, double > InterpolationFactory;
75 if( !supports< Topology >( key ) )
77 typename InterpolationFactory::Object *interpolation = InterpolationFactory::template create< Topology >( key );
78 Object *localKeys =
new Object( *interpolation );
79 InterpolationFactory::release( interpolation );
83 template<
class Topology >
84 static bool supports (
const Key &key )
86 return Impl::IsSimplex< Topology >::value;
88 static void release( Object *
object ) {
delete object; }
102 template<
unsigned int dim,
class Field >
103 struct RTL2InterpolationBuilder
105 static const unsigned int dimension = dim;
106 typedef OrthonormalBasisFactory< dimension, Field > TestBasisFactory;
107 typedef typename TestBasisFactory::Object TestBasis;
108 typedef OrthonormalBasisFactory< dimension-1, Field > TestFaceBasisFactory;
109 typedef typename TestFaceBasisFactory::Object TestFaceBasis;
112 RTL2InterpolationBuilder () =
default;
114 RTL2InterpolationBuilder (
const RTL2InterpolationBuilder & ) =
delete;
115 RTL2InterpolationBuilder ( RTL2InterpolationBuilder && ) =
delete;
117 ~RTL2InterpolationBuilder ()
119 TestBasisFactory::release( testBasis_ );
120 for( FaceStructure &f : faceStructure_ )
121 TestFaceBasisFactory::release( f.basis_ );
124 unsigned int topologyId ()
const {
return topologyId_; }
128 std::size_t order ()
const {
return order_; }
130 unsigned int faceSize ()
const {
return faceSize_; }
132 TestBasis *testBasis ()
const {
return testBasis_; }
133 TestFaceBasis *testFaceBasis (
unsigned int f )
const { assert( f < faceSize() );
return faceStructure_[ f ].basis_; }
135 const Normal &normal (
unsigned int f )
const { assert( f < faceSize() );
return *(faceStructure_[ f ].normal_); }
137 template<
class Topology >
138 void build ( std::size_t order )
141 topologyId_ = Topology::id;
143 testBasis_ = (order > 0 ? TestBasisFactory::template create< Topology >( order-1 ) : nullptr);
146 faceSize_ = refElement.size( 1 );
147 faceStructure_.reserve( faceSize_ );
148 for(
unsigned int face = 0; face < faceSize_; ++face )
151 faceStructure_.emplace_back( faceBasis, refElement.integrationOuterNormal( face ) );
153 assert( faceStructure_.size() == faceSize_ );
159 FaceStructure( TestFaceBasis *tfb,
const Normal &n )
160 : basis_( tfb ), normal_( &n )
163 TestFaceBasis *basis_;
167 template<
class FaceTopology >
168 struct CreateFaceBasis
170 static TestFaceBasis *
apply ( std::size_t order ) {
return TestFaceBasisFactory::template create< FaceTopology >( order ); }
173 std::vector< FaceStructure > faceStructure_;
174 TestBasis *testBasis_ =
nullptr;
175 unsigned int topologyId_, faceSize_;
189 template<
unsigned int dimension,
class F>
191 :
public InterpolationHelper< F ,dimension >
194 typedef InterpolationHelper<F,dimension> Base;
198 typedef RTL2InterpolationBuilder<dimension,Field> Builder;
204 template<
class Function,
class Vector >
205 auto interpolate (
const Function &function, Vector &coefficients )
const
206 -> std::enable_if_t< std::is_same< decltype(std::declval<Vector>().resize(1) ),
void >::value,
void>
208 coefficients.resize(size());
209 typename Base::template Helper<Function,Vector,true> func( function,coefficients );
213 template<
class Basis,
class Matrix >
214 auto interpolate (
const Basis &basis,
Matrix &matrix )
const
215 -> std::enable_if_t< std::is_same<
216 decltype(std::declval<Matrix>().rowPtr(0)),
typename Matrix::Field* >::value,
void>
218 matrix.resize( size(), basis.size() );
219 typename Base::template Helper<Basis,Matrix,false> func( basis,matrix );
223 std::size_t order()
const
227 std::size_t size()
const
231 template <
class Topology>
232 void build( std::size_t order )
236 builder_.template build<Topology>(order_);
237 if (builder_.testBasis())
238 size_ += dimension*builder_.testBasis()->size();
239 for (
unsigned int f=0; f<builder_.faceSize(); ++f )
240 if (builder_.testFaceBasis(f))
241 size_ += builder_.testFaceBasis(f)->size();
244 void setLocalKeys(std::vector< LocalKey > &keys)
const
247 unsigned int row = 0;
248 for (
unsigned int f=0; f<builder_.faceSize(); ++f)
250 if (builder_.faceSize())
251 for (
unsigned int i=0; i<builder_.testFaceBasis(f)->size(); ++i,++row)
254 if (builder_.testBasis())
255 for (
unsigned int i=0; i<builder_.testBasis()->size()*dimension; ++i,++row)
257 assert( row == size() );
261 template<
class Func,
class Container,
bool type >
262 void interpolate (
typename Base::template Helper<Func,Container,type> &func )
const
266 std::vector< Field > testBasisVal;
268 for (
unsigned int i=0; i<size(); ++i)
269 for (
unsigned int j=0; j<func.size(); ++j)
272 unsigned int row = 0;
280 for (
unsigned int f=0; f<builder_.faceSize(); ++f)
282 if (!builder_.testFaceBasis(f))
284 testBasisVal.resize(builder_.testFaceBasis(f)->size());
286 const auto &geometry = refElement.template geometry< 1 >( f );
288 const FaceQuadrature &faceQuad = FaceQuadratureRules::rule( subGeoType, 2*order_+2 );
290 const unsigned int quadratureSize = faceQuad.size();
291 for(
unsigned int qi = 0; qi < quadratureSize; ++qi )
294 builder_.testFaceBasis(f)->template evaluate<0>(faceQuad[qi].position(),testBasisVal);
296 testBasisVal[0] = 1.;
297 fillBnd( row, testBasisVal,
298 func.evaluate( geometry.global( faceQuad[qi].position() ) ),
299 builder_.normal(f), faceQuad[qi].weight(),
303 row += builder_.testFaceBasis(f)->size();
306 if (builder_.testBasis())
308 testBasisVal.resize(builder_.testBasis()->size());
314 const unsigned int quadratureSize = elemQuad.size();
315 for(
unsigned int qi = 0; qi < quadratureSize; ++qi )
317 builder_.testBasis()->template evaluate<0>(elemQuad[qi].position(),testBasisVal);
318 fillInterior( row, testBasisVal,
319 func.evaluate(elemQuad[qi].position()),
320 elemQuad[qi].weight(),
324 row += builder_.testBasis()->size()*dimension;
331 template <
class MVal,
class RTVal,
class Matrix>
332 void fillBnd (
unsigned int startRow,
339 const unsigned int endRow = startRow+mVal.size();
340 typename RTVal::const_iterator rtiter = rtVal.begin();
341 for (
unsigned int col = 0; col < rtVal.size() ; ++rtiter,++col)
343 Field cFactor = (*rtiter)*normal;
344 typename MVal::const_iterator miter = mVal.
begin();
345 for (
unsigned int row = startRow;
346 row!=endRow; ++miter, ++row )
348 matrix.add(row,col, (weight*cFactor)*(*miter) );
350 assert( miter == mVal.end() );
353 template <
class MVal,
class RTVal,
class Matrix>
354 void fillInterior (
unsigned int startRow,
360 const unsigned int endRow = startRow+mVal.size()*dimension;
361 typename RTVal::const_iterator rtiter = rtVal.begin();
362 for (
unsigned int col = 0; col < rtVal.size() ; ++rtiter,++col)
364 typename MVal::const_iterator miter = mVal.begin();
365 for (
unsigned int row = startRow;
366 row!=endRow; ++miter,row+=dimension )
368 for (
unsigned int i=0; i<dimension; ++i)
370 matrix.add(row+i,col, (weight*(*miter))*(*rtiter)[i] );
373 assert( miter == mVal.end() );
382 template <
unsigned int dim,
class Field >
383 struct RaviartThomasL2InterpolationFactory
385 typedef RTL2InterpolationBuilder<dim,Field> Builder;
387 typedef std::size_t Key;
388 typedef typename std::remove_const<Object>::type NonConstObject;
389 template <
class Topology>
390 static Object *create(
const Key &key )
392 if ( !supports<Topology>(key) )
394 NonConstObject *interpol =
new NonConstObject();
395 interpol->template build<Topology>(key);
398 template<
class Topology >
399 static bool supports (
const Key &key )
401 return Impl::IsSimplex<Topology>::value;
403 static void release( Object *
object ) {
delete object; }
Iterator begin()
begin iterator
Definition: densevector.hh:348
Base class template for function classes.
Definition: function.hh:29
Unique label for each type of entities that can occur in DUNE grids.
Definition: type.hh:279
Describe position of one degree of freedom.
Definition: localkey.hh:21
A generic dynamic dense matrix.
Definition: matrix.hh:558
Abstract base class for quadrature rules.
Definition: quadraturerules.hh:126
A container for all quadrature rules of dimension dim
Definition: quadraturerules.hh:172
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:254
An L2-based interpolation for Raviart Thomas.
Definition: raviartthomassimplexinterpolation.hh:192
A few common exception classes.
GeometryType
Type representing VTK's entity geometry types.
Definition: common.hh:180
decltype(auto) apply(F &&f, ArgTuple &&args)
Apply function with arguments given as tuple.
Definition: apply.hh:46
Dune namespace.
Definition: alignedallocator.hh:14
static const ReferenceElement & general(const GeometryType &type)
get general reference elements
Definition: referenceelements.hh:196
A unique label for each type of element that can occur in a grid.