3#ifndef DUNE_GEOMETRY_REFERENCEELEMENTS_HH
4#define DUNE_GEOMETRY_REFERENCEELEMENTS_HH
31 template<
class ctype,
int dim >
32 class ReferenceElementContainer;
34 template<
class ctype,
int dim >
35 struct ReferenceElements;
43 unsigned int size (
unsigned int topologyId,
int dim,
int codim );
54 unsigned int subTopologyId (
unsigned int topologyId,
int dim,
int codim,
unsigned int i );
61 void subTopologyNumbering (
unsigned int topologyId,
int dim,
int codim,
unsigned int i,
int subcodim,
62 unsigned int *beginOut,
unsigned int *endOut );
70 template<
class ct,
int cdim >
72 checkInside (
unsigned int topologyId,
int dim,
const FieldVector< ct, cdim > &x, ct tolerance, ct factor = ct( 1 ) )
74 assert( (dim >= 0) && (dim <= cdim) );
75 assert( topologyId < numTopologies( dim ) );
79 const ct baseFactor = (isPrism( topologyId, dim ) ? factor : factor - x[ dim-1 ]);
80 if( (x[ dim-1 ] > -tolerance) && (factor - x[ dim-1 ] > -tolerance) )
81 return checkInside< ct, cdim >( baseTopologyId( topologyId, dim ), dim-1, x, tolerance, baseFactor );
94 template<
class ct,
int cdim >
96 referenceCorners (
unsigned int topologyId,
int dim, FieldVector< ct, cdim > *corners )
98 assert( (dim >= 0) && (dim <= cdim) );
99 assert( topologyId < numTopologies( dim ) );
103 const unsigned int nBaseCorners
104 = referenceCorners( baseTopologyId( topologyId, dim ), dim-1, corners );
105 assert( nBaseCorners == size( baseTopologyId( topologyId, dim ), dim-1, dim-1 ) );
106 if( isPrism( topologyId, dim ) )
108 std::copy( corners, corners + nBaseCorners, corners + nBaseCorners );
109 for(
unsigned int i = 0; i < nBaseCorners; ++i )
110 corners[ i+nBaseCorners ][ dim-1 ] = ct( 1 );
111 return 2*nBaseCorners;
115 corners[ nBaseCorners ] = FieldVector< ct, cdim >( ct( 0 ) );
116 corners[ nBaseCorners ][ dim-1 ] = ct( 1 );
117 return nBaseCorners+1;
122 *corners = FieldVector< ct, cdim >( ct( 0 ) );
132 unsigned long referenceVolumeInverse (
unsigned int topologyId,
int dim );
135 inline ct referenceVolume (
unsigned int topologyId,
int dim )
137 return ct( 1 ) / ct( referenceVolumeInverse( topologyId, dim ) );
145 template<
class ct,
int cdim >
147 referenceOrigins (
unsigned int topologyId,
int dim,
int codim, FieldVector< ct, cdim > *origins )
149 assert( (dim >= 0) && (dim <= cdim) );
150 assert( topologyId < numTopologies( dim ) );
151 assert( (codim >= 0) && (codim <= dim) );
155 const unsigned int baseId = baseTopologyId( topologyId, dim );
156 if( isPrism( topologyId, dim ) )
158 const unsigned int n = (codim < dim ? referenceOrigins( baseId, dim-1, codim, origins ) : 0);
159 const unsigned int m = referenceOrigins( baseId, dim-1, codim-1, origins+n );
160 for(
unsigned int i = 0; i < m; ++i )
162 origins[ n+m+i ] = origins[ n+i ];
163 origins[ n+m+i ][ dim-1 ] = ct( 1 );
169 const unsigned int m = referenceOrigins( baseId, dim-1, codim-1, origins );
172 origins[ m ] = FieldVector< ct, cdim >( ct( 0 ) );
173 origins[ m ][ dim-1 ] = ct( 1 );
177 return m+referenceOrigins( baseId, dim-1, codim, origins+m );
182 origins[ 0 ] = FieldVector< ct, cdim >( ct( 0 ) );
192 template<
class ct,
int cdim,
int mydim >
194 referenceEmbeddings (
unsigned int topologyId,
int dim,
int codim,
195 FieldVector< ct, cdim > *origins,
196 FieldMatrix< ct, mydim, cdim > *jacobianTransposeds )
198 assert( (0 <= codim) && (codim <= dim) && (dim <= cdim) );
199 assert( (dim - codim <= mydim) && (mydim <= cdim) );
200 assert( topologyId < numTopologies( dim ) );
204 const unsigned int baseId = baseTopologyId( topologyId, dim );
205 if( isPrism( topologyId, dim ) )
207 const unsigned int n = (codim < dim ? referenceEmbeddings( baseId, dim-1, codim, origins, jacobianTransposeds ) : 0);
208 for(
unsigned int i = 0; i < n; ++i )
209 jacobianTransposeds[ i ][ dim-codim-1 ][ dim-1 ] = ct( 1 );
211 const unsigned int m = referenceEmbeddings( baseId, dim-1, codim-1, origins+n, jacobianTransposeds+n );
212 std::copy( origins+n, origins+n+m, origins+n+m );
213 std::copy( jacobianTransposeds+n, jacobianTransposeds+n+m, jacobianTransposeds+n+m );
214 for(
unsigned int i = 0; i < m; ++i )
215 origins[ n+m+i ][ dim-1 ] = ct( 1 );
221 const unsigned int m = referenceEmbeddings( baseId, dim-1, codim-1, origins, jacobianTransposeds );
224 origins[ m ] = FieldVector< ct, cdim >( ct( 0 ) );
225 origins[ m ][ dim-1 ] = ct( 1 );
226 jacobianTransposeds[ m ] = FieldMatrix< ct, mydim, cdim >( ct( 0 ) );
231 const unsigned int n = referenceEmbeddings( baseId, dim-1, codim, origins+m, jacobianTransposeds+m );
232 for(
unsigned int i = 0; i < n; ++i )
234 for(
int k = 0; k < dim-1; ++k )
235 jacobianTransposeds[ m+i ][ dim-codim-1 ][ k ] = -origins[ m+i ][ k ];
236 jacobianTransposeds[ m+i ][ dim-codim-1 ][ dim-1 ] = ct( 1 );
244 origins[ 0 ] = FieldVector< ct, cdim >( ct( 0 ) );
245 jacobianTransposeds[ 0 ] = FieldMatrix< ct, mydim, cdim >( ct( 0 ) );
246 for(
int k = 0; k < dim; ++k )
247 jacobianTransposeds[ 0 ][ k ][ k ] = ct( 1 );
257 template<
class ct,
int cdim >
259 referenceIntegrationOuterNormals (
unsigned int topologyId,
int dim,
260 const FieldVector< ct, cdim > *origins,
261 FieldVector< ct, cdim > *normals )
263 assert( (dim > 0) && (dim <= cdim) );
264 assert( topologyId < numTopologies( dim ) );
268 const unsigned int baseId = baseTopologyId( topologyId, dim );
269 if( isPrism( topologyId, dim ) )
271 const unsigned int numBaseFaces
272 = referenceIntegrationOuterNormals( baseId, dim-1, origins, normals );
274 for(
unsigned int i = 0; i < 2; ++i )
276 normals[ numBaseFaces+i ] = FieldVector< ct, cdim >( ct( 0 ) );
277 normals[ numBaseFaces+i ][ dim-1 ] = ct( 2*
int( i )-1 );
280 return numBaseFaces+2;
284 normals[ 0 ] = FieldVector< ct, cdim >( ct( 0 ) );
285 normals[ 0 ][ dim-1 ] = ct( -1 );
287 const unsigned int numBaseFaces
288 = referenceIntegrationOuterNormals( baseId, dim-1, origins+1, normals+1 );
289 for(
unsigned int i = 1; i <= numBaseFaces; ++i )
290 normals[ i ][ dim-1 ] = normals[ i ]*origins[ i ];
292 return numBaseFaces+1;
297 for(
unsigned int i = 0; i < 2; ++i )
299 normals[ i ] = FieldVector< ct, cdim >( ct( 0 ) );
300 normals[ i ][ 0 ] = ct( 2*
int( i )-1 );
307 template<
class ct,
int cdim >
309 referenceIntegrationOuterNormals (
unsigned int topologyId,
int dim,
310 FieldVector< ct, cdim > *normals )
312 assert( (dim > 0) && (dim <= cdim) );
314 FieldVector< ct, cdim > *origins
315 =
new FieldVector< ct, cdim >[ size( topologyId, dim, 1 ) ];
316 referenceOrigins( topologyId, dim, 1, origins );
318 const unsigned int numFaces
319 = referenceIntegrationOuterNormals( topologyId, dim, origins, normals );
320 assert( numFaces == size( topologyId, dim, 1 ) );
352 template<
class ctype,
int dim >
357 friend class ReferenceElementContainer< ctype, dim >;
366 template<
int codim >
struct CreateGeometries;
370 template<
int codim >
383 assert( (c >= 0) && (c <= dim) );
384 return info_[ c ].size();
398 int size (
int i,
int c,
int cc )
const
400 assert( (i >= 0) && (i <
size( c )) );
401 return info_[ c ][ i ].size( cc );
419 assert( (i >= 0) && (i <
size( c )) );
420 return info_[ c ][ i ].number( ii, cc );
433 assert( (i >= 0) && (i <
size( c )) );
434 return info_[ c ][ i ].type();
451 assert( (c >= 0) && (c <= dim) );
452 return baryCenters_[ c ][ i ];
464 const ctype tolerance = ctype( 64 ) * std::numeric_limits< ctype >::epsilon();
465 return Impl::template checkInside< ctype, dim >(
type().
id(), dim, local, tolerance );
479 template<
int codim >
482 return std::get< codim >( geometries_ )[ i ];
500 assert( (face >= 0) && (face <
int( integrationNormals_.size() )) );
501 return integrationNormals_[ face ];
505 void initialize (
unsigned int topologyId )
507 assert( topologyId < Impl::numTopologies( dim ) );
510 for(
int codim = 0; codim <= dim; ++codim )
512 const unsigned int size = Impl::size( topologyId, dim, codim );
513 info_[ codim ].resize(
size );
514 for(
unsigned int i = 0; i <
size; ++i )
515 info_[ codim ][ i ].initialize( topologyId, codim, i );
519 const unsigned int numVertices =
size( dim );
520 baryCenters_[ dim ].resize( numVertices );
521 Impl::referenceCorners( topologyId, dim, &(baryCenters_[ dim ][ 0 ]) );
524 for(
int codim = 0; codim < dim; ++codim )
526 baryCenters_[ codim ].resize(
size(codim) );
527 for(
int i = 0; i <
size( codim ); ++i )
530 const unsigned int numCorners =
size( i, codim, dim );
531 for(
unsigned int j = 0; j < numCorners; ++j )
532 baryCenters_[ codim ][ i ] += baryCenters_[ dim ][
subEntity( i, codim, j, dim ) ];
533 baryCenters_[ codim ][ i ] *= ctype( 1 ) / ctype( numCorners );
538 volume_ = Impl::template referenceVolume< ctype >( topologyId, dim );
543 integrationNormals_.resize(
size( 1 ) );
544 Impl::referenceIntegrationOuterNormals( topologyId, dim, &(integrationNormals_[ 0 ]) );
551 template<
int... codim >
552 static std::tuple< std::vector< typename Codim< codim >::Geometry >... >
553 makeGeometryTable ( std::integer_sequence< int, codim... > );
556 typedef decltype( makeGeometryTable( std::make_integer_sequence< int, dim+1 >() ) ) GeometryTable;
561 std::vector< FieldVector< ctype, dim > > baryCenters_[ dim+1 ];
562 std::vector< FieldVector< ctype, dim > > integrationNormals_;
565 GeometryTable geometries_;
567 std::vector< SubEntityInfo > info_[ dim+1 ];
571 template<
class ctype,
int dim >
575 : numbering_(
nullptr )
577 std::fill( offset_.begin(), offset_.end(), 0 );
581 : offset_( other.offset_ ),
584 numbering_ = allocate();
585 std::copy( other.numbering_, other.numbering_ + capacity(), numbering_ );
593 offset_ = other.offset_;
595 deallocate( numbering_ );
596 numbering_ = allocate();
597 std::copy( other.numbering_, other.numbering_ + capacity(), numbering_ );
602 int size (
int cc )
const
604 assert( (cc >= codim()) && (cc <= dim) );
605 return (offset_[ cc+1 ] - offset_[ cc ]);
608 int number (
int ii,
int cc )
const
610 assert( (ii >= 0) && (ii <
size( cc )) );
611 return numbering_[ offset_[ cc ] + ii ];
616 void initialize (
unsigned int topologyId,
int codim,
unsigned int i )
618 const unsigned int subId = Impl::subTopologyId( topologyId, dim, codim, i );
622 for(
int cc = 0; cc <= codim; ++cc )
624 for(
int cc = codim; cc <= dim; ++cc )
625 offset_[ cc+1 ] = offset_[ cc ] + Impl::size( subId, dim-codim, cc-codim );
628 deallocate( numbering_ );
629 numbering_ = allocate();
630 for(
int cc = codim; cc <= dim; ++cc )
631 Impl::subTopologyNumbering( topologyId, dim, codim, i, cc-codim, numbering_+offset_[ cc ], numbering_+offset_[ cc+1 ] );
635 int codim ()
const {
return dim -
type().
dim(); }
637 unsigned int *allocate () {
return (capacity() != 0 ?
new unsigned int[ capacity() ] :
nullptr); }
638 void deallocate (
unsigned int *ptr ) {
delete[] ptr; }
639 unsigned int capacity ()
const {
return offset_[ dim+1 ]; }
642 unsigned int *numbering_;
643 std::array< unsigned int, dim+2 > offset_;
648 template<
class ctype,
int dim >
649 template<
int codim >
666 static void apply (
const ReferenceElement< ctype, dim > &refElement, GeometryTable &geometries )
668 const int size = refElement.size( codim );
669 std::vector< FieldVector< ctype, dim > > origins(
size );
670 std::vector< FieldMatrix< ctype, dim - codim, dim > > jacobianTransposeds(
size );
671 Impl::referenceEmbeddings( refElement.type().id(), dim, codim, &(origins[ 0 ]), &(jacobianTransposeds[ 0 ]) );
673 std::get< codim >( geometries ).reserve(
size );
674 for(
int i = 0; i <
size; ++i )
676 typename Codim< codim >::Geometry geometry( subRefElement( refElement, i, std::integral_constant< int, codim >() ), origins[ i ], jacobianTransposeds[ i ] );
677 std::get< codim >( geometries ).push_back(
geometry );
687 template<
class ctype,
int dim >
688 class ReferenceElementContainer
690 static const unsigned int numTopologies = (1u << dim);
693 typedef ReferenceElement< ctype, dim > value_type;
694 typedef const value_type *const_iterator;
696 ReferenceElementContainer ()
698 for(
unsigned int topologyId = 0; topologyId < numTopologies; ++topologyId )
699 values_[ topologyId ].initialize( topologyId );
702 const value_type &operator() (
const GeometryType &type )
const
704 assert( type.dim() == dim );
705 return values_[ type.id() ];
708 const value_type &simplex ()
const
710 return values_[ Impl::SimplexTopology< dim >::type::id ];
713 const value_type &cube ()
const
715 return values_[ Impl::CubeTopology< dim >::type::id ];
718 const value_type &pyramid ()
const
720 return values_[ Impl::PyramidTopology< dim >::type::id ];
723 const value_type &prism ()
const
725 return values_[ Impl::PrismTopology< dim >::type::id ];
728 const_iterator begin ()
const {
return values_; }
729 const_iterator end ()
const {
return values_ + numTopologies; }
732 value_type values_[ numTopologies ];
750 template<
class ctype,
int dim >
759 return container() ( type );
765 return container().simplex();
771 return container().cube();
774 static Iterator begin () {
return container().begin(); }
775 static Iterator end () {
return container().end(); }
778 DUNE_EXPORT static const ReferenceElementContainer< ctype, dim > &container ()
780 static ReferenceElementContainer< ctype, dim > container;
An implementation of the Geometry interface for affine geometries.
Fallback implementation of the std::array class (a static array)
Implementation of the Geometry interface for affine geometries.
Definition: affinegeometry.hh:461
Unique label for each type of entities that can occur in DUNE grids.
Definition: type.hh:268
unsigned int dim() const
Return dimension of the type.
Definition: type.hh:565
This class provides access to geometric and topological properties of a reference element.
Definition: referenceelements.hh:354
int size(int c) const
number of subentities of codimension c
Definition: referenceelements.hh:381
int size(int i, int c, int cc) const
number of subentities of codimension cc of subentity (i,c)
Definition: referenceelements.hh:398
const FieldVector< ctype, dim > & position(int i, int c) const
position of the barycenter of entity (i,c)
Definition: referenceelements.hh:449
ctype volume() const
obtain the volume of the reference element
Definition: referenceelements.hh:486
int subEntity(int i, int c, int ii, int cc) const
obtain number of ii-th subentity with codim cc of (i,c)
Definition: referenceelements.hh:417
bool checkInside(const FieldVector< ctype, dim > &local) const
check if a coordinate is in the reference element
Definition: referenceelements.hh:462
const GeometryType & type(int i, int c) const
obtain the type of subentity (i,c)
Definition: referenceelements.hh:431
Codim< codim >::Geometry geometry(int i) const
obtain the embedding of subentity (i,codim) into the reference element
Definition: referenceelements.hh:480
const FieldVector< ctype, dim > & integrationOuterNormal(int face) const
obtain the integration outer normal of the reference element
Definition: referenceelements.hh:498
const GeometryType & type() const
obtain the type of this reference element
Definition: referenceelements.hh:438
Implements a matrix constructed from a given type representing a field and compile-time given number ...
A static for loop for template meta-programming.
Implements a vector constructed from a given type representing a field and a compile-time given size.
Dune namespace.
Definition: alignment.hh:11
A static loop using TMP.
Definition: forloop.hh:67
Collection of types depending on the codimension.
Definition: referenceelements.hh:372
AffineGeometry< ctype, dim-codim, dim > Geometry
type of geometry embedding a subentity into the reference element
Definition: referenceelements.hh:374
topological information about the subentities of a reference element
Definition: referenceelements.hh:573
Class providing access to the singletons of the reference elements.
Definition: referenceelements.hh:752
static const ReferenceElement< ctype, dim > & simplex()
get simplex reference elements
Definition: referenceelements.hh:763
static const ReferenceElement< ctype, dim > & cube()
get hypercube reference elements
Definition: referenceelements.hh:769
static const ReferenceElement< ctype, dim > & general(const GeometryType &type)
get general reference elements
Definition: referenceelements.hh:757
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.
#define DUNE_UNUSED_PARAMETER(parm)
A macro to mark intentionally unused function parameters with.
Definition: unused.hh:18
Definition of macros controlling symbol visibility at the ABI level.
#define DUNE_EXPORT
Export a symbol as part of the public ABI.
Definition: visibility.hh:18