DUNE-FEM (unstable)

nedelecsimplexinterpolation.hh
1// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2// vi: set et ts=4 sw=2 sts=2:
3// SPDX-FileCopyrightInfo: Copyright © DUNE Project contributors, see file LICENSE.md in module root
4// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
5#ifndef DUNE_LOCALFUNCTIONS_NEDELEC_NEDELECSIMPLEX_NEDELECSIMPLEXINTERPOLATION_HH
6#define DUNE_LOCALFUNCTIONS_NEDELEC_NEDELECSIMPLEX_NEDELECSIMPLEXINTERPOLATION_HH
7
8#include <fstream>
9#include <utility>
10#include <numeric>
11
13
15#include <dune/geometry/referenceelements.hh>
16#include <dune/geometry/type.hh>
17
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>
22
23namespace Dune
24{
25
26 // Internal Forward Declarations
27 // -----------------------------
28
32 template < unsigned int dim, class Field >
33 struct NedelecL2InterpolationFactory;
34
35
36
37 // LocalCoefficientsContainer
38 // --------------------------
39
40 class LocalCoefficientsContainer
41 {
42 typedef LocalCoefficientsContainer This;
43
44 public:
45 template <class Setter>
46 LocalCoefficientsContainer ( const Setter &setter )
47 {
48 setter.setLocalKeys(localKey_);
49 }
50
51 const LocalKey &localKey ( const unsigned int i ) const
52 {
53 assert( i < size() );
54 return localKey_[ i ];
55 }
56
57 std::size_t size () const
58 {
59 return localKey_.size();
60 }
61
62 private:
63 std::vector< LocalKey > localKey_;
64 };
65
66
67
68 // NedelecCoefficientsFactory
69 // --------------------------------
70
71 template < unsigned int dim >
72 struct NedelecCoefficientsFactory
73 {
74 typedef std::size_t Key;
75 typedef const LocalCoefficientsContainer Object;
76
77 template< GeometryType::Id geometryId >
78 static Object *create( const Key &key )
79 {
80 typedef NedelecL2InterpolationFactory< dim, double > InterpolationFactory;
81 if( !supports< geometryId >( key ) )
82 return nullptr;
83 typename InterpolationFactory::Object *interpolation = InterpolationFactory::template create< geometryId >( key );
84 Object *localKeys = new Object( *interpolation );
85 InterpolationFactory::release( interpolation );
86 return localKeys;
87 }
88
89 template< GeometryType::Id geometryId >
90 static bool supports ( const Key &key )
91 {
92 GeometryType gt = geometryId;
93 return gt.isTriangle() || gt.isTetrahedron() ;
94 }
95 static void release( Object *object ) { delete object; }
96 };
97
98
99
100 // NedelecL2InterpolationBuilder
101 // ------------------------
102
103 // L2 Interpolation requires:
104 // - for element
105 // - test basis
106 // - for each face (dynamic)
107 // - test basis
108 // - tangents
109 // - for each edge (dynamic)
110 // - test basis
111 // - tangent
112 template< unsigned int dim, class Field >
113 struct NedelecL2InterpolationBuilder
114 {
115 static const unsigned int dimension = dim;
116
117 // for the dofs associated to the element
118 typedef OrthonormalBasisFactory< dimension, Field > TestBasisFactory;
119 typedef typename TestBasisFactory::Object TestBasis;
120
121 // for the dofs associated to the faces
122 typedef OrthonormalBasisFactory< dimension-1, Field > TestFaceBasisFactory;
123 typedef typename TestFaceBasisFactory::Object TestFaceBasis;
124
125 // for the dofs associated to the edges
126 typedef OrthonormalBasisFactory< 1, Field > TestEdgeBasisFactory;
127 typedef typename TestEdgeBasisFactory::Object TestEdgeBasis;
128
129 // the tangent of the edges
130 typedef FieldVector< Field, dimension > Tangent;
131
132 // the normal and the tangents of the faces
133 typedef FieldVector< Field, dimension > Normal;
134 typedef std::array<FieldVector< Field, dimension >,dim-1> FaceTangents;
135
136 NedelecL2InterpolationBuilder () = default;
137
138 NedelecL2InterpolationBuilder ( const NedelecL2InterpolationBuilder & ) = delete;
139 NedelecL2InterpolationBuilder ( NedelecL2InterpolationBuilder && ) = delete;
140
141 ~NedelecL2InterpolationBuilder ()
142 {
143 TestBasisFactory::release( testBasis_ );
144 for( FaceStructure &f : faceStructure_ )
145 TestFaceBasisFactory::release( f.basis_ );
146 for( EdgeStructure& e : edgeStructure_ )
147 TestEdgeBasisFactory::release( e.basis_ );
148 }
149
150 unsigned int topologyId () const
151 {
152 return geometry_.id();
153 }
154
155 GeometryType type () const
156 {
157 return geometry_;
158 }
159
160 std::size_t order () const
161 {
162 return order_;
163 }
164
165 // number of faces
166 unsigned int faceSize () const
167 {
168 return numberOfFaces_;
169 }
170
171 // number of edges
172 unsigned int edgeSize () const
173 {
174 return numberOfEdges_;
175 }
176
177 // basis associated to the element
178 TestBasis *testBasis () const
179 {
180 return testBasis_;
181 }
182
183 // basis associated to face f
184 TestFaceBasis *testFaceBasis ( unsigned int f ) const
185 {
186 assert( f < faceSize() );
187 return faceStructure_[ f ].basis_;
188 }
189
190 // basis associated to edge e
191 TestEdgeBasis *testEdgeBasis ( unsigned int e ) const
192 {
193 assert( e < edgeSize() );
194 return edgeStructure_[ e ].basis_;
195 }
196
197 const Tangent& edgeTangent ( unsigned int e ) const
198 {
199 assert( e < edgeSize() );
200 return edgeStructure_[ e ].tangent_;
201 }
202
203 const FaceTangents& faceTangents ( unsigned int f ) const
204 {
205 assert( f < faceSize() );
206 return faceStructure_[ f ].faceTangents_;
207 }
208
209 const Normal &normal ( unsigned int f ) const
210 {
211 assert( f < faceSize() );
212 return faceStructure_[ f ].normal_;
213 }
214
215 template< GeometryType::Id geometryId >
216 void build ( std::size_t order )
217 {
218 constexpr GeometryType geometry = geometryId;
219 order_ = order;
220 geometry_ = geometry;
221
222 /*
223 * The Nedelec parameter begins at 1.
224 * This is the numbering used by J.C. Nedelec himself.
225 * See "Mixed Finite Elements in \R^3" published in 1980.
226 *
227 * This construction is based on the construction of Raviart-Thomas elements.
228 * There the numbering starts at 0.
229 * Because of this we reduce the order internally by 1.
230 */
231 order--;
232
233 // if dimension == 2: order-1 on element
234 // if dimension == 3: order-2 on element
235 int requiredOrder = static_cast<int>(dimension==3);
236 testBasis_ = (order > requiredOrder ? TestBasisFactory::template create< geometry >( order-1-requiredOrder ) : nullptr);
237
238 const auto &refElement = ReferenceElements< Field, dimension >::general( type() );
239
240 numberOfFaces_ = refElement.size( 1 );
241 faceStructure_.reserve( numberOfFaces_ );
242
243 // compute the basis, tangents and normals of each face
244 for (std::size_t i=0; i<numberOfFaces_; i++)
245 {
246 FieldVector<Field,dimension> zero(0);
247 FaceTangents faceTangents;
248 faceTangents.fill(zero);
249
250 // use the first dim-1 vertices of a face to compute the tangents
251 auto vertices = refElement.subEntities(i,1,dim).begin();
252 auto vertex1 = *vertices;
253 for(int j=1; j<dim;j++)
254 {
255 auto vertex2 = vertices[j];
256
257 faceTangents[j-1] = refElement.position(vertex2,dim) - refElement.position(vertex1,dim);
258
259 // By default, edges point from the vertex with the smaller index
260 // to the vertex with the larger index.
261 if (vertex1>vertex2)
262 faceTangents[j-1] *=-1;
263
264 vertex1 = vertex2;
265 }
266
267 /* For simplices or cubes of arbitrary dimension you could just use
268 *
269 * ```
270 * GeometryType faceGeometry = Impl::getBase(geometry_);
271 * TestFaceBasis *faceBasis = ( dim == 3 && order > 0 ? TestFaceBasisFactory::template create< faceGeometry >( order-1 ) : nullptr);
272 * ```
273 *
274 * For i.e. Prisms and Pyramids in 3d this does not work because they contain squares and triangles as faces.
275 * And depending on the dynamic face index a different face geometry is needed.
276 *
277 */
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 );
280 }
281 assert( faceStructure_.size() == numberOfFaces_ );
282
283 numberOfEdges_ = refElement.size( dim-1 );
284 edgeStructure_.reserve( numberOfEdges_ );
285
286 // compute the basis and tangent of each edge
287 for (std::size_t i=0; i<numberOfEdges_; i++)
288 {
289 auto vertexIterator = refElement.subEntities(i,dim-1,dim).begin();
290 auto v0 = *vertexIterator;
291 auto v1 = *(++vertexIterator);
292
293 // By default, edges point from the vertex with the smaller index
294 // to the vertex with the larger index.
295 if (v0>v1)
296 std::swap(v0,v1);
297 auto tangent = refElement.position(v1,dim) - refElement.position(v0,dim);
298
299 TestEdgeBasis *edgeBasis = Impl::IfGeometryType< CreateEdgeBasis, 1 >::apply( refElement.type( i, dim-1 ), order );
300 edgeStructure_.emplace_back( edgeBasis, tangent );
301 }
302 assert( edgeStructure_.size() == numberOfEdges_ );
303 }
304
305 private:
306
307 // helper struct for edges
308 struct EdgeStructure
309 {
310 EdgeStructure( TestEdgeBasis *teb, const Tangent &t )
311 : basis_( teb ), tangent_( t )
312 {}
313
314 TestEdgeBasis *basis_;
316 };
317
318 template< GeometryType::Id edgeGeometryId >
319 struct CreateEdgeBasis
320 {
321 static TestEdgeBasis *apply ( std::size_t order ) { return TestEdgeBasisFactory::template create< edgeGeometryId >( order ); }
322 };
323
324 // helper struct for faces
325 struct FaceStructure
326 {
327 FaceStructure( TestFaceBasis *tfb, const Normal& normal, const FaceTangents& faceTangents )
328 : basis_( tfb ), normal_(normal), faceTangents_( faceTangents )
329 {}
330
331 TestFaceBasis *basis_;
333 const FaceTangents faceTangents_;
334 };
335
336 template< GeometryType::Id faceGeometryId >
337 struct CreateFaceBasis
338 {
339 static TestFaceBasis *apply ( std::size_t order ) { return TestFaceBasisFactory::template create< faceGeometryId >( order ); }
340 };
341
342 TestBasis *testBasis_ = nullptr;
343 std::vector< FaceStructure > faceStructure_;
344 unsigned int numberOfFaces_;
345 std::vector< EdgeStructure > edgeStructure_;
346 unsigned int numberOfEdges_;
347 GeometryType geometry_;
348 std::size_t order_;
349 };
350
351
352
353 // NedelecL2Interpolation
354 // ----------------------------
355
361 template< unsigned int dimension, class F>
363 : public InterpolationHelper< F ,dimension >
364 {
366 typedef InterpolationHelper<F,dimension> Base;
367
368 public:
369 typedef F Field;
370 typedef NedelecL2InterpolationBuilder<dimension,Field> Builder;
371 typedef typename Builder::FaceTangents FaceTangents;
372
374 : order_(0),
375 size_(0)
376 {}
377
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
382 {
383 coefficients.resize(size());
384 typename Base::template Helper<Function,Vector,true> func( function,coefficients );
385 interpolate(func);
386 }
387
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
393 {
394 matrix.resize( size(), basis.size() );
395 typename Base::template Helper<Basis,Matrix,false> func( basis,matrix );
396 interpolate(func);
397 }
398
399 std::size_t order() const
400 {
401 return order_;
402 }
403 std::size_t size() const
404 {
405 return size_;
406 }
407
408 template <GeometryType::Id geometryId>
409 void build( std::size_t order )
410 {
411 size_ = 0;
412 order_ = order;
413 builder_.template build<geometryId>(order_);
414 if (builder_.testBasis())
415 size_ += dimension*builder_.testBasis()->size();
416
417 for ( unsigned int f=0; f<builder_.faceSize(); ++f )
418 if (builder_.testFaceBasis(f))
419 size_ += (dimension-1)*builder_.testFaceBasis(f)->size();
420
421 for ( unsigned int e=0; e<builder_.edgeSize(); ++e )
422 if (builder_.testEdgeBasis(e))
423 size_ += builder_.testEdgeBasis(e)->size();
424 }
425
426 void setLocalKeys(std::vector< LocalKey > &keys) const
427 {
428 keys.resize(size());
429 unsigned int row = 0;
430 for (unsigned int e=0; e<builder_.edgeSize(); ++e)
431 {
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);
435 }
436 for (unsigned int f=0; f<builder_.faceSize(); ++f)
437 {
438 if (builder_.faceSize())
439 for (unsigned int i=0; i<builder_.testFaceBasis(f)->size()*(dimension-1); ++i,++row)
440 keys[row] = LocalKey(f,1,i);
441 }
442
443 if (builder_.testBasis())
444 for (unsigned int i=0; i<builder_.testBasis()->size()*dimension; ++i,++row)
445 keys[row] = LocalKey(0,0,i);
446 assert( row == size() );
447 }
448
449 protected:
450 template< class Func, class Container, bool type >
451 void interpolate ( typename Base::template Helper<Func,Container,type> &func ) const
452 {
453 const Dune::GeometryType geoType( builder_.topologyId(), dimension );
454
455 std::vector<Field> testBasisVal;
456
457 for (unsigned int i=0; i<size(); ++i)
458 for (unsigned int j=0; j<func.size(); ++j)
459 func.set(i,j,0);
460
461 unsigned int row = 0;
462
463 // edge dofs:
464 typedef Dune::QuadratureRule<Field, 1> EdgeQuadrature;
465 typedef Dune::QuadratureRules<Field, 1> EdgeQuadratureRules;
466
467 const auto &refElement = Dune::ReferenceElements< Field, dimension >::general( geoType );
468
469 for (unsigned int e=0; e<builder_.edgeSize(); ++e)
470 {
471 if (!builder_.testEdgeBasis(e))
472 continue;
473 testBasisVal.resize(builder_.testEdgeBasis(e)->size());
474
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 );
478
479 const unsigned int quadratureSize = edgeQuad.size();
480 for( unsigned int qi = 0; qi < quadratureSize; ++qi )
481 {
482 if (dimension>1)
483 builder_.testEdgeBasis(e)->template evaluate<0>(edgeQuad[qi].position(),testBasisVal);
484 else
485 testBasisVal[0] = 1.;
486 computeEdgeDofs(row,
487 testBasisVal,
488 func.evaluate( geometry.global( edgeQuad[qi].position() ) ),
489 builder_.edgeTangent(e),
490 edgeQuad[qi].weight(),
491 func);
492 }
493
494 row += builder_.testEdgeBasis(e)->size();
495 }
496
497 // face dofs:
498 typedef Dune::QuadratureRule<Field, dimension-1> FaceQuadrature;
499 typedef Dune::QuadratureRules<Field, dimension-1> FaceQuadratureRules;
500
501 for (unsigned int f=0; f<builder_.faceSize(); ++f)
502 {
503 if (builder_.testFaceBasis(f))
504 {
505 testBasisVal.resize(builder_.testFaceBasis(f)->size());
506
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 );
510
511 const unsigned int quadratureSize = faceQuad.size();
512 for( unsigned int qi = 0; qi < quadratureSize; ++qi )
513 {
514 if (dimension>1)
515 builder_.testFaceBasis(f)->template evaluate<0>(faceQuad[qi].position(),testBasisVal);
516 else
517 testBasisVal[0] = 1.;
518
519 computeFaceDofs( row,
520 testBasisVal,
521 func.evaluate( geometry.global( faceQuad[qi].position() ) ),
522 builder_.faceTangents(f),
523 builder_.normal(f),
524 faceQuad[qi].weight(),
525 func);
526 }
527
528 row += builder_.testFaceBasis(f)->size()*(dimension-1);
529 }
530 }
531
532 // element dofs
533 if (builder_.testBasis())
534 {
535 testBasisVal.resize(builder_.testBasis()->size());
536
539 const Quadrature &elemQuad = QuadratureRules::rule( geoType, 2*order_+1 );
540
541 const unsigned int quadratureSize = elemQuad.size();
542 for( unsigned int qi = 0; qi < quadratureSize; ++qi )
543 {
544 builder_.testBasis()->template evaluate<0>(elemQuad[qi].position(),testBasisVal);
545 computeInteriorDofs(row,
546 testBasisVal,
547 func.evaluate(elemQuad[qi].position()),
548 elemQuad[qi].weight(),
549 func );
550 }
551
552 row += builder_.testBasis()->size()*dimension;
553 }
554 assert(row==size());
555 }
556
557 private:
567 template <class MVal, class NedVal,class Matrix>
568 void computeEdgeDofs (unsigned int startRow,
569 const MVal &mVal,
570 const NedVal &nedVal,
571 const FieldVector<Field,dimension> &tangent,
572 const Field &weight,
573 Matrix &matrix) const
574 {
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)
578 {
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) );
583
584 assert( mIter == mVal.end() );
585 }
586 }
587
598 template <class MVal, class NedVal,class Matrix>
599 void computeFaceDofs (unsigned int startRow,
600 const MVal &mVal,
601 const NedVal &nedVal,
602 const FaceTangents& faceTangents,
603 const FieldVector<Field,dimension> &normal,
604 const Field &weight,
605 Matrix &matrix) const
606 {
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)
610 {
611 auto const& u=*nedIter;
612 auto const& n=normal;
613 FieldVector<Field,dimension> nedTimesNormal = { u[1]*n[2]-u[2]*n[1],
614 u[2]*n[0]-u[0]*n[2],
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)
618 {
619 for(int i=0; i<dimension-1;i++)
620 {
621 auto test = *mIter*faceTangents[i];
622 matrix.add(row+i,col, weight*(nedTimesNormal*test) );
623 }
624 row += dimension-1;
625 }
626
627 assert( mIter == mVal.end() );
628 }
629 }
630
639 template <class MVal, class NedVal,class Matrix>
640 void computeInteriorDofs (unsigned int startRow,
641 const MVal &mVal,
642 const NedVal &nedVal,
643 Field weight,
644 Matrix &matrix) const
645 {
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)
649 {
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] );
654
655 assert( mIter == mVal.end() );
656 }
657 }
658
659 public:
660 Builder builder_;
661 std::size_t order_;
662 std::size_t size_;
663 };
664
665 template < unsigned int dim, class Field >
666 struct NedelecL2InterpolationFactory
667 {
668 typedef NedelecL2InterpolationBuilder<dim,Field> Builder;
669 typedef const NedelecL2Interpolation<dim,Field> Object;
670 typedef std::size_t Key;
671 typedef typename std::remove_const<Object>::type NonConstObject;
672
673 template <GeometryType::Id geometryId>
674 static Object *create( const Key &key )
675 {
676 if ( !supports<geometryId>(key) )
677 return 0;
678 NonConstObject *interpol = new NonConstObject();
679 interpol->template build<geometryId>(key);
680 return interpol;
681 }
682
683 template <GeometryType::Id geometryId>
684 static bool supports( const Key &key )
685 {
686 GeometryType gt = geometryId;
687 return gt.isTriangle() || gt.isTetrahedron() ;
688 }
689 static void release( Object *object ) { delete object; }
690 };
691
692} // namespace Dune
693
694#endif // #ifndef DUNE_LOCALFUNCTIONS_NEDELEC_NEDELECSIMPLEX_NEDELECSIMPLEXINTERPOLATION_HH
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.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 21, 23:30, 2024)