3 #ifndef DUNE_LOCALFUNCTIONS_RAVIARTTHOMAS_RAVIARTTHOMASSIMPLEX_RAVIARTTHOMASSIMPLEXINTERPOLATION_HH
4 #define DUNE_LOCALFUNCTIONS_RAVIARTTHOMAS_RAVIARTTHOMASSIMPLEX_RAVIARTTHOMASSIMPLEXINTERPOLATION_HH
9 #include <dune/common/exceptions.hh>
11 #include <dune/geometry/quadraturerules.hh>
12 #include <dune/geometry/referenceelements.hh>
13 #include <dune/geometry/type.hh>
14 #include <dune/geometry/typeindex.hh>
27 template <
unsigned int dim,
class Field >
28 struct RaviartThomasL2InterpolationFactory;
35 class LocalCoefficientsContainer
40 template <
class Setter>
43 setter.setLocalKeys(localKey_);
49 return localKey_[ i ];
54 return localKey_.size();
58 std::vector< LocalKey > localKey_;
66 template <
unsigned int dim >
69 typedef std::size_t
Key;
72 template< GeometryType::Id geometryId >
76 if( !supports< geometryId >( key ) )
78 typename InterpolationFactory::Object *interpolation = InterpolationFactory::template create< geometryId >( key );
80 InterpolationFactory::release( interpolation );
84 template< GeometryType::Id geometryId >
87 return GeometryType(geometryId).isSimplex();
103 template<
unsigned int dim,
class Field >
117 typedef FieldVector< Field, dimension >
Normal;
127 for( FaceStructure &f : faceStructure_ )
131 [[deprecated(
"Use type().id() instead.")]]
134 GeometryType
type ()
const {
return geometry_; }
136 std::size_t
order ()
const {
return order_; }
139 unsigned int faceSize ()
const {
return faceSize_; }
148 const Normal &
normal (
unsigned int f )
const { assert( f <
faceSize() );
return *(faceStructure_[ f ].normal_); }
150 template< GeometryType::Id geometryId >
153 constexpr GeometryType geometry = geometryId;
154 geometry_ = geometry;
157 testBasis_ = (
order > 0 ? TestBasisFactory::template create< geometry >(
order-1 ) :
nullptr);
159 const auto &refElement = ReferenceElements< Field, dimension >::general(
type() );
160 faceSize_ = refElement.size( 1 );
161 faceStructure_.reserve( faceSize_ );
162 for(
unsigned int face = 0; face < faceSize_; ++face )
175 TestFaceBasis *faceBasis = Impl::toGeometryTypeIdConstant<
dimension-1>(refElement.type( face, 1 ), [&](
auto faceGeometryTypeId) {
176 return TestFaceBasisFactory::template create< decltype(faceGeometryTypeId)::value >(
order );
178 faceStructure_.emplace_back( faceBasis, refElement.integrationOuterNormal( face ) );
180 assert( faceStructure_.size() == faceSize_ );
187 : basis_( tfb ), normal_( &n )
191 const Dune::FieldVector< Field, dimension > *normal_;
194 std::vector< FaceStructure > faceStructure_;
196 GeometryType geometry_;
197 unsigned int faceSize_;
211 template<
unsigned int dimension,
class F>
226 template<
class Function,
class Vector >
227 auto interpolate (
const Function &
function, Vector &coefficients )
const
228 -> std::enable_if_t< std::is_same< decltype(std::declval<Vector>().resize(1) ),
void >::
value,
void>
230 coefficients.resize(
size());
235 template<
class Basis,
class Matrix >
237 -> std::enable_if_t< std::is_same<
238 decltype(std::declval<Matrix>().rowPtr(0)),
typename Matrix::Field* >
::value,
void>
240 matrix.resize(
size(), basis.size() );
253 template <GeometryType::Id geometryId>
258 builder_.template build<geometryId>(order_);
261 for (
unsigned int f=0; f<builder_.
faceSize(); ++f )
269 unsigned int row = 0;
270 for (
unsigned int f=0; f<builder_.
faceSize(); ++f)
277 for (
unsigned int i=0; i<builder_.
testBasis()->
size()*dimension; ++i,++row)
279 assert( row ==
size() );
283 template<
class Func,
class Container,
bool type >
286 const Dune::GeometryType geoType = builder_.
type();
288 std::vector< Field > testBasisVal;
290 for (
unsigned int i=0; i<
size(); ++i)
291 for (
unsigned int j=0; j<func.size(); ++j)
294 unsigned int row = 0;
297 typedef Dune::QuadratureRule<
Field, dimension-1> FaceQuadrature;
298 typedef Dune::QuadratureRules<
Field, dimension-1> FaceQuadratureRules;
300 const auto &refElement = Dune::ReferenceElements< Field, dimension >::general( geoType );
302 for (
unsigned int f=0; f<builder_.
faceSize(); ++f)
308 const auto &geometry = refElement.template geometry< 1 >( f );
309 const Dune::GeometryType subGeoType( geometry.type().id(), dimension-1 );
310 const FaceQuadrature &faceQuad = FaceQuadratureRules::rule( subGeoType, 2*order_+2 );
312 const unsigned int quadratureSize = faceQuad.size();
313 for(
unsigned int qi = 0; qi < quadratureSize; ++qi )
316 builder_.
testFaceBasis(f)->template evaluate<0>(faceQuad[qi].position(),testBasisVal);
318 testBasisVal[0] = 1.;
319 fillBnd( row, testBasisVal,
320 func.evaluate( geometry.global( faceQuad[qi].position() ) ),
321 builder_.
normal(f), faceQuad[qi].weight(),
332 typedef Dune::QuadratureRule<Field, dimension> Quadrature;
333 typedef Dune::QuadratureRules<Field, dimension> QuadratureRules;
334 const Quadrature &elemQuad = QuadratureRules::rule( geoType, 2*order_+1 );
336 const unsigned int quadratureSize = elemQuad.size();
337 for(
unsigned int qi = 0; qi < quadratureSize; ++qi )
339 builder_.
testBasis()->template evaluate<0>(elemQuad[qi].position(),testBasisVal);
340 fillInterior( row, testBasisVal,
341 func.evaluate(elemQuad[qi].position()),
342 elemQuad[qi].weight(),
361 template <
class MVal,
class RTVal,
class Matrix>
362 void fillBnd (
unsigned int startRow,
365 const FieldVector<Field,dimension> &normal,
367 Matrix &matrix)
const
369 const unsigned int endRow = startRow+mVal.size();
370 typename RTVal::const_iterator rtiter = rtVal.begin();
371 for (
unsigned int col = 0; col < rtVal.size() ; ++rtiter,++col)
373 Field cFactor = (*rtiter)*normal;
374 typename MVal::const_iterator miter = mVal.begin();
375 for (
unsigned int row = startRow;
376 row!=endRow; ++miter, ++row )
378 matrix.add(row,col, (weight*cFactor)*(*miter) );
380 assert( miter == mVal.end() );
391 template <
class MVal,
class RTVal,
class Matrix>
392 void fillInterior (
unsigned int startRow,
396 Matrix &matrix)
const
398 const unsigned int endRow = startRow+mVal.size()*dimension;
399 typename RTVal::const_iterator rtiter = rtVal.begin();
400 for (
unsigned int col = 0; col < rtVal.size() ; ++rtiter,++col)
402 typename MVal::const_iterator miter = mVal.begin();
403 for (
unsigned int row = startRow;
404 row!=endRow; ++miter,row+=dimension )
406 for (
unsigned int i=0; i<dimension; ++i)
408 matrix.add(row+i,col, (weight*(*miter))*(*rtiter)[i] );
411 assert( miter == mVal.end() );
420 template <
unsigned int dim,
class Field >
428 template <GeometryType::Id geometryId>
431 if ( !supports<geometryId>(key) )
434 interpol->template build<geometryId>(key);
437 template< GeometryType::Id geometryId >
440 return GeometryType(geometryId).isSimplex();
Definition: bdfmcube.hh:16
@ value
Definition: tensor.hh:166
Describe position of one degree of freedom.
Definition: localkey.hh:21
Definition: nedelecsimplexinterpolation.hh:36
const LocalKey & localKey(const unsigned int i) const
Definition: raviartthomassimplexinterpolation.hh:46
LocalCoefficientsContainer(const Setter &setter)
Definition: nedelecsimplexinterpolation.hh:41
std::size_t size() const
Definition: nedelecsimplexinterpolation.hh:52
Definition: orthonormalbasis.hh:18
static void release(Object *object)
Definition: orthonormalbasis.hh:55
Definition: raviartthomassimplexinterpolation.hh:422
std::remove_const< Object >::type NonConstObject
Definition: raviartthomassimplexinterpolation.hh:426
static Object * create(const Key &key)
Definition: raviartthomassimplexinterpolation.hh:429
static void release(Object *object)
Definition: raviartthomassimplexinterpolation.hh:442
static bool supports(const Key &key)
Definition: raviartthomassimplexinterpolation.hh:438
RTL2InterpolationBuilder< dim, Field > Builder
Definition: raviartthomassimplexinterpolation.hh:423
const RaviartThomasL2Interpolation< dim, Field > Object
Definition: raviartthomassimplexinterpolation.hh:424
std::size_t Key
Definition: raviartthomassimplexinterpolation.hh:425
Definition: raviartthomassimplexinterpolation.hh:68
std::size_t Key
Definition: raviartthomassimplexinterpolation.hh:69
static Object * create(const Key &key)
Definition: raviartthomassimplexinterpolation.hh:73
static void release(Object *object)
Definition: raviartthomassimplexinterpolation.hh:89
static bool supports(const Key &key)
Definition: raviartthomassimplexinterpolation.hh:85
const LocalCoefficientsContainer Object
Definition: raviartthomassimplexinterpolation.hh:70
Definition: raviartthomassimplexinterpolation.hh:105
FieldVector< Field, dimension > Normal
Definition: raviartthomassimplexinterpolation.hh:117
TestBasis * testBasis() const
Definition: raviartthomassimplexinterpolation.hh:142
TestBasisFactory::Object TestBasis
Definition: raviartthomassimplexinterpolation.hh:110
TestFaceBasisFactory::Object TestFaceBasis
Definition: raviartthomassimplexinterpolation.hh:114
unsigned int faceSize() const
Definition: raviartthomassimplexinterpolation.hh:139
void build(std::size_t order)
Definition: raviartthomassimplexinterpolation.hh:151
RTL2InterpolationBuilder()=default
TestFaceBasis * testFaceBasis(unsigned int f) const
Definition: raviartthomassimplexinterpolation.hh:145
GeometryType type() const
Definition: raviartthomassimplexinterpolation.hh:134
RTL2InterpolationBuilder(const RTL2InterpolationBuilder &)=delete
OrthonormalBasisFactory< dimension-1, Field > TestFaceBasisFactory
Definition: raviartthomassimplexinterpolation.hh:113
RTL2InterpolationBuilder(RTL2InterpolationBuilder &&)=delete
OrthonormalBasisFactory< dimension, Field > TestBasisFactory
Definition: raviartthomassimplexinterpolation.hh:109
std::size_t order() const
Definition: raviartthomassimplexinterpolation.hh:136
static const unsigned int dimension
Definition: raviartthomassimplexinterpolation.hh:106
const Normal & normal(unsigned int f) const
Definition: raviartthomassimplexinterpolation.hh:148
unsigned int topologyId() const
Definition: raviartthomassimplexinterpolation.hh:132
~RTL2InterpolationBuilder()
Definition: raviartthomassimplexinterpolation.hh:124
An L2-based interpolation for Raviart Thomas.
Definition: raviartthomassimplexinterpolation.hh:214
std::size_t order() const
Definition: raviartthomassimplexinterpolation.hh:245
RaviartThomasL2Interpolation()
Definition: raviartthomassimplexinterpolation.hh:221
void interpolate(typename Base::template Helper< Func, Container, type > &func) const
Definition: raviartthomassimplexinterpolation.hh:284
auto interpolate(const Basis &basis, Matrix &matrix) const -> std::enable_if_t< std::is_same< decltype(std::declval< Matrix >().rowPtr(0)), typename Matrix::Field * >::value, void >
Definition: raviartthomassimplexinterpolation.hh:236
RTL2InterpolationBuilder< dimension, Field > Builder
Definition: raviartthomassimplexinterpolation.hh:220
F Field
Definition: raviartthomassimplexinterpolation.hh:219
void build(std::size_t order)
Definition: raviartthomassimplexinterpolation.hh:254
auto interpolate(const Function &function, Vector &coefficients) const -> std::enable_if_t< std::is_same< decltype(std::declval< Vector >().resize(1)), void >::value, void >
Definition: raviartthomassimplexinterpolation.hh:227
std::size_t size() const
Definition: raviartthomassimplexinterpolation.hh:249
void setLocalKeys(std::vector< LocalKey > &keys) const
Definition: raviartthomassimplexinterpolation.hh:266
Definition: interpolationhelper.hh:20
Definition: interpolationhelper.hh:22
Definition: interpolationhelper.hh:85
Definition: polynomialbasis.hh:63
unsigned int size() const
Definition: polynomialbasis.hh:111