5#ifndef DUNE_LOCALFUNCTIONS_NEDELEC_NEDELECSIMPLEX_NEDELECSIMPLEXINTERPOLATION_HH
6#define DUNE_LOCALFUNCTIONS_NEDELEC_NEDELECSIMPLEX_NEDELECSIMPLEXINTERPOLATION_HH
12#include <dune/common/exceptions.hh>
14#include <dune/geometry/quadraturerules.hh>
15#include <dune/geometry/referenceelements.hh>
16#include <dune/geometry/type.hh>
32 template <
unsigned int dim,
class Field >
33 struct NedelecL2InterpolationFactory;
45 template <
class Setter>
48 setter.setLocalKeys(localKey_);
54 return localKey_[ i ];
59 return localKey_.size();
63 std::vector< LocalKey > localKey_;
71 template <
unsigned int dim >
74 typedef std::size_t
Key;
77 template< GeometryType::Id geometryId >
81 if( !supports< geometryId >( key ) )
83 typename InterpolationFactory::Object *interpolation = InterpolationFactory::template create< geometryId >( key );
85 InterpolationFactory::release( interpolation );
89 template< GeometryType::Id geometryId >
92 GeometryType gt = geometryId;
93 return gt.isTriangle() || gt.isTetrahedron() ;
112 template<
unsigned int dim,
class Field >
130 typedef FieldVector< Field, dimension >
Tangent;
133 typedef FieldVector< Field, dimension >
Normal;
134 typedef std::array<FieldVector< Field, dimension >,dim-1>
FaceTangents;
144 for( FaceStructure &f : faceStructure_ )
146 for( EdgeStructure& e : edgeStructure_ )
152 return geometry_.id();
168 return numberOfFaces_;
174 return numberOfEdges_;
187 return faceStructure_[ f ].basis_;
194 return edgeStructure_[ e ].basis_;
200 return edgeStructure_[ e ].tangent_;
206 return faceStructure_[ f ].faceTangents_;
212 return faceStructure_[ f ].normal_;
215 template< GeometryType::Id geometryId >
218 constexpr GeometryType geometry = geometryId;
220 geometry_ = geometry;
235 int requiredOrder =
static_cast<int>(
dimension==3);
236 testBasis_ = (
order > requiredOrder ? TestBasisFactory::template create< geometry >(
order-1-requiredOrder ) :
nullptr);
238 const auto &refElement = ReferenceElements< Field, dimension >::general(
type() );
240 numberOfFaces_ = refElement.size( 1 );
241 faceStructure_.reserve( numberOfFaces_ );
244 for (std::size_t i=0; i<numberOfFaces_; i++)
246 FieldVector<Field,dimension> zero(0);
251 auto vertices = refElement.subEntities(i,1,dim).begin();
252 auto vertex1 = *vertices;
253 for(
int j=1; j<dim;j++)
255 auto vertex2 = vertices[j];
257 faceTangents[j-1] = refElement.position(vertex2,dim) - refElement.position(vertex1,dim);
278 TestFaceBasis *faceBasis = ( dim == 3 &&
order > 0 ? Impl::IfGeometryType< CreateFaceBasis, dimension-1 >::apply( refElement.type( i, 1 ),
order-1 ) :
nullptr);
279 faceStructure_.emplace_back( faceBasis, refElement.integrationOuterNormal(i),
faceTangents );
281 assert( faceStructure_.size() == numberOfFaces_ );
283 numberOfEdges_ = refElement.size( dim-1 );
284 edgeStructure_.reserve( numberOfEdges_ );
287 for (std::size_t i=0; i<numberOfEdges_; i++)
289 auto vertexIterator = refElement.subEntities(i,dim-1,dim).begin();
290 auto v0 = *vertexIterator;
291 auto v1 = *(++vertexIterator);
297 auto tangent = refElement.position(v1,dim) - refElement.position(v0,dim);
299 TestEdgeBasis *edgeBasis = Impl::IfGeometryType< CreateEdgeBasis, 1 >::apply( refElement.type( i, dim-1 ),
order );
300 edgeStructure_.emplace_back( edgeBasis, tangent );
302 assert( edgeStructure_.size() == numberOfEdges_ );
311 : basis_( teb ), tangent_( t )
315 const Dune::FieldVector< Field, dimension > tangent_;
318 template< GeometryType::Id edgeGeometryId >
319 struct CreateEdgeBasis
321 static TestEdgeBasis *apply ( std::size_t
order ) {
return TestEdgeBasisFactory::template create< edgeGeometryId >(
order ); }
332 const Dune::FieldVector< Field, dimension > normal_;
336 template< GeometryType::Id faceGeometryId >
337 struct CreateFaceBasis
339 static TestFaceBasis *apply ( std::size_t
order ) {
return TestFaceBasisFactory::template create< faceGeometryId >(
order ); }
343 std::vector< FaceStructure > faceStructure_;
344 unsigned int numberOfFaces_;
345 std::vector< EdgeStructure > edgeStructure_;
346 unsigned int numberOfEdges_;
347 GeometryType geometry_;
361 template<
unsigned int dimension,
class F>
378 template<
class Function,
class Vector,
379 decltype(std::declval<Vector>().size(),
bool{}) =
true,
380 decltype(std::declval<Vector>().resize(0u),
bool{}) =
true>
381 void interpolate (
const Function &function, Vector &coefficients )
const
383 coefficients.resize(
size());
388 template<
class Basis,
class Matrix,
389 decltype(std::declval<Matrix>().rows(),
bool{}) =
true,
390 decltype(std::declval<Matrix>().cols(),
bool{}) =
true,
391 decltype(std::declval<Matrix>().resize(0u,0u),
bool{}) =
true>
394 matrix.resize(
size(), basis.size() );
408 template <GeometryType::Id geometryId>
429 unsigned int row = 0;
434 keys[row] =
LocalKey(e,dimension-1,i);
446 assert( row ==
size() );
450 template<
class Func,
class Container,
bool type >
455 std::vector<Field> testBasisVal;
457 for (
unsigned int i=0; i<
size(); ++i)
458 for (
unsigned int j=0; j<func.size(); ++j)
461 unsigned int row = 0;
464 typedef Dune::QuadratureRule<Field, 1> EdgeQuadrature;
465 typedef Dune::QuadratureRules<Field, 1> EdgeQuadratureRules;
467 const auto &refElement = Dune::ReferenceElements< Field, dimension >::general( geoType );
475 const auto &geometry = refElement.template geometry< dimension-1 >( e );
476 const Dune::GeometryType subGeoType( geometry.type().id(), 1 );
477 const EdgeQuadrature &edgeQuad = EdgeQuadratureRules::rule( subGeoType, 2*
order_+2 );
479 const unsigned int quadratureSize = edgeQuad.size();
480 for(
unsigned int qi = 0; qi < quadratureSize; ++qi )
485 testBasisVal[0] = 1.;
488 func.evaluate( geometry.global( edgeQuad[qi].position() ) ),
490 edgeQuad[qi].weight(),
498 typedef Dune::QuadratureRule<
Field, dimension-1> FaceQuadrature;
499 typedef Dune::QuadratureRules<
Field, dimension-1> FaceQuadratureRules;
507 const auto &geometry = refElement.template geometry< 1 >( f );
508 const Dune::GeometryType subGeoType( geometry.type().id(), dimension-1 );
509 const FaceQuadrature &faceQuad = FaceQuadratureRules::rule( subGeoType, 2*
order_+2 );
511 const unsigned int quadratureSize = faceQuad.size();
512 for(
unsigned int qi = 0; qi < quadratureSize; ++qi )
517 testBasisVal[0] = 1.;
519 computeFaceDofs( row,
521 func.evaluate( geometry.global( faceQuad[qi].position() ) ),
524 faceQuad[qi].weight(),
537 typedef Dune::QuadratureRule<Field, dimension> Quadrature;
538 typedef Dune::QuadratureRules<Field, dimension> QuadratureRules;
539 const Quadrature &elemQuad = QuadratureRules::rule( geoType, 2*
order_+1 );
541 const unsigned int quadratureSize = elemQuad.size();
542 for(
unsigned int qi = 0; qi < quadratureSize; ++qi )
544 builder_.
testBasis()->template evaluate<0>(elemQuad[qi].position(),testBasisVal);
545 computeInteriorDofs(row,
547 func.evaluate(elemQuad[qi].position()),
548 elemQuad[qi].weight(),
567 template <
class MVal,
class NedVal,
class Matrix>
568 void computeEdgeDofs (
unsigned int startRow,
570 const NedVal &nedVal,
571 const FieldVector<Field,dimension> &tangent,
573 Matrix &matrix)
const
575 const unsigned int endRow = startRow+mVal.size();
576 typename NedVal::const_iterator nedIter = nedVal.begin();
577 for (
unsigned int col = 0; col < nedVal.size() ; ++nedIter,++col)
579 Field cFactor = (*nedIter)*tangent;
580 typename MVal::const_iterator mIter = mVal.begin();
581 for (
unsigned int row = startRow; row!=endRow; ++mIter, ++row )
582 matrix.add(row,col, (weight*cFactor)*(*mIter) );
584 assert( mIter == mVal.end() );
598 template <
class MVal,
class NedVal,
class Matrix>
599 void computeFaceDofs (
unsigned int startRow,
601 const NedVal &nedVal,
603 const FieldVector<Field,dimension> &normal,
605 Matrix &matrix)
const
607 const unsigned int endRow = startRow+mVal.size()*(dimension-1);
608 typename NedVal::const_iterator nedIter = nedVal.begin();
609 for (
unsigned int col = 0; col < nedVal.size() ; ++nedIter,++col)
611 auto const& u=*nedIter;
612 auto const& n=normal;
613 FieldVector<Field,dimension> nedTimesNormal = { u[1]*n[2]-u[2]*n[1],
615 u[0]*n[1]-u[1]*n[0]};
616 typename MVal::const_iterator mIter = mVal.begin();
617 for (
unsigned int row = startRow; row!=endRow; ++mIter)
619 for(
int i=0; i<dimension-1;i++)
621 auto test = *mIter*faceTangents[i];
622 matrix.add(row+i,col, weight*(nedTimesNormal*test) );
627 assert( mIter == mVal.end() );
639 template <
class MVal,
class NedVal,
class Matrix>
640 void computeInteriorDofs (
unsigned int startRow,
642 const NedVal &nedVal,
644 Matrix &matrix)
const
646 const unsigned int endRow = startRow+mVal.size()*dimension;
647 typename NedVal::const_iterator nedIter = nedVal.begin();
648 for (
unsigned int col = 0; col < nedVal.size() ; ++nedIter,++col)
650 typename MVal::const_iterator mIter = mVal.begin();
651 for (
unsigned int row = startRow; row!=endRow; ++mIter,row+=dimension )
652 for (
unsigned int i=0; i<dimension; ++i)
653 matrix.add(row+i,col, (weight*(*mIter))*(*nedIter)[i] );
655 assert( mIter == mVal.end() );
665 template <
unsigned int dim,
class Field >
673 template <GeometryType::Id geometryId>
676 if ( !supports<geometryId>(key) )
679 interpol->template build<geometryId>(key);
683 template <GeometryType::Id geometryId>
686 GeometryType gt = geometryId;
687 return gt.isTriangle() || gt.isTetrahedron() ;
Definition bdfmcube.hh:18
Describe position of one degree of freedom.
Definition localkey.hh:24
Definition nedelecsimplexinterpolation.hh:667
static Object * create(const Key &key)
Definition nedelecsimplexinterpolation.hh:674
const NedelecL2Interpolation< dim, Field > Object
Definition nedelecsimplexinterpolation.hh:669
NedelecL2InterpolationBuilder< dim, Field > Builder
Definition nedelecsimplexinterpolation.hh:668
std::size_t Key
Definition nedelecsimplexinterpolation.hh:670
static bool supports(const Key &key)
Definition nedelecsimplexinterpolation.hh:684
std::remove_const< Object >::type NonConstObject
Definition nedelecsimplexinterpolation.hh:671
static void release(Object *object)
Definition nedelecsimplexinterpolation.hh:689
Definition nedelecsimplexinterpolation.hh:41
LocalCoefficientsContainer(const Setter &setter)
Definition nedelecsimplexinterpolation.hh:46
const LocalKey & localKey(const unsigned int i) const
Definition nedelecsimplexinterpolation.hh:51
std::size_t size() const
Definition nedelecsimplexinterpolation.hh:57
Definition nedelecsimplexinterpolation.hh:73
static Object * create(const Key &key)
Definition nedelecsimplexinterpolation.hh:78
static bool supports(const Key &key)
Definition nedelecsimplexinterpolation.hh:90
const LocalCoefficientsContainer Object
Definition nedelecsimplexinterpolation.hh:75
std::size_t Key
Definition nedelecsimplexinterpolation.hh:74
static void release(Object *object)
Definition nedelecsimplexinterpolation.hh:95
Definition nedelecsimplexinterpolation.hh:114
TestEdgeBasis * testEdgeBasis(unsigned int e) const
Definition nedelecsimplexinterpolation.hh:191
~NedelecL2InterpolationBuilder()
Definition nedelecsimplexinterpolation.hh:141
GeometryType type() const
Definition nedelecsimplexinterpolation.hh:155
TestBasisFactory::Object TestBasis
Definition nedelecsimplexinterpolation.hh:119
FieldVector< Field, dimension > Tangent
Definition nedelecsimplexinterpolation.hh:130
TestFaceBasisFactory::Object TestFaceBasis
Definition nedelecsimplexinterpolation.hh:123
TestFaceBasis * testFaceBasis(unsigned int f) const
Definition nedelecsimplexinterpolation.hh:184
TestEdgeBasisFactory::Object TestEdgeBasis
Definition nedelecsimplexinterpolation.hh:127
FieldVector< Field, dimension > Normal
Definition nedelecsimplexinterpolation.hh:133
void build(std::size_t order)
Definition nedelecsimplexinterpolation.hh:216
OrthonormalBasisFactory< dimension, Field > TestBasisFactory
Definition nedelecsimplexinterpolation.hh:118
OrthonormalBasisFactory< dimension-1, Field > TestFaceBasisFactory
Definition nedelecsimplexinterpolation.hh:122
const FaceTangents & faceTangents(unsigned int f) const
Definition nedelecsimplexinterpolation.hh:203
unsigned int faceSize() const
Definition nedelecsimplexinterpolation.hh:166
TestBasis * testBasis() const
Definition nedelecsimplexinterpolation.hh:178
std::array< FieldVector< Field, dimension >, dim-1 > FaceTangents
Definition nedelecsimplexinterpolation.hh:134
OrthonormalBasisFactory< 1, Field > TestEdgeBasisFactory
Definition nedelecsimplexinterpolation.hh:126
const Tangent & edgeTangent(unsigned int e) const
Definition nedelecsimplexinterpolation.hh:197
NedelecL2InterpolationBuilder(NedelecL2InterpolationBuilder &&)=delete
std::size_t order() const
Definition nedelecsimplexinterpolation.hh:160
unsigned int edgeSize() const
Definition nedelecsimplexinterpolation.hh:172
unsigned int topologyId() const
Definition nedelecsimplexinterpolation.hh:150
NedelecL2InterpolationBuilder(const NedelecL2InterpolationBuilder &)=delete
static const unsigned int dimension
Definition nedelecsimplexinterpolation.hh:115
NedelecL2InterpolationBuilder()=default
const Normal & normal(unsigned int f) const
Definition nedelecsimplexinterpolation.hh:209
An L2-based interpolation for Nedelec.
Definition nedelecsimplexinterpolation.hh:364
Builder::FaceTangents FaceTangents
Definition nedelecsimplexinterpolation.hh:371
F Field
Definition nedelecsimplexinterpolation.hh:369
std::size_t size() const
Definition nedelecsimplexinterpolation.hh:403
void interpolate(typename Base::template Helper< Func, Container, type > &func) const
Definition nedelecsimplexinterpolation.hh:451
std::size_t order_
Definition nedelecsimplexinterpolation.hh:661
NedelecL2InterpolationBuilder< dimension, Field > Builder
Definition nedelecsimplexinterpolation.hh:370
std::size_t size_
Definition nedelecsimplexinterpolation.hh:662
NedelecL2Interpolation()
Definition nedelecsimplexinterpolation.hh:373
void interpolate(const Function &function, Vector &coefficients) const
Definition nedelecsimplexinterpolation.hh:381
void build(std::size_t order)
Definition nedelecsimplexinterpolation.hh:409
std::size_t order() const
Definition nedelecsimplexinterpolation.hh:399
void setLocalKeys(std::vector< LocalKey > &keys) const
Definition nedelecsimplexinterpolation.hh:426
Builder builder_
Definition nedelecsimplexinterpolation.hh:660
Definition orthonormalbasis.hh:20
static void release(Object *object)
Definition orthonormalbasis.hh:57
Definition interpolationhelper.hh:21
Definition interpolationhelper.hh:23
Definition polynomialbasis.hh:65
unsigned int size() const
Definition polynomialbasis.hh:113