3 #ifndef DUNE_GEOMETRY_REFERENCEELEMENTS_HH 4 #define DUNE_GEOMETRY_REFERENCEELEMENTS_HH 14 #include <dune/common/array.hh> 15 #include <dune/common/fmatrix.hh> 16 #include <dune/common/fvector.hh> 17 #include <dune/common/hybridutilities.hh> 18 #include <dune/common/std/utility.hh> 19 #include <dune/common/typetraits.hh> 20 #include <dune/common/visibility.hh> 21 #include <dune/common/unused.hh> 32 template<
class ctype,
int dim >
35 template<
class ctype,
int dim >
44 unsigned int size (
unsigned int topologyId,
int dim,
int codim );
55 unsigned int subTopologyId (
unsigned int topologyId,
int dim,
int codim,
unsigned int i );
62 void subTopologyNumbering (
unsigned int topologyId,
int dim,
int codim,
unsigned int i,
int subcodim,
63 unsigned int *beginOut,
unsigned int *endOut );
71 template<
class ct,
int cdim >
73 checkInside (
unsigned int topologyId,
int dim,
const FieldVector< ct, cdim > &x, ct tolerance, ct factor = ct( 1 ) )
75 assert( (dim >= 0) && (dim <= cdim) );
80 const ct baseFactor = (
isPrism( topologyId, dim ) ? factor : factor - x[ dim-1 ]);
81 if( (x[ dim-1 ] > -tolerance) && (factor - x[ dim-1 ] > -tolerance) )
82 return checkInside< ct, cdim >(
baseTopologyId( topologyId, dim ), dim-1, x, tolerance, baseFactor );
95 template<
class ct,
int cdim >
97 referenceCorners (
unsigned int topologyId,
int dim, FieldVector< ct, cdim > *corners )
99 assert( (dim >= 0) && (dim <= cdim) );
104 const unsigned int nBaseCorners
107 if(
isPrism( topologyId, dim ) )
109 std::copy( corners, corners + nBaseCorners, corners + nBaseCorners );
110 for(
unsigned int i = 0; i < nBaseCorners; ++i )
111 corners[ i+nBaseCorners ][ dim-1 ] = ct( 1 );
112 return 2*nBaseCorners;
116 corners[ nBaseCorners ] = FieldVector< ct, cdim >( ct( 0 ) );
117 corners[ nBaseCorners ][ dim-1 ] = ct( 1 );
118 return nBaseCorners+1;
123 *corners = FieldVector< ct, cdim >( ct( 0 ) );
146 template<
class ct,
int cdim >
148 referenceOrigins (
unsigned int topologyId,
int dim,
int codim, FieldVector< ct, cdim > *origins )
150 assert( (dim >= 0) && (dim <= cdim) );
152 assert( (codim >= 0) && (codim <= dim) );
157 if(
isPrism( topologyId, dim ) )
159 const unsigned int n = (codim < dim ?
referenceOrigins( baseId, dim-1, codim, origins ) : 0);
160 const unsigned int m =
referenceOrigins( baseId, dim-1, codim-1, origins+n );
161 for(
unsigned int i = 0; i < m; ++i )
163 origins[ n+m+i ] = origins[ n+i ];
164 origins[ n+m+i ][ dim-1 ] = ct( 1 );
173 origins[ m ] = FieldVector< ct, cdim >( ct( 0 ) );
174 origins[ m ][ dim-1 ] = ct( 1 );
183 origins[ 0 ] = FieldVector< ct, cdim >( ct( 0 ) );
193 template<
class ct,
int cdim,
int mydim >
196 FieldVector< ct, cdim > *origins,
197 FieldMatrix< ct, mydim, cdim > *jacobianTransposeds )
199 assert( (0 <= codim) && (codim <= dim) && (dim <= cdim) );
200 assert( (dim - codim <= mydim) && (mydim <= cdim) );
206 if(
isPrism( topologyId, dim ) )
208 const unsigned int n = (codim < dim ?
referenceEmbeddings( baseId, dim-1, codim, origins, jacobianTransposeds ) : 0);
209 for(
unsigned int i = 0; i < n; ++i )
210 jacobianTransposeds[ i ][ dim-codim-1 ][ dim-1 ] = ct( 1 );
212 const unsigned int m =
referenceEmbeddings( baseId, dim-1, codim-1, origins+n, jacobianTransposeds+n );
213 std::copy( origins+n, origins+n+m, origins+n+m );
214 std::copy( jacobianTransposeds+n, jacobianTransposeds+n+m, jacobianTransposeds+n+m );
215 for(
unsigned int i = 0; i < m; ++i )
216 origins[ n+m+i ][ dim-1 ] = ct( 1 );
222 const unsigned int m =
referenceEmbeddings( baseId, dim-1, codim-1, origins, jacobianTransposeds );
225 origins[ m ] = FieldVector< ct, cdim >( ct( 0 ) );
226 origins[ m ][ dim-1 ] = ct( 1 );
227 jacobianTransposeds[ m ] = FieldMatrix< ct, mydim, cdim >( ct( 0 ) );
232 const unsigned int n =
referenceEmbeddings( baseId, dim-1, codim, origins+m, jacobianTransposeds+m );
233 for(
unsigned int i = 0; i < n; ++i )
235 for(
int k = 0; k < dim-1; ++k )
236 jacobianTransposeds[ m+i ][ dim-codim-1 ][ k ] = -origins[ m+i ][ k ];
237 jacobianTransposeds[ m+i ][ dim-codim-1 ][ dim-1 ] = ct( 1 );
245 origins[ 0 ] = FieldVector< ct, cdim >( ct( 0 ) );
246 jacobianTransposeds[ 0 ] = FieldMatrix< ct, mydim, cdim >( ct( 0 ) );
247 for(
int k = 0; k < dim; ++k )
248 jacobianTransposeds[ 0 ][ k ][ k ] = ct( 1 );
258 template<
class ct,
int cdim >
261 const FieldVector< ct, cdim > *origins,
262 FieldVector< ct, cdim > *normals )
264 assert( (dim > 0) && (dim <= cdim) );
270 if(
isPrism( topologyId, dim ) )
272 const unsigned int numBaseFaces
275 for(
unsigned int i = 0; i < 2; ++i )
277 normals[ numBaseFaces+i ] = FieldVector< ct, cdim >( ct( 0 ) );
278 normals[ numBaseFaces+i ][ dim-1 ] = ct( 2*
int( i )-1 );
281 return numBaseFaces+2;
285 normals[ 0 ] = FieldVector< ct, cdim >( ct( 0 ) );
286 normals[ 0 ][ dim-1 ] = ct( -1 );
288 const unsigned int numBaseFaces
290 for(
unsigned int i = 1; i <= numBaseFaces; ++i )
291 normals[ i ][ dim-1 ] = normals[ i ]*origins[ i ];
293 return numBaseFaces+1;
298 for(
unsigned int i = 0; i < 2; ++i )
300 normals[ i ] = FieldVector< ct, cdim >( ct( 0 ) );
301 normals[ i ][ 0 ] = ct( 2*
int( i )-1 );
308 template<
class ct,
int cdim >
311 FieldVector< ct, cdim > *normals )
313 assert( (dim > 0) && (dim <= cdim) );
315 FieldVector< ct, cdim > *origins
316 =
new FieldVector< ct, cdim >[
size( topologyId, dim, 1 ) ];
319 const unsigned int numFaces
321 assert( numFaces ==
size( topologyId, dim, 1 ) );
353 template<
class ctype,
int dim >
367 template<
int codim >
struct CreateGeometries;
371 template<
int codim >
384 assert( (c >= 0) && (c <= dim) );
385 return info_[ c ].size();
399 int size (
int i,
int c,
int cc )
const 401 assert( (i >= 0) && (i <
size( c )) );
402 return info_[ c ][ i ].size( cc );
420 assert( (i >= 0) && (i <
size( c )) );
421 return info_[ c ][ i ].number( ii, cc );
434 assert( (i >= 0) && (i <
size( c )) );
435 return info_[ c ][ i ].type();
450 const FieldVector< ctype, dim > &
position(
int i,
int c )
const 452 assert( (c >= 0) && (c <= dim) );
453 return baryCenters_[ c ][ i ];
465 const ctype tolerance = ctype( 64 ) * std::numeric_limits< ctype >::epsilon();
466 return Impl::template checkInside< ctype, dim >( type().id(), dim, local, tolerance );
480 template<
int codim >
483 return std::get< codim >( geometries_ )[ i ];
501 assert( (face >= 0) && (face <
int( integrationNormals_.size() )) );
502 return integrationNormals_[ face ];
506 void initialize (
unsigned int topologyId )
511 for(
int codim = 0; codim <= dim; ++codim )
514 info_[ codim ].resize( size );
515 for(
unsigned int i = 0; i <
size; ++i )
516 info_[ codim ][ i ].initialize( topologyId, codim, i );
520 const unsigned int numVertices =
size( dim );
521 baryCenters_[ dim ].resize( numVertices );
525 for(
int codim = 0; codim < dim; ++codim )
527 baryCenters_[ codim ].resize(
size(codim) );
528 for(
int i = 0; i <
size( codim ); ++i )
530 baryCenters_[ codim ][ i ] = FieldVector< ctype, dim >( ctype( 0 ) );
531 const unsigned int numCorners =
size( i, codim, dim );
532 for(
unsigned int j = 0; j < numCorners; ++j )
533 baryCenters_[ codim ][ i ] += baryCenters_[ dim ][ subEntity( i, codim, j, dim ) ];
534 baryCenters_[ codim ][ i ] *= ctype( 1 ) / ctype( numCorners );
539 volume_ = Impl::template referenceVolume< ctype >( topologyId, dim );
544 integrationNormals_.resize(
size( 1 ) );
549 Hybrid::forEach( Std::make_index_sequence< dim+1 >{}, [ & ](
auto i ){ CreateGeometries< i >::apply( *
this, geometries_ ); } );
552 template<
int... codim >
553 static std::tuple< std::vector< typename Codim< codim >::Geometry >... >
554 makeGeometryTable ( std::integer_sequence< int, codim... > );
557 typedef decltype( makeGeometryTable( std::make_integer_sequence< int, dim+1 >() ) ) GeometryTable;
562 std::vector< FieldVector< ctype, dim > > baryCenters_[ dim+1 ];
563 std::vector< FieldVector< ctype, dim > > integrationNormals_;
566 GeometryTable geometries_;
568 std::vector< SubEntityInfo > info_[ dim+1 ];
572 template<
class ctype,
int dim >
576 : numbering_( nullptr )
578 std::fill( offset_.begin(), offset_.end(), 0 );
582 : offset_( other.offset_ ),
585 numbering_ = allocate();
586 std::copy( other.numbering_, other.numbering_ + capacity(), numbering_ );
594 offset_ = other.offset_;
596 deallocate( numbering_ );
597 numbering_ = allocate();
598 std::copy( other.numbering_, other.numbering_ + capacity(), numbering_ );
605 assert( (cc >= codim()) && (cc <= dim) );
606 return (offset_[ cc+1 ] - offset_[ cc ]);
611 assert( (ii >= 0) && (ii <
size( cc )) );
612 return numbering_[ offset_[ cc ] + ii ];
617 void initialize (
unsigned int topologyId,
int codim,
unsigned int i )
623 for(
int cc = 0; cc <= codim; ++cc )
625 for(
int cc = codim; cc <= dim; ++cc )
626 offset_[ cc+1 ] = offset_[ cc ] +
Impl::size( subId, dim-codim, cc-codim );
629 deallocate( numbering_ );
630 numbering_ = allocate();
631 for(
int cc = codim; cc <= dim; ++cc )
636 int codim ()
const {
return dim - type().dim(); }
638 unsigned int *
allocate () {
return (capacity() != 0 ?
new unsigned int[ capacity() ] :
nullptr); }
640 unsigned int capacity ()
const {
return offset_[ dim+1 ]; }
643 unsigned int *numbering_;
644 std::array< unsigned int, dim+2 > offset_;
649 template<
class ctype,
int dim >
650 template<
int codim >
663 DUNE_UNUSED_PARAMETER(i);
669 const int size = refElement.
size( codim );
670 std::vector< FieldVector< ctype, dim > > origins( size );
671 std::vector< FieldMatrix< ctype, dim - codim, dim > > jacobianTransposeds( size );
674 std::get< codim >( geometries ).reserve( size );
675 for(
int i = 0; i <
size; ++i )
677 typename Codim< codim >::Geometry geometry( subRefElement( refElement, i, std::integral_constant< int, codim >() ), origins[ i ], jacobianTransposeds[ i ] );
678 std::get< codim >( geometries ).push_back( geometry );
688 template<
class ctype,
int dim >
699 for(
unsigned int topologyId = 0; topologyId <
numTopologies; ++topologyId )
700 values_[ topologyId ].initialize( topologyId );
705 assert( type.
dim() == dim );
706 return values_[ type.
id() ];
714 const value_type &
cube ()
const 729 const_iterator
begin ()
const {
return values_; }
751 template<
class ctype,
int dim >
760 return container() ( type );
766 return container().simplex();
772 return container().cube();
775 static Iterator
begin () {
return container().begin(); }
776 static Iterator
end () {
return container().end(); }
788 #endif // #ifndef DUNE_GEOMETRY_REFERENCEELEMENTS_HH ReferenceElementContainer< ctype, dim >::const_iterator Iterator
Definition: referenceelements.hh:754
Definition: referenceelements.hh:33
Codim< codim >::Geometry geometry(int i) const
obtain the embedding of subentity (i,codim) into the reference element
Definition: referenceelements.hh:481
unsigned int referenceIntegrationOuterNormals(unsigned int topologyId, int dim, const FieldVector< ct, cdim > *origins, FieldVector< ct, cdim > *normals)
Definition: referenceelements.hh:260
Class providing access to the singletons of the reference elements.
Definition: affinegeometry.hh:28
const FieldVector< ctype, dim > & integrationOuterNormal(int face) const
obtain the integration outer normal of the reference element
Definition: referenceelements.hh:499
Unique label for each type of entities that can occur in DUNE grids.
Definition: type.hh:267
unsigned int dim() const
Return dimension of the type.
Definition: type.hh:565
Collection of types depending on the codimension.
Definition: referenceelements.hh:372
const_iterator end() const
Definition: referenceelements.hh:730
static const ReferenceElement< ctype, dim > & simplex()
get simplex reference elements
Definition: referenceelements.hh:764
const GeometryType & type() const
Definition: referenceelements.hh:615
const value_type & pyramid() const
Definition: referenceelements.hh:719
bool checkInside(unsigned int topologyId, int dim, const FieldVector< ct, cdim > &x, ct tolerance, ct factor=ct(1))
Definition: referenceelements.hh:73
Implementation of the Geometry interface for affine geometries.
Definition: affinegeometry.hh:460
This class provides access to geometric and topological properties of a reference element...
Definition: affinegeometry.hh:25
ReferenceElement< ctype, dim > value_type
Definition: referenceelements.hh:694
int size(int c) const
number of subentities of codimension c
Definition: referenceelements.hh:382
void subTopologyNumbering(unsigned int topologyId, int dim, int codim, unsigned int i, int subcodim, unsigned int *beginOut, unsigned int *endOut)
Definition: referenceelements.cc:85
unsigned int size(unsigned int topologyId, int dim, int codim)
Compute the number of subentities of a given codimension.
Definition: referenceelements.cc:16
Definition: affinegeometry.hh:18
static unsigned int numTopologies(int dim) noexcept
obtain the number of topologies of a given dimension
Definition: type.hh:93
const FieldVector< ctype, dim > & position(int i, int c) const
position of the barycenter of entity (i,c)
Definition: referenceelements.hh:450
A unique label for each type of element that can occur in a grid.
const_iterator begin() const
Definition: referenceelements.hh:729
const value_type & cube() const
Definition: referenceelements.hh:714
int size(int i, int c, int cc) const
number of subentities of codimension cc of subentity (i,c)
Definition: referenceelements.hh:399
void initialize(unsigned int topologyId, int codim, unsigned int i)
Definition: referenceelements.hh:617
const value_type & prism() const
Definition: referenceelements.hh:724
const GeometryType & type() const
obtain the type of this reference element
Definition: referenceelements.hh:439
An implementation of the Geometry interface for affine geometries.
int size(int cc) const
Definition: referenceelements.hh:603
static Iterator begin()
Definition: referenceelements.hh:775
static unsigned int baseTopologyId(unsigned int topologyId, int dim, int codim=1) noexcept
obtain the base topology of a given codimension
Definition: type.hh:160
AffineGeometry< ctype, dim-codim, dim > Geometry
type of geometry embedding a subentity into the reference element
Definition: referenceelements.hh:375
static const ReferenceElement< ctype, dim > & general(const GeometryType &type)
get general reference elements
Definition: referenceelements.hh:758
static const ReferenceElement< ctype, dim > & cube()
get hypercube reference elements
Definition: referenceelements.hh:770
unsigned int referenceCorners(unsigned int topologyId, int dim, FieldVector< ct, cdim > *corners)
Definition: referenceelements.hh:97
unsigned int capacity() const
Definition: referenceelements.hh:640
static bool isPrism(unsigned int topologyId, int dim, int codim=0) noexcept
check whether a prism construction was used to create a given codimension
Definition: type.hh:127
const value_type * const_iterator
Definition: referenceelements.hh:695
unsigned int * allocate()
Definition: referenceelements.hh:638
const value_type & simplex() const
Definition: referenceelements.hh:709
int number(int ii, int cc) const
Definition: referenceelements.hh:609
~SubEntityInfo()
Definition: referenceelements.hh:589
SubEntityInfo()
Definition: referenceelements.hh:575
unsigned int referenceOrigins(unsigned int topologyId, int dim, int codim, FieldVector< ct, cdim > *origins)
Definition: referenceelements.hh:148
unsigned long referenceVolumeInverse(unsigned int topologyId, int dim)
Definition: referenceelements.cc:171
static Iterator end()
Definition: referenceelements.hh:776
unsigned int id() const
Return the topology id the type.
Definition: type.hh:570
bool checkInside(const FieldVector< ctype, dim > &local) const
check if a coordinate is in the reference element
Definition: referenceelements.hh:463
ct referenceVolume(unsigned int topologyId, int dim)
Definition: referenceelements.hh:136
int codim() const
Definition: referenceelements.hh:636
const GeometryType & type(int i, int c) const
obtain the type of subentity (i,c)
Definition: referenceelements.hh:432
ReferenceElementContainer()
Definition: referenceelements.hh:697
unsigned int subTopologyId(unsigned int topologyId, int dim, int codim, unsigned int i)
Compute the topology id of a given subentity.
Definition: referenceelements.cc:47
ctype volume() const
obtain the volume of the reference element
Definition: referenceelements.hh:487
void deallocate(unsigned int *ptr)
Definition: referenceelements.hh:639
unsigned int referenceEmbeddings(unsigned int topologyId, int dim, int codim, FieldVector< ct, cdim > *origins, FieldMatrix< ct, mydim, cdim > *jacobianTransposeds)
Definition: referenceelements.hh:195
SubEntityInfo(const SubEntityInfo &other)
Definition: referenceelements.hh:581
topological information about the subentities of a reference element
Definition: referenceelements.hh:573
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:418