5#ifndef DUNE_LOCALFUNCTIONS_NEDELEC_NEDELECSIMPLEX_NEDELECSIMPLEXINTERPOLATION_HH
6#define DUNE_LOCALFUNCTIONS_NEDELEC_NEDELECSIMPLEX_NEDELECSIMPLEXINTERPOLATION_HH
15#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>
32 template <
unsigned int dim,
class Field >
33 struct NedelecL2InterpolationFactory;
40 class LocalCoefficientsContainer
42 typedef LocalCoefficientsContainer This;
45 template <
class Setter>
46 LocalCoefficientsContainer (
const Setter &setter )
48 setter.setLocalKeys(localKey_);
51 const LocalKey &localKey (
const unsigned int i )
const
54 return localKey_[ i ];
57 std::size_t
size ()
const
59 return localKey_.size();
63 std::vector< LocalKey > localKey_;
71 template <
unsigned int dim >
72 struct NedelecCoefficientsFactory
74 typedef std::size_t Key;
75 typedef const LocalCoefficientsContainer Object;
77 template< GeometryType::Id geometryId >
78 static Object *create(
const Key &key )
80 typedef NedelecL2InterpolationFactory< dim, double > InterpolationFactory;
81 if( !supports< geometryId >( key ) )
83 typename InterpolationFactory::Object *interpolation = InterpolationFactory::template create< geometryId >( key );
84 Object *localKeys =
new Object( *interpolation );
85 InterpolationFactory::release( interpolation );
89 template< GeometryType::Id geometryId >
90 static bool supports (
const Key &key )
93 return gt.isTriangle() ||
gt.isTetrahedron() ;
95 static void release( Object *
object ) {
delete object; }
112 template<
unsigned int dim,
class Field >
113 struct NedelecL2InterpolationBuilder
115 static const unsigned int dimension = dim;
118 typedef OrthonormalBasisFactory< dimension, Field > TestBasisFactory;
119 typedef typename TestBasisFactory::Object TestBasis;
122 typedef OrthonormalBasisFactory< dimension-1, Field > TestFaceBasisFactory;
123 typedef typename TestFaceBasisFactory::Object TestFaceBasis;
126 typedef OrthonormalBasisFactory< 1, Field > TestEdgeBasisFactory;
127 typedef typename TestEdgeBasisFactory::Object TestEdgeBasis;
134 typedef std::array<FieldVector< Field, dimension >,dim-1> FaceTangents;
136 NedelecL2InterpolationBuilder () =
default;
138 NedelecL2InterpolationBuilder (
const NedelecL2InterpolationBuilder & ) =
delete;
139 NedelecL2InterpolationBuilder ( NedelecL2InterpolationBuilder && ) =
delete;
141 ~NedelecL2InterpolationBuilder ()
143 TestBasisFactory::release( testBasis_ );
144 for( FaceStructure &f : faceStructure_ )
145 TestFaceBasisFactory::release( f.basis_ );
146 for( EdgeStructure& e : edgeStructure_ )
147 TestEdgeBasisFactory::release( e.basis_ );
150 unsigned int topologyId ()
const
152 return geometry_.id();
160 std::size_t order ()
const
166 unsigned int faceSize ()
const
168 return numberOfFaces_;
172 unsigned int edgeSize ()
const
174 return numberOfEdges_;
178 TestBasis *testBasis ()
const
184 TestFaceBasis *testFaceBasis (
unsigned int f )
const
186 assert( f < faceSize() );
187 return faceStructure_[ f ].basis_;
191 TestEdgeBasis *testEdgeBasis (
unsigned int e )
const
193 assert( e < edgeSize() );
194 return edgeStructure_[ e ].basis_;
197 const Tangent& edgeTangent (
unsigned int e )
const
199 assert( e < edgeSize() );
200 return edgeStructure_[ e ].tangent_;
203 const FaceTangents& faceTangents (
unsigned int f )
const
205 assert( f < faceSize() );
206 return faceStructure_[ f ].faceTangents_;
209 const Normal &normal (
unsigned int f )
const
211 assert( f < faceSize() );
212 return faceStructure_[ f ].normal_;
215 template< GeometryType::Id geometryId >
216 void build ( std::size_t order )
220 geometry_ = geometry;
235 int requiredOrder =
static_cast<int>(dimension==3);
236 testBasis_ = (order > requiredOrder ? TestBasisFactory::template create< geometry >( order-1-requiredOrder ) : nullptr);
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);
247 FaceTangents faceTangents;
248 faceTangents.fill(zero);
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);
262 faceTangents[j-1] *=-1;
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_ );
310 EdgeStructure( TestEdgeBasis *teb,
const Tangent &t )
311 : basis_( teb ), tangent_( t )
314 TestEdgeBasis *basis_;
318 template< GeometryType::Id edgeGeometryId >
319 struct CreateEdgeBasis
321 static TestEdgeBasis *apply ( std::size_t order ) {
return TestEdgeBasisFactory::template create< edgeGeometryId >( order ); }
327 FaceStructure( TestFaceBasis *tfb,
const Normal& normal,
const FaceTangents& faceTangents )
328 : basis_( tfb ), normal_(normal), faceTangents_( faceTangents )
331 TestFaceBasis *basis_;
333 const FaceTangents faceTangents_;
336 template< GeometryType::Id faceGeometryId >
337 struct CreateFaceBasis
339 static TestFaceBasis *apply ( std::size_t order ) {
return TestFaceBasisFactory::template create< faceGeometryId >( order ); }
342 TestBasis *testBasis_ =
nullptr;
343 std::vector< FaceStructure > faceStructure_;
344 unsigned int numberOfFaces_;
345 std::vector< EdgeStructure > edgeStructure_;
346 unsigned int numberOfEdges_;
361 template<
unsigned int dimension,
class F>
363 :
public InterpolationHelper< F ,dimension >
366 typedef InterpolationHelper<F,dimension> Base;
370 typedef NedelecL2InterpolationBuilder<dimension,Field> Builder;
371 typedef typename Builder::FaceTangents FaceTangents;
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());
384 typename Base::template Helper<Function,Vector,true> func( function,coefficients );
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>
392 void interpolate (
const Basis &basis,
Matrix &matrix )
const
394 matrix.resize( size(), basis.size() );
395 typename Base::template Helper<Basis,Matrix,false> func( basis,matrix );
399 std::size_t order()
const
403 std::size_t size()
const
408 template <GeometryType::Id geometryId>
409 void build( std::size_t order )
413 builder_.template build<geometryId>(order_);
414 if (builder_.testBasis())
415 size_ += dimension*builder_.testBasis()->size();
417 for (
unsigned int f=0; f<builder_.faceSize(); ++f )
418 if (builder_.testFaceBasis(f))
419 size_ += (dimension-1)*builder_.testFaceBasis(f)->size();
421 for (
unsigned int e=0; e<builder_.edgeSize(); ++e )
422 if (builder_.testEdgeBasis(e))
423 size_ += builder_.testEdgeBasis(e)->size();
426 void setLocalKeys(std::vector< LocalKey > &keys)
const
429 unsigned int row = 0;
430 for (
unsigned int e=0; e<builder_.edgeSize(); ++e)
432 if (builder_.edgeSize())
433 for (
unsigned int i=0; i<builder_.testEdgeBasis(e)->size(); ++i,++row)
434 keys[row] =
LocalKey(e,dimension-1,i);
436 for (
unsigned int f=0; f<builder_.faceSize(); ++f)
438 if (builder_.faceSize())
439 for (
unsigned int i=0; i<builder_.testFaceBasis(f)->size()*(dimension-1); ++i,++row)
443 if (builder_.testBasis())
444 for (
unsigned int i=0; i<builder_.testBasis()->size()*dimension; ++i,++row)
446 assert( row == size() );
450 template<
class Func,
class Container,
bool type >
451 void interpolate (
typename Base::template Helper<Func,Container,type> &func )
const
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;
469 for (
unsigned int e=0; e<builder_.edgeSize(); ++e)
471 if (!builder_.testEdgeBasis(e))
473 testBasisVal.resize(builder_.testEdgeBasis(e)->size());
475 const auto &geometry = refElement.template geometry< dimension-1 >( e );
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 )
483 builder_.testEdgeBasis(e)->template evaluate<0>(edgeQuad[qi].position(),testBasisVal);
485 testBasisVal[0] = 1.;
488 func.evaluate( geometry.global( edgeQuad[qi].position() ) ),
489 builder_.edgeTangent(e),
490 edgeQuad[qi].weight(),
494 row += builder_.testEdgeBasis(e)->size();
501 for (
unsigned int f=0; f<builder_.faceSize(); ++f)
503 if (builder_.testFaceBasis(f))
505 testBasisVal.resize(builder_.testFaceBasis(f)->size());
507 const auto &geometry = refElement.template geometry< 1 >( f );
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 )
515 builder_.testFaceBasis(f)->template evaluate<0>(faceQuad[qi].position(),testBasisVal);
517 testBasisVal[0] = 1.;
519 computeFaceDofs( row,
521 func.evaluate( geometry.global( faceQuad[qi].position() ) ),
522 builder_.faceTangents(f),
524 faceQuad[qi].weight(),
528 row += builder_.testFaceBasis(f)->size()*(dimension-1);
533 if (builder_.testBasis())
535 testBasisVal.resize(builder_.testBasis()->size());
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(),
552 row += builder_.testBasis()->size()*dimension;
567 template <
class MVal,
class NedVal,
class Matrix>
568 void computeEdgeDofs (
unsigned int startRow,
570 const NedVal &nedVal,
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,
602 const FaceTangents& faceTangents,
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;
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,
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 >
666 struct NedelecL2InterpolationFactory
668 typedef NedelecL2InterpolationBuilder<dim,Field> Builder;
670 typedef std::size_t Key;
671 typedef typename std::remove_const<Object>::type NonConstObject;
673 template <GeometryType::Id geometryId>
674 static Object *create(
const Key &key )
676 if ( !supports<geometryId>(key) )
678 NonConstObject *interpol =
new NonConstObject();
679 interpol->template build<geometryId>(key);
683 template <GeometryType::Id geometryId>
684 static bool supports(
const Key &key )
686 GeometryType
gt = geometryId;
687 return gt.isTriangle() ||
gt.isTetrahedron() ;
689 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
An L2-based interpolation for Nedelec.
Definition: nedelecsimplexinterpolation.hh:364
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
actual interface class for quadratures
A few common exception classes.
GeometryType
Type representing VTK's entity geometry types.
Definition: common.hh:132
bool gt(const T &first, const T &second, typename EpsilonType< T >::Type epsilon)
test if first greater than second
Definition: float_cmp.cc:158
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.