Dune Core Modules (2.6.0)

referenceelementimplementation.hh
1// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2// vi: set et ts=4 sw=2 sts=2:
3#ifndef DUNE_GEOMETRY_REFERENCEELEMENTIMPLEMENTATION_HH
4#define DUNE_GEOMETRY_REFERENCEELEMENTIMPLEMENTATION_HH
5
6#include <cassert>
7
8#include <algorithm>
9#include <limits>
10#include <tuple>
11#include <utility>
12#include <vector>
13#include <array>
14
17#include <dune/common/hybridutilities.hh>
19#include <dune/common/unused.hh>
20
21#include <dune/geometry/referenceelement.hh>
23#include <dune/geometry/type.hh>
24
25namespace Dune
26{
27
28 namespace Geo
29 {
30
31#ifndef DOXYGEN
32
33 // Internal Forward Declarations
34 // -----------------------------
35
36 namespace Impl
37 {
38 template< class ctype, int dim >
39 class ReferenceElementContainer;
40 }
41
42 template< class ctype, int dim >
43 struct ReferenceElements;
44
45
46
47 namespace Impl
48 {
49
50 using Dune::Impl::isPrism;
51 using Dune::Impl::isPyramid;
52 using Dune::Impl::baseTopologyId;
53 using Dune::Impl::prismConstruction;
54 using Dune::Impl::pyramidConstruction;
55 using Dune::Impl::numTopologies;
56
58 unsigned int size ( unsigned int topologyId, int dim, int codim );
59
60
61
69 unsigned int subTopologyId ( unsigned int topologyId, int dim, int codim, unsigned int i );
70
71
72
73 // subTopologyNumbering
74 // --------------------
75
76 void subTopologyNumbering ( unsigned int topologyId, int dim, int codim, unsigned int i, int subcodim,
77 unsigned int *beginOut, unsigned int *endOut );
78
79
80
81
82 // checkInside
83 // -----------
84
85 template< class ct, int cdim >
86 inline bool
87 checkInside ( unsigned int topologyId, int dim, const FieldVector< ct, cdim > &x, ct tolerance, ct factor = ct( 1 ) )
88 {
89 assert( (dim >= 0) && (dim <= cdim) );
90 assert( topologyId < numTopologies( dim ) );
91
92 if( dim > 0 )
93 {
94 const ct baseFactor = (isPrism( topologyId, dim ) ? factor : factor - x[ dim-1 ]);
95 if( (x[ dim-1 ] > -tolerance) && (factor - x[ dim-1 ] > -tolerance) )
96 return checkInside< ct, cdim >( baseTopologyId( topologyId, dim ), dim-1, x, tolerance, baseFactor );
97 else
98 return false;
99 }
100 else
101 return true;
102 }
103
104
105
106 // referenceCorners
107 // ----------------
108
109 template< class ct, int cdim >
110 inline unsigned int
111 referenceCorners ( unsigned int topologyId, int dim, FieldVector< ct, cdim > *corners )
112 {
113 assert( (dim >= 0) && (dim <= cdim) );
114 assert( topologyId < numTopologies( dim ) );
115
116 if( dim > 0 )
117 {
118 const unsigned int nBaseCorners
119 = referenceCorners( baseTopologyId( topologyId, dim ), dim-1, corners );
120 assert( nBaseCorners == size( baseTopologyId( topologyId, dim ), dim-1, dim-1 ) );
121 if( isPrism( topologyId, dim ) )
122 {
123 std::copy( corners, corners + nBaseCorners, corners + nBaseCorners );
124 for( unsigned int i = 0; i < nBaseCorners; ++i )
125 corners[ i+nBaseCorners ][ dim-1 ] = ct( 1 );
126 return 2*nBaseCorners;
127 }
128 else
129 {
130 corners[ nBaseCorners ] = FieldVector< ct, cdim >( ct( 0 ) );
131 corners[ nBaseCorners ][ dim-1 ] = ct( 1 );
132 return nBaseCorners+1;
133 }
134 }
135 else
136 {
137 *corners = FieldVector< ct, cdim >( ct( 0 ) );
138 return 1;
139 }
140 }
141
142
143
144 // referenceVolume
145 // ---------------
146
147 unsigned long referenceVolumeInverse ( unsigned int topologyId, int dim );
148
149 template< class ct >
150 inline ct referenceVolume ( unsigned int topologyId, int dim )
151 {
152 return ct( 1 ) / ct( referenceVolumeInverse( topologyId, dim ) );
153 }
154
155
156
157 // referenceOrigins
158 // ----------------
159
160 template< class ct, int cdim >
161 inline unsigned int
162 referenceOrigins ( unsigned int topologyId, int dim, int codim, FieldVector< ct, cdim > *origins )
163 {
164 assert( (dim >= 0) && (dim <= cdim) );
165 assert( topologyId < numTopologies( dim ) );
166 assert( (codim >= 0) && (codim <= dim) );
167
168 if( codim > 0 )
169 {
170 const unsigned int baseId = baseTopologyId( topologyId, dim );
171 if( isPrism( topologyId, dim ) )
172 {
173 const unsigned int n = (codim < dim ? referenceOrigins( baseId, dim-1, codim, origins ) : 0);
174 const unsigned int m = referenceOrigins( baseId, dim-1, codim-1, origins+n );
175 for( unsigned int i = 0; i < m; ++i )
176 {
177 origins[ n+m+i ] = origins[ n+i ];
178 origins[ n+m+i ][ dim-1 ] = ct( 1 );
179 }
180 return n+2*m;
181 }
182 else
183 {
184 const unsigned int m = referenceOrigins( baseId, dim-1, codim-1, origins );
185 if( codim == dim )
186 {
187 origins[ m ] = FieldVector< ct, cdim >( ct( 0 ) );
188 origins[ m ][ dim-1 ] = ct( 1 );
189 return m+1;
190 }
191 else
192 return m+referenceOrigins( baseId, dim-1, codim, origins+m );
193 }
194 }
195 else
196 {
197 origins[ 0 ] = FieldVector< ct, cdim >( ct( 0 ) );
198 return 1;
199 }
200 }
201
202
203
204 // referenceEmbeddings
205 // -------------------
206
207 template< class ct, int cdim, int mydim >
208 inline unsigned int
209 referenceEmbeddings ( unsigned int topologyId, int dim, int codim,
210 FieldVector< ct, cdim > *origins,
211 FieldMatrix< ct, mydim, cdim > *jacobianTransposeds )
212 {
213 assert( (0 <= codim) && (codim <= dim) && (dim <= cdim) );
214 assert( (dim - codim <= mydim) && (mydim <= cdim) );
215 assert( topologyId < numTopologies( dim ) );
216
217 if( codim > 0 )
218 {
219 const unsigned int baseId = baseTopologyId( topologyId, dim );
220 if( isPrism( topologyId, dim ) )
221 {
222 const unsigned int n = (codim < dim ? referenceEmbeddings( baseId, dim-1, codim, origins, jacobianTransposeds ) : 0);
223 for( unsigned int i = 0; i < n; ++i )
224 jacobianTransposeds[ i ][ dim-codim-1 ][ dim-1 ] = ct( 1 );
225
226 const unsigned int m = referenceEmbeddings( baseId, dim-1, codim-1, origins+n, jacobianTransposeds+n );
227 std::copy( origins+n, origins+n+m, origins+n+m );
228 std::copy( jacobianTransposeds+n, jacobianTransposeds+n+m, jacobianTransposeds+n+m );
229 for( unsigned int i = 0; i < m; ++i )
230 origins[ n+m+i ][ dim-1 ] = ct( 1 );
231
232 return n+2*m;
233 }
234 else
235 {
236 const unsigned int m = referenceEmbeddings( baseId, dim-1, codim-1, origins, jacobianTransposeds );
237 if( codim == dim )
238 {
239 origins[ m ] = FieldVector< ct, cdim >( ct( 0 ) );
240 origins[ m ][ dim-1 ] = ct( 1 );
241 jacobianTransposeds[ m ] = FieldMatrix< ct, mydim, cdim >( ct( 0 ) );
242 return m+1;
243 }
244 else
245 {
246 const unsigned int n = referenceEmbeddings( baseId, dim-1, codim, origins+m, jacobianTransposeds+m );
247 for( unsigned int i = 0; i < n; ++i )
248 {
249 for( int k = 0; k < dim-1; ++k )
250 jacobianTransposeds[ m+i ][ dim-codim-1 ][ k ] = -origins[ m+i ][ k ];
251 jacobianTransposeds[ m+i ][ dim-codim-1 ][ dim-1 ] = ct( 1 );
252 }
253 return m+n;
254 }
255 }
256 }
257 else
258 {
259 origins[ 0 ] = FieldVector< ct, cdim >( ct( 0 ) );
260 jacobianTransposeds[ 0 ] = FieldMatrix< ct, mydim, cdim >( ct( 0 ) );
261 for( int k = 0; k < dim; ++k )
262 jacobianTransposeds[ 0 ][ k ][ k ] = ct( 1 );
263 return 1;
264 }
265 }
266
267
268
269 // referenceIntegrationOuterNormals
270 // --------------------------------
271
272 template< class ct, int cdim >
273 inline unsigned int
274 referenceIntegrationOuterNormals ( unsigned int topologyId, int dim,
275 const FieldVector< ct, cdim > *origins,
276 FieldVector< ct, cdim > *normals )
277 {
278 assert( (dim > 0) && (dim <= cdim) );
279 assert( topologyId < numTopologies( dim ) );
280
281 if( dim > 1 )
282 {
283 const unsigned int baseId = baseTopologyId( topologyId, dim );
284 if( isPrism( topologyId, dim ) )
285 {
286 const unsigned int numBaseFaces
287 = referenceIntegrationOuterNormals( baseId, dim-1, origins, normals );
288
289 for( unsigned int i = 0; i < 2; ++i )
290 {
291 normals[ numBaseFaces+i ] = FieldVector< ct, cdim >( ct( 0 ) );
292 normals[ numBaseFaces+i ][ dim-1 ] = ct( 2*int( i )-1 );
293 }
294
295 return numBaseFaces+2;
296 }
297 else
298 {
299 normals[ 0 ] = FieldVector< ct, cdim >( ct( 0 ) );
300 normals[ 0 ][ dim-1 ] = ct( -1 );
301
302 const unsigned int numBaseFaces
303 = referenceIntegrationOuterNormals( baseId, dim-1, origins+1, normals+1 );
304 for( unsigned int i = 1; i <= numBaseFaces; ++i )
305 normals[ i ][ dim-1 ] = normals[ i ]*origins[ i ];
306
307 return numBaseFaces+1;
308 }
309 }
310 else
311 {
312 for( unsigned int i = 0; i < 2; ++i )
313 {
314 normals[ i ] = FieldVector< ct, cdim >( ct( 0 ) );
315 normals[ i ][ 0 ] = ct( 2*int( i )-1 );
316 }
317
318 return 2;
319 }
320 }
321
322 template< class ct, int cdim >
323 inline unsigned int
324 referenceIntegrationOuterNormals ( unsigned int topologyId, int dim,
325 FieldVector< ct, cdim > *normals )
326 {
327 assert( (dim > 0) && (dim <= cdim) );
328
329 FieldVector< ct, cdim > *origins
330 = new FieldVector< ct, cdim >[ size( topologyId, dim, 1 ) ];
331 referenceOrigins( topologyId, dim, 1, origins );
332
333 const unsigned int numFaces
334 = referenceIntegrationOuterNormals( topologyId, dim, origins, normals );
335 assert( numFaces == size( topologyId, dim, 1 ) );
336
337 delete[] origins;
338
339 return numFaces;
340 }
341
342 } // namespace Impl
343
344
345
346 // ReferenceElement
347 // ----------------
348
367 template< class ctype_, int dim >
368 class ReferenceElementImplementation
369 {
370
371 public:
372
374 using ctype = ctype_;
375
377 using CoordinateField = ctype;
378
380 using Coordinate = Dune::FieldVector<ctype,dim>;
381
383 static constexpr int dimension = dim;
384
385 private:
386
387 friend class Impl::ReferenceElementContainer< ctype, dim >;
388
389 struct SubEntityInfo;
390
391 template< int codim > struct CreateGeometries;
392
393 public:
395 template< int codim >
396 struct Codim
397 {
399 typedef AffineGeometry< ctype, dim-codim, dim > Geometry;
400 };
401
402 // ReferenceElement cannot be copied.
403 ReferenceElementImplementation ( const ReferenceElementImplementation& ) = delete;
404
405 // ReferenceElementImplementation cannot be copied.
406 ReferenceElementImplementation& operator= ( const ReferenceElementImplementation& ) = delete;
407
408 // ReferenceElementImplementation is default-constructible (required for storage in std::array)
409 ReferenceElementImplementation () = default;
410
415 int size ( int c ) const
416 {
417 assert( (c >= 0) && (c <= dim) );
418 return info_[ c ].size();
419 }
420
432 int size ( int i, int c, int cc ) const
433 {
434 assert( (i >= 0) && (i < size( c )) );
435 return info_[ c ][ i ].size( cc );
436 }
437
451 int subEntity ( int i, int c, int ii, int cc ) const
452 {
453 assert( (i >= 0) && (i < size( c )) );
454 return info_[ c ][ i ].number( ii, cc );
455 }
456
465 const GeometryType &type ( int i, int c ) const
466 {
467 assert( (i >= 0) && (i < size( c )) );
468 return info_[ c ][ i ].type();
469 }
470
472 const GeometryType &type () const { return type( 0, 0 ); }
473
483 const Coordinate &position( int i, int c ) const
484 {
485 assert( (c >= 0) && (c <= dim) );
486 return baryCenters_[ c ][ i ];
487 }
488
496 bool checkInside ( const Coordinate &local ) const
497 {
498 const ctype tolerance = ctype( 64 ) * std::numeric_limits< ctype >::epsilon();
499 return Impl::template checkInside< ctype, dim >( type().id(), dim, local, tolerance );
500 }
501
513 template< int codim >
514 typename Codim< codim >::Geometry geometry ( int i ) const
515 {
516 return std::get< codim >( geometries_ )[ i ];
517 }
518
520 ctype volume () const
521 {
522 return volume_;
523 }
524
532 const Coordinate &integrationOuterNormal ( int face ) const
533 {
534 assert( (face >= 0) && (face < int( integrationNormals_.size() )) );
535 return integrationNormals_[ face ];
536 }
537
538 private:
539 void initialize ( unsigned int topologyId )
540 {
541 assert( topologyId < Impl::numTopologies( dim ) );
542
543 // set up subentities
544 for( int codim = 0; codim <= dim; ++codim )
545 {
546 const unsigned int size = Impl::size( topologyId, dim, codim );
547 info_[ codim ].resize( size );
548 for( unsigned int i = 0; i < size; ++i )
549 info_[ codim ][ i ].initialize( topologyId, codim, i );
550 }
551
552 // compute corners
553 const unsigned int numVertices = size( dim );
554 baryCenters_[ dim ].resize( numVertices );
555 Impl::referenceCorners( topologyId, dim, &(baryCenters_[ dim ][ 0 ]) );
556
557 // compute barycenters
558 for( int codim = 0; codim < dim; ++codim )
559 {
560 baryCenters_[ codim ].resize( size(codim) );
561 for( int i = 0; i < size( codim ); ++i )
562 {
563 baryCenters_[ codim ][ i ] = Coordinate( ctype( 0 ) );
564 const unsigned int numCorners = size( i, codim, dim );
565 for( unsigned int j = 0; j < numCorners; ++j )
566 baryCenters_[ codim ][ i ] += baryCenters_[ dim ][ subEntity( i, codim, j, dim ) ];
567 baryCenters_[ codim ][ i ] *= ctype( 1 ) / ctype( numCorners );
568 }
569 }
570
571 // compute reference element volume
572 volume_ = Impl::template referenceVolume< ctype >( topologyId, dim );
573
574 // compute integration outer normals
575 if( dim > 0 )
576 {
577 integrationNormals_.resize( size( 1 ) );
578 Impl::referenceIntegrationOuterNormals( topologyId, dim, &(integrationNormals_[ 0 ]) );
579 }
580
581 // set up geometries
582 Hybrid::forEach( std::make_index_sequence< dim+1 >{}, [ & ]( auto i ){ CreateGeometries< i >::apply( *this, geometries_ ); } );
583 }
584
585 template< int... codim >
586 static std::tuple< std::vector< typename Codim< codim >::Geometry >... >
587 makeGeometryTable ( std::integer_sequence< int, codim... > );
588
590 typedef decltype( makeGeometryTable( std::make_integer_sequence< int, dim+1 >() ) ) GeometryTable;
591
593 ctype volume_;
594
595 std::vector< Coordinate > baryCenters_[ dim+1 ];
596 std::vector< Coordinate > integrationNormals_;
597
599 GeometryTable geometries_;
600
601 std::vector< SubEntityInfo > info_[ dim+1 ];
602 };
603
605 template< class ctype, int dim >
606 struct ReferenceElementImplementation< ctype, dim >::SubEntityInfo
607 {
608 SubEntityInfo ()
609 : numbering_( nullptr )
610 {
611 std::fill( offset_.begin(), offset_.end(), 0 );
612 }
613
614 SubEntityInfo ( const SubEntityInfo &other )
615 : offset_( other.offset_ ),
616 type_( other.type_ )
617 {
618 numbering_ = allocate();
619 std::copy( other.numbering_, other.numbering_ + capacity(), numbering_ );
620 }
621
622 ~SubEntityInfo () { deallocate( numbering_ ); }
623
624 const SubEntityInfo &operator= ( const SubEntityInfo &other )
625 {
626 type_ = other.type_;
627 offset_ = other.offset_;
628
629 deallocate( numbering_ );
630 numbering_ = allocate();
631 std::copy( other.numbering_, other.numbering_ + capacity(), numbering_ );
632
633 return *this;
634 }
635
636 int size ( int cc ) const
637 {
638 assert( (cc >= codim()) && (cc <= dim) );
639 return (offset_[ cc+1 ] - offset_[ cc ]);
640 }
641
642 int number ( int ii, int cc ) const
643 {
644 assert( (ii >= 0) && (ii < size( cc )) );
645 return numbering_[ offset_[ cc ] + ii ];
646 }
647
648 const GeometryType &type () const { return type_; }
649
650 void initialize ( unsigned int topologyId, int codim, unsigned int i )
651 {
652 const unsigned int subId = Impl::subTopologyId( topologyId, dim, codim, i );
653 type_ = GeometryType( subId, dim-codim );
654
655 // compute offsets
656 for( int cc = 0; cc <= codim; ++cc )
657 offset_[ cc ] = 0;
658 for( int cc = codim; cc <= dim; ++cc )
659 offset_[ cc+1 ] = offset_[ cc ] + Impl::size( subId, dim-codim, cc-codim );
660
661 // compute subnumbering
662 deallocate( numbering_ );
663 numbering_ = allocate();
664 for( int cc = codim; cc <= dim; ++cc )
665 Impl::subTopologyNumbering( topologyId, dim, codim, i, cc-codim, numbering_+offset_[ cc ], numbering_+offset_[ cc+1 ] );
666 }
667
668 protected:
669 int codim () const { return dim - type().dim(); }
670
671 unsigned int *allocate () { return (capacity() != 0 ? new unsigned int[ capacity() ] : nullptr); }
672 void deallocate ( unsigned int *ptr ) { delete[] ptr; }
673 unsigned int capacity () const { return offset_[ dim+1 ]; }
674
675 private:
676 unsigned int *numbering_;
677 std::array< unsigned int, dim+2 > offset_;
678 GeometryType type_;
679 };
680
681
682 template< class ctype, int dim >
683 template< int codim >
684 struct ReferenceElementImplementation< ctype, dim >::CreateGeometries
685 {
686 template< int cc >
687 static typename ReferenceElements< ctype, dim-cc >::ReferenceElement
688 subRefElement( const ReferenceElementImplementation< ctype, dim > &refElement, int i, std::integral_constant< int, cc > )
689 {
690 return ReferenceElements< ctype, dim-cc >::general( refElement.type( i, cc ) );
691 }
692
694 subRefElement( const ReferenceElementImplementation< ctype, dim > &refElement, int i, std::integral_constant< int, 0 > )
695 {
697 return refElement;
698 }
699
700 static void apply ( const ReferenceElementImplementation< ctype, dim > &refElement, GeometryTable &geometries )
701 {
702 const int size = refElement.size( codim );
703 std::vector< FieldVector< ctype, dim > > origins( size );
704 std::vector< FieldMatrix< ctype, dim - codim, dim > > jacobianTransposeds( size );
705 Impl::referenceEmbeddings( refElement.type().id(), dim, codim, &(origins[ 0 ]), &(jacobianTransposeds[ 0 ]) );
706
707 std::get< codim >( geometries ).reserve( size );
708 for( int i = 0; i < size; ++i )
709 {
710 typename Codim< codim >::Geometry geometry( subRefElement( refElement, i, std::integral_constant< int, codim >() ), origins[ i ], jacobianTransposeds[ i ] );
711 std::get< codim >( geometries ).push_back( geometry );
712 }
713 }
714 };
715
716#endif // DOXYGEN
717
718 } // namespace Geo
719
720} // namespace Dune
721
722#endif // #ifndef DUNE_GEOMETRY_REFERENCEELEMENTIMPLEMENTATION_HH
An implementation of the Geometry interface for affine geometries.
GeometryType
Type representing VTK's entity geometry types.
Definition: common.hh:178
Implements a matrix constructed from a given type representing a field and compile-time given number ...
Implements a vector constructed from a given type representing a field and a compile-time given size.
decltype(auto) apply(F &&f, ArgTuple &&args)
Apply function with arguments given as tuple.
Definition: apply.hh:58
#define DUNE_UNUSED_PARAMETER(parm)
A macro to mark intentionally unused function parameters with.
Definition: unused.hh:25
unspecified-type ReferenceElement
Returns the type of reference element for the argument type T.
Definition: referenceelements.hh:495
constexpr void forEach(Range &&range, F &&f)
Range based for loop.
Definition: hybridutilities.hh:308
Dune namespace.
Definition: alignedallocator.hh:10
typename Container::ReferenceElement ReferenceElement
The reference element type.
Definition: referenceelements.hh:186
static const ReferenceElement & general(const GeometryType &type)
get general reference elements
Definition: referenceelements.hh:196
A unique label for each type of element that can occur in a grid.
Traits for type conversions and type information.
Definition of the DUNE_UNUSED macro for the case that config.h is not available.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 13, 23:29, 2024)