3#ifndef DUNE_LOCALFUNCTIONS_NEDELEC_NEDELECSIMPLEX_NEDELECSIMPLEXINTERPOLATION_HH
4#define DUNE_LOCALFUNCTIONS_NEDELEC_NEDELECSIMPLEX_NEDELECSIMPLEXINTERPOLATION_HH
10#include <dune/common/exceptions.hh>
12#include <dune/geometry/quadraturerules.hh>
13#include <dune/geometry/referenceelements.hh>
14#include <dune/geometry/type.hh>
27 template <
unsigned int dim,
class Field >
28 struct NedelecL2InterpolationFactory;
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 GeometryType gt = geometryId;
88 return gt.isTriangle() || gt.isTetrahedron() ;
107 template<
unsigned int dim,
class Field >
125 typedef FieldVector< Field, dimension >
Tangent;
128 typedef FieldVector< Field, dimension >
Normal;
129 typedef std::array<FieldVector< Field, dimension >,dim-1>
FaceTangents;
139 for( FaceStructure &f : faceStructure_ )
141 for( EdgeStructure& e : edgeStructure_ )
147 return geometry_.id();
163 return numberOfFaces_;
169 return numberOfEdges_;
182 return faceStructure_[ f ].basis_;
189 return edgeStructure_[ e ].basis_;
195 return edgeStructure_[ e ].tangent_;
201 return faceStructure_[ f ].faceTangents_;
207 return faceStructure_[ f ].normal_;
210 template< GeometryType::Id geometryId >
213 constexpr GeometryType geometry = geometryId;
215 geometry_ = geometry;
230 int requiredOrder =
static_cast<int>(
dimension==3);
231 testBasis_ = (
order > requiredOrder ? TestBasisFactory::template create< geometry >(
order-1-requiredOrder ) :
nullptr);
233 const auto &refElement = ReferenceElements< Field, dimension >::general(
type() );
235 numberOfFaces_ = refElement.size( 1 );
236 faceStructure_.reserve( numberOfFaces_ );
239 for (std::size_t i=0; i<numberOfFaces_; i++)
241 FieldVector<Field,dimension> zero(0);
246 auto vertices = refElement.subEntities(i,1,dim).begin();
247 auto vertex1 = *vertices;
248 for(
int j=1; j<dim;j++)
250 auto vertex2 = vertices[j];
252 faceTangents[j-1] = refElement.position(vertex2,dim) - refElement.position(vertex1,dim);
273 TestFaceBasis *faceBasis = ( dim == 3 &&
order > 0 ? Impl::IfGeometryType< CreateFaceBasis, dimension-1 >::apply( refElement.type( i, 1 ),
order-1 ) :
nullptr);
274 faceStructure_.emplace_back( faceBasis, refElement.integrationOuterNormal(i),
faceTangents );
276 assert( faceStructure_.size() == numberOfFaces_ );
278 numberOfEdges_ = refElement.size( dim-1 );
279 edgeStructure_.reserve( numberOfEdges_ );
282 for (std::size_t i=0; i<numberOfEdges_; i++)
284 auto vertexIterator = refElement.subEntities(i,dim-1,dim).begin();
285 auto v0 = *vertexIterator;
286 auto v1 = *(++vertexIterator);
292 auto tangent = std::move(refElement.position(v1,dim) - refElement.position(v0,dim));
294 TestEdgeBasis *edgeBasis = Impl::IfGeometryType< CreateEdgeBasis, 1 >::apply( refElement.type( i, dim-1 ),
order );
295 edgeStructure_.emplace_back( edgeBasis, tangent );
297 assert( edgeStructure_.size() == numberOfEdges_ );
306 : basis_( teb ), tangent_( t )
310 const Dune::FieldVector< Field, dimension > tangent_;
313 template< GeometryType::Id edgeGeometryId >
314 struct CreateEdgeBasis
316 static TestEdgeBasis *apply ( std::size_t
order ) {
return TestEdgeBasisFactory::template create< edgeGeometryId >(
order ); }
327 const Dune::FieldVector< Field, dimension > normal_;
331 template< GeometryType::Id faceGeometryId >
332 struct CreateFaceBasis
334 static TestFaceBasis *apply ( std::size_t
order ) {
return TestFaceBasisFactory::template create< faceGeometryId >(
order ); }
338 std::vector< FaceStructure > faceStructure_;
339 unsigned int numberOfFaces_;
340 std::vector< EdgeStructure > edgeStructure_;
341 unsigned int numberOfEdges_;
342 GeometryType geometry_;
356 template<
unsigned int dimension,
class F>
373 template<
class Function,
class Vector >
374 auto interpolate (
const Function &function, Vector &coefficients )
const
375 -> std::enable_if_t< std::is_same< decltype(std::declval<Vector>().resize(1) ),
void >::value,
void>
377 coefficients.resize(
size());
382 template<
class Basis,
class Matrix >
384 -> std::enable_if_t< std::is_same<
385 decltype(std::declval<Matrix>().rowPtr(0)),
typename Matrix::Field* >::value,
void>
387 matrix.resize(
size(), basis.size() );
401 template <GeometryType::Id geometryId>
422 unsigned int row = 0;
427 keys[row] =
LocalKey(e,dimension-1,i);
439 assert( row ==
size() );
443 template<
class Func,
class Container,
bool type >
448 std::vector<Field> testBasisVal;
450 for (
unsigned int i=0; i<
size(); ++i)
451 for (
unsigned int j=0; j<func.size(); ++j)
454 unsigned int row = 0;
457 typedef Dune::QuadratureRule<Field, 1> EdgeQuadrature;
458 typedef Dune::QuadratureRules<Field, 1> EdgeQuadratureRules;
460 const auto &refElement = Dune::ReferenceElements< Field, dimension >::general( geoType );
468 const auto &geometry = refElement.template geometry< dimension-1 >( e );
469 const Dune::GeometryType subGeoType( geometry.type().id(), 1 );
470 const EdgeQuadrature &edgeQuad = EdgeQuadratureRules::rule( subGeoType, 2*
order_+2 );
472 const unsigned int quadratureSize = edgeQuad.size();
473 for(
unsigned int qi = 0; qi < quadratureSize; ++qi )
478 testBasisVal[0] = 1.;
481 func.evaluate( geometry.global( edgeQuad[qi].position() ) ),
483 edgeQuad[qi].weight(),
491 typedef Dune::QuadratureRule<
Field, dimension-1> FaceQuadrature;
492 typedef Dune::QuadratureRules<
Field, dimension-1> FaceQuadratureRules;
500 const auto &geometry = refElement.template geometry< 1 >( f );
501 const Dune::GeometryType subGeoType( geometry.type().id(), dimension-1 );
502 const FaceQuadrature &faceQuad = FaceQuadratureRules::rule( subGeoType, 2*
order_+2 );
504 const unsigned int quadratureSize = faceQuad.size();
505 for(
unsigned int qi = 0; qi < quadratureSize; ++qi )
510 testBasisVal[0] = 1.;
512 computeFaceDofs( row,
514 func.evaluate( geometry.global( faceQuad[qi].position() ) ),
517 faceQuad[qi].weight(),
530 typedef Dune::QuadratureRule<Field, dimension> Quadrature;
531 typedef Dune::QuadratureRules<Field, dimension> QuadratureRules;
532 const Quadrature &elemQuad = QuadratureRules::rule( geoType, 2*
order_+1 );
534 const unsigned int quadratureSize = elemQuad.size();
535 for(
unsigned int qi = 0; qi < quadratureSize; ++qi )
537 builder_.
testBasis()->template evaluate<0>(elemQuad[qi].position(),testBasisVal);
538 computeInteriorDofs(row,
540 func.evaluate(elemQuad[qi].position()),
541 elemQuad[qi].weight(),
560 template <
class MVal,
class NedVal,
class Matrix>
561 void computeEdgeDofs (
unsigned int startRow,
563 const NedVal &nedVal,
564 const FieldVector<Field,dimension> &tangent,
566 Matrix &matrix)
const
568 const unsigned int endRow = startRow+mVal.size();
569 typename NedVal::const_iterator nedIter = nedVal.begin();
570 for (
unsigned int col = 0; col < nedVal.size() ; ++nedIter,++col)
572 Field cFactor = (*nedIter)*tangent;
573 typename MVal::const_iterator mIter = mVal.begin();
574 for (
unsigned int row = startRow; row!=endRow; ++mIter, ++row )
575 matrix.add(row,col, (weight*cFactor)*(*mIter) );
577 assert( mIter == mVal.end() );
591 template <
class MVal,
class NedVal,
class Matrix>
592 void computeFaceDofs (
unsigned int startRow,
594 const NedVal &nedVal,
596 const FieldVector<Field,dimension> &normal,
598 Matrix &matrix)
const
600 const unsigned int endRow = startRow+mVal.size()*(dimension-1);
601 typename NedVal::const_iterator nedIter = nedVal.begin();
602 for (
unsigned int col = 0; col < nedVal.size() ; ++nedIter,++col)
604 auto const& u=*nedIter;
605 auto const& n=normal;
606 FieldVector<Field,dimension> nedTimesNormal = { u[1]*n[2]-u[2]*n[1],
608 u[0]*n[1]-u[1]*n[0]};
609 typename MVal::const_iterator mIter = mVal.begin();
610 for (
unsigned int row = startRow; row!=endRow; ++mIter)
612 for(
int i=0; i<dimension-1;i++)
614 auto test = *mIter*faceTangents[i];
615 matrix.add(row+i,col, weight*(nedTimesNormal*test) );
620 assert( mIter == mVal.end() );
632 template <
class MVal,
class NedVal,
class Matrix>
633 void computeInteriorDofs (
unsigned int startRow,
635 const NedVal &nedVal,
637 Matrix &matrix)
const
639 const unsigned int endRow = startRow+mVal.size()*dimension;
640 typename NedVal::const_iterator nedIter = nedVal.begin();
641 for (
unsigned int col = 0; col < nedVal.size() ; ++nedIter,++col)
643 typename MVal::const_iterator mIter = mVal.begin();
644 for (
unsigned int row = startRow; row!=endRow; ++mIter,row+=dimension )
645 for (
unsigned int i=0; i<dimension; ++i)
646 matrix.add(row+i,col, (weight*(*mIter))*(*nedIter)[i] );
648 assert( mIter == mVal.end() );
658 template <
unsigned int dim,
class Field >
666 template <GeometryType::Id geometryId>
669 if ( !supports<geometryId>(key) )
672 interpol->template build<geometryId>(key);
676 template <GeometryType::Id geometryId>
679 GeometryType gt = geometryId;
680 return gt.isTriangle() || gt.isTetrahedron() ;
Definition: bdfmcube.hh:16
Describe position of one degree of freedom.
Definition: localkey.hh:21
Definition: nedelecsimplexinterpolation.hh:660
static Object * create(const Key &key)
Definition: nedelecsimplexinterpolation.hh:667
const NedelecL2Interpolation< dim, Field > Object
Definition: nedelecsimplexinterpolation.hh:662
NedelecL2InterpolationBuilder< dim, Field > Builder
Definition: nedelecsimplexinterpolation.hh:661
std::size_t Key
Definition: nedelecsimplexinterpolation.hh:663
static bool supports(const Key &key)
Definition: nedelecsimplexinterpolation.hh:677
std::remove_const< Object >::type NonConstObject
Definition: nedelecsimplexinterpolation.hh:664
static void release(Object *object)
Definition: nedelecsimplexinterpolation.hh:682
Definition: nedelecsimplexinterpolation.hh:36
LocalCoefficientsContainer(const Setter &setter)
Definition: nedelecsimplexinterpolation.hh:41
const LocalKey & localKey(const unsigned int i) const
Definition: nedelecsimplexinterpolation.hh:46
std::size_t size() const
Definition: nedelecsimplexinterpolation.hh:52
Definition: nedelecsimplexinterpolation.hh:68
static Object * create(const Key &key)
Definition: nedelecsimplexinterpolation.hh:73
static bool supports(const Key &key)
Definition: nedelecsimplexinterpolation.hh:85
const LocalCoefficientsContainer Object
Definition: nedelecsimplexinterpolation.hh:70
std::size_t Key
Definition: nedelecsimplexinterpolation.hh:69
static void release(Object *object)
Definition: nedelecsimplexinterpolation.hh:90
Definition: nedelecsimplexinterpolation.hh:109
TestEdgeBasis * testEdgeBasis(unsigned int e) const
Definition: nedelecsimplexinterpolation.hh:186
~NedelecL2InterpolationBuilder()
Definition: nedelecsimplexinterpolation.hh:136
GeometryType type() const
Definition: nedelecsimplexinterpolation.hh:150
TestBasisFactory::Object TestBasis
Definition: nedelecsimplexinterpolation.hh:114
FieldVector< Field, dimension > Tangent
Definition: nedelecsimplexinterpolation.hh:125
TestFaceBasisFactory::Object TestFaceBasis
Definition: nedelecsimplexinterpolation.hh:118
TestFaceBasis * testFaceBasis(unsigned int f) const
Definition: nedelecsimplexinterpolation.hh:179
TestEdgeBasisFactory::Object TestEdgeBasis
Definition: nedelecsimplexinterpolation.hh:122
FieldVector< Field, dimension > Normal
Definition: nedelecsimplexinterpolation.hh:128
void build(std::size_t order)
Definition: nedelecsimplexinterpolation.hh:211
OrthonormalBasisFactory< dimension, Field > TestBasisFactory
Definition: nedelecsimplexinterpolation.hh:113
OrthonormalBasisFactory< dimension-1, Field > TestFaceBasisFactory
Definition: nedelecsimplexinterpolation.hh:117
const FaceTangents & faceTangents(unsigned int f) const
Definition: nedelecsimplexinterpolation.hh:198
unsigned int faceSize() const
Definition: nedelecsimplexinterpolation.hh:161
TestBasis * testBasis() const
Definition: nedelecsimplexinterpolation.hh:173
std::array< FieldVector< Field, dimension >, dim-1 > FaceTangents
Definition: nedelecsimplexinterpolation.hh:129
OrthonormalBasisFactory< 1, Field > TestEdgeBasisFactory
Definition: nedelecsimplexinterpolation.hh:121
const Tangent & edgeTangent(unsigned int e) const
Definition: nedelecsimplexinterpolation.hh:192
NedelecL2InterpolationBuilder(NedelecL2InterpolationBuilder &&)=delete
std::size_t order() const
Definition: nedelecsimplexinterpolation.hh:155
unsigned int edgeSize() const
Definition: nedelecsimplexinterpolation.hh:167
unsigned int topologyId() const
Definition: nedelecsimplexinterpolation.hh:145
NedelecL2InterpolationBuilder(const NedelecL2InterpolationBuilder &)=delete
static const unsigned int dimension
Definition: nedelecsimplexinterpolation.hh:110
NedelecL2InterpolationBuilder()=default
const Normal & normal(unsigned int f) const
Definition: nedelecsimplexinterpolation.hh:204
An L2-based interpolation for Nedelec.
Definition: nedelecsimplexinterpolation.hh:359
Builder::FaceTangents FaceTangents
Definition: nedelecsimplexinterpolation.hh:366
F Field
Definition: nedelecsimplexinterpolation.hh:364
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: nedelecsimplexinterpolation.hh:374
std::size_t size() const
Definition: nedelecsimplexinterpolation.hh:396
void interpolate(typename Base::template Helper< Func, Container, type > &func) const
Definition: nedelecsimplexinterpolation.hh:444
std::size_t order_
Definition: nedelecsimplexinterpolation.hh:654
NedelecL2InterpolationBuilder< dimension, Field > Builder
Definition: nedelecsimplexinterpolation.hh:365
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: nedelecsimplexinterpolation.hh:383
std::size_t size_
Definition: nedelecsimplexinterpolation.hh:655
NedelecL2Interpolation()
Definition: nedelecsimplexinterpolation.hh:368
void build(std::size_t order)
Definition: nedelecsimplexinterpolation.hh:402
std::size_t order() const
Definition: nedelecsimplexinterpolation.hh:392
void setLocalKeys(std::vector< LocalKey > &keys) const
Definition: nedelecsimplexinterpolation.hh:419
Builder builder_
Definition: nedelecsimplexinterpolation.hh:653
Definition: orthonormalbasis.hh:18
static void release(Object *object)
Definition: orthonormalbasis.hh:55
Definition: interpolationhelper.hh:20
Definition: interpolationhelper.hh:22
Definition: polynomialbasis.hh:63
unsigned int size() const
Definition: polynomialbasis.hh:111