3#ifndef DUNE_GEOMETRY_REFERENCEELEMENTIMPLEMENTATION_HH
4#define DUNE_GEOMETRY_REFERENCEELEMENTIMPLEMENTATION_HH
18#include <dune/common/hybridutilities.hh>
20#include <dune/common/iteratorrange.hh>
24#include <dune/geometry/referenceelement.hh>
41 template<
class ctype,
int dim >
42 class ReferenceElementContainer;
45 template<
class ctype,
int dim >
46 struct ReferenceElements;
53 using Dune::Impl::isPrism;
54 using Dune::Impl::isPyramid;
55 using Dune::Impl::baseTopologyId;
56 using Dune::Impl::prismConstruction;
57 using Dune::Impl::pyramidConstruction;
58 using Dune::Impl::numTopologies;
61 unsigned int size (
unsigned int topologyId,
int dim,
int codim );
72 unsigned int subTopologyId (
unsigned int topologyId,
int dim,
int codim,
unsigned int i );
79 void subTopologyNumbering (
unsigned int topologyId,
int dim,
int codim,
unsigned int i,
int subcodim,
80 unsigned int *beginOut,
unsigned int *endOut );
88 template<
class ct,
int cdim >
90 checkInside (
unsigned int topologyId,
int dim,
const FieldVector< ct, cdim > &x, ct tolerance, ct factor = ct( 1 ) )
92 assert( (dim >= 0) && (dim <= cdim) );
93 assert( topologyId < numTopologies( dim ) );
97 const ct baseFactor = (isPrism( topologyId, dim ) ? factor : factor - x[ dim-1 ]);
98 if( (x[ dim-1 ] > -tolerance) && (factor - x[ dim-1 ] > -tolerance) )
99 return checkInside< ct, cdim >( baseTopologyId( topologyId, dim ), dim-1, x, tolerance, baseFactor );
112 template<
class ct,
int cdim >
114 referenceCorners (
unsigned int topologyId,
int dim, FieldVector< ct, cdim > *corners )
116 assert( (dim >= 0) && (dim <= cdim) );
117 assert( topologyId < numTopologies( dim ) );
121 const unsigned int nBaseCorners
122 = referenceCorners( baseTopologyId( topologyId, dim ), dim-1, corners );
123 assert( nBaseCorners == size( baseTopologyId( topologyId, dim ), dim-1, dim-1 ) );
124 if( isPrism( topologyId, dim ) )
126 std::copy( corners, corners + nBaseCorners, corners + nBaseCorners );
127 for(
unsigned int i = 0; i < nBaseCorners; ++i )
128 corners[ i+nBaseCorners ][ dim-1 ] = ct( 1 );
129 return 2*nBaseCorners;
133 corners[ nBaseCorners ] = FieldVector< ct, cdim >( ct( 0 ) );
134 corners[ nBaseCorners ][ dim-1 ] = ct( 1 );
135 return nBaseCorners+1;
140 *corners = FieldVector< ct, cdim >( ct( 0 ) );
150 unsigned long referenceVolumeInverse (
unsigned int topologyId,
int dim );
153 inline ct referenceVolume (
unsigned int topologyId,
int dim )
155 return ct( 1 ) / ct( referenceVolumeInverse( topologyId, dim ) );
163 template<
class ct,
int cdim >
165 referenceOrigins (
unsigned int topologyId,
int dim,
int codim, FieldVector< ct, cdim > *origins )
167 assert( (dim >= 0) && (dim <= cdim) );
168 assert( topologyId < numTopologies( dim ) );
169 assert( (codim >= 0) && (codim <= dim) );
173 const unsigned int baseId = baseTopologyId( topologyId, dim );
174 if( isPrism( topologyId, dim ) )
176 const unsigned int n = (codim < dim ? referenceOrigins( baseId, dim-1, codim, origins ) : 0);
177 const unsigned int m = referenceOrigins( baseId, dim-1, codim-1, origins+n );
178 for(
unsigned int i = 0; i < m; ++i )
180 origins[ n+m+i ] = origins[ n+i ];
181 origins[ n+m+i ][ dim-1 ] = ct( 1 );
187 const unsigned int m = referenceOrigins( baseId, dim-1, codim-1, origins );
190 origins[ m ] = FieldVector< ct, cdim >( ct( 0 ) );
191 origins[ m ][ dim-1 ] = ct( 1 );
195 return m+referenceOrigins( baseId, dim-1, codim, origins+m );
200 origins[ 0 ] = FieldVector< ct, cdim >( ct( 0 ) );
210 template<
class ct,
int cdim,
int mydim >
212 referenceEmbeddings (
unsigned int topologyId,
int dim,
int codim,
213 FieldVector< ct, cdim > *origins,
214 FieldMatrix< ct, mydim, cdim > *jacobianTransposeds )
216 assert( (0 <= codim) && (codim <= dim) && (dim <= cdim) );
217 assert( (dim - codim <= mydim) && (mydim <= cdim) );
218 assert( topologyId < numTopologies( dim ) );
222 const unsigned int baseId = baseTopologyId( topologyId, dim );
223 if( isPrism( topologyId, dim ) )
225 const unsigned int n = (codim < dim ? referenceEmbeddings( baseId, dim-1, codim, origins, jacobianTransposeds ) : 0);
226 for(
unsigned int i = 0; i < n; ++i )
227 jacobianTransposeds[ i ][ dim-codim-1 ][ dim-1 ] = ct( 1 );
229 const unsigned int m = referenceEmbeddings( baseId, dim-1, codim-1, origins+n, jacobianTransposeds+n );
230 std::copy( origins+n, origins+n+m, origins+n+m );
231 std::copy( jacobianTransposeds+n, jacobianTransposeds+n+m, jacobianTransposeds+n+m );
232 for(
unsigned int i = 0; i < m; ++i )
233 origins[ n+m+i ][ dim-1 ] = ct( 1 );
239 const unsigned int m = referenceEmbeddings( baseId, dim-1, codim-1, origins, jacobianTransposeds );
242 origins[ m ] = FieldVector< ct, cdim >( ct( 0 ) );
243 origins[ m ][ dim-1 ] = ct( 1 );
244 jacobianTransposeds[ m ] = FieldMatrix< ct, mydim, cdim >( ct( 0 ) );
249 const unsigned int n = referenceEmbeddings( baseId, dim-1, codim, origins+m, jacobianTransposeds+m );
250 for(
unsigned int i = 0; i < n; ++i )
252 for(
int k = 0; k < dim-1; ++k )
253 jacobianTransposeds[ m+i ][ dim-codim-1 ][ k ] = -origins[ m+i ][ k ];
254 jacobianTransposeds[ m+i ][ dim-codim-1 ][ dim-1 ] = ct( 1 );
262 origins[ 0 ] = FieldVector< ct, cdim >( ct( 0 ) );
263 jacobianTransposeds[ 0 ] = FieldMatrix< ct, mydim, cdim >( ct( 0 ) );
264 for(
int k = 0; k < dim; ++k )
265 jacobianTransposeds[ 0 ][ k ][ k ] = ct( 1 );
275 template<
class ct,
int cdim >
277 referenceIntegrationOuterNormals (
unsigned int topologyId,
int dim,
278 const FieldVector< ct, cdim > *origins,
279 FieldVector< ct, cdim > *normals )
281 assert( (dim > 0) && (dim <= cdim) );
282 assert( topologyId < numTopologies( dim ) );
286 const unsigned int baseId = baseTopologyId( topologyId, dim );
287 if( isPrism( topologyId, dim ) )
289 const unsigned int numBaseFaces
290 = referenceIntegrationOuterNormals( baseId, dim-1, origins, normals );
292 for(
unsigned int i = 0; i < 2; ++i )
294 normals[ numBaseFaces+i ] = FieldVector< ct, cdim >( ct( 0 ) );
295 normals[ numBaseFaces+i ][ dim-1 ] = ct( 2*
int( i )-1 );
298 return numBaseFaces+2;
302 normals[ 0 ] = FieldVector< ct, cdim >( ct( 0 ) );
303 normals[ 0 ][ dim-1 ] = ct( -1 );
305 const unsigned int numBaseFaces
306 = referenceIntegrationOuterNormals( baseId, dim-1, origins+1, normals+1 );
307 for(
unsigned int i = 1; i <= numBaseFaces; ++i )
308 normals[ i ][ dim-1 ] = normals[ i ]*origins[ i ];
310 return numBaseFaces+1;
315 for(
unsigned int i = 0; i < 2; ++i )
317 normals[ i ] = FieldVector< ct, cdim >( ct( 0 ) );
318 normals[ i ][ 0 ] = ct( 2*
int( i )-1 );
325 template<
class ct,
int cdim >
327 referenceIntegrationOuterNormals (
unsigned int topologyId,
int dim,
328 FieldVector< ct, cdim > *normals )
330 assert( (dim > 0) && (dim <= cdim) );
332 FieldVector< ct, cdim > *origins
333 =
new FieldVector< ct, cdim >[ size( topologyId, dim, 1 ) ];
334 referenceOrigins( topologyId, dim, 1, origins );
336 const unsigned int numFaces
337 = referenceIntegrationOuterNormals( topologyId, dim, origins, normals );
338 assert( numFaces == size( topologyId, dim, 1 ) );
370 template<
class ctype_,
int dim >
371 class ReferenceElementImplementation
377 using ctype = ctype_;
380 using CoordinateField = ctype;
386 static constexpr int dimension = dim;
389 typedef ctype Volume;
393 friend class Impl::ReferenceElementContainer< ctype, dim >;
395 struct SubEntityInfo;
397 template<
int codim >
struct CreateGeometries;
401 template<
int codim >
405 typedef AffineGeometry< ctype, dim-codim, dim > Geometry;
409 ReferenceElementImplementation (
const ReferenceElementImplementation& ) =
delete;
412 ReferenceElementImplementation& operator= (
const ReferenceElementImplementation& ) =
delete;
415 ReferenceElementImplementation () =
default;
421 int size (
int c )
const
423 assert( (c >= 0) && (c <= dim) );
424 return info_[ c ].size();
438 int size (
int i,
int c,
int cc )
const
440 assert( (i >= 0) && (i < size( c )) );
441 return info_[ c ][ i ].size( cc );
457 int subEntity (
int i,
int c,
int ii,
int cc )
const
459 assert( (i >= 0) && (i < size( c )) );
460 return info_[ c ][ i ].number( ii, cc );
478 auto subEntities (
int i,
int c,
int cc )
const
480 assert( (i >= 0) && (i < size( c )) );
481 return info_[ c ][ i ].numbers( cc );
494 assert( (i >= 0) && (i < size( c )) );
495 return info_[ c ][ i ].type();
499 const GeometryType &type ()
const {
return type( 0, 0 ); }
510 const Coordinate &position(
int i,
int c )
const
512 assert( (c >= 0) && (c <= dim) );
513 return baryCenters_[ c ][ i ];
523 bool checkInside (
const Coordinate &local )
const
525 const ctype tolerance = ctype( 64 ) * std::numeric_limits< ctype >::epsilon();
526 return Impl::template checkInside< ctype, dim >( type().
id(), dim, local, tolerance );
540 template<
int codim >
541 typename Codim< codim >::Geometry geometry (
int i )
const
543 return std::get< codim >( geometries_ )[ i ];
547 Volume volume ()
const
559 const Coordinate &integrationOuterNormal (
int face )
const
561 assert( (face >= 0) && (face <
int( integrationNormals_.size() )) );
562 return integrationNormals_[ face ];
566 void initialize (
unsigned int topologyId )
568 assert( topologyId < Impl::numTopologies( dim ) );
571 for(
int codim = 0; codim <= dim; ++codim )
573 const unsigned int size = Impl::size( topologyId, dim, codim );
574 info_[ codim ].resize( size );
575 for(
unsigned int i = 0; i < size; ++i )
576 info_[ codim ][ i ].initialize( topologyId, codim, i );
580 const unsigned int numVertices = size( dim );
581 baryCenters_[ dim ].resize( numVertices );
582 Impl::referenceCorners( topologyId, dim, &(baryCenters_[ dim ][ 0 ]) );
585 for(
int codim = 0; codim < dim; ++codim )
587 baryCenters_[ codim ].resize( size(codim) );
588 for(
int i = 0; i < size( codim ); ++i )
590 baryCenters_[ codim ][ i ] = Coordinate( ctype( 0 ) );
591 const unsigned int numCorners = size( i, codim, dim );
592 for(
unsigned int j = 0; j < numCorners; ++j )
593 baryCenters_[ codim ][ i ] += baryCenters_[ dim ][ subEntity( i, codim, j, dim ) ];
594 baryCenters_[ codim ][ i ] *= ctype( 1 ) / ctype( numCorners );
599 volume_ = Impl::template referenceVolume< ctype >( topologyId, dim );
604 integrationNormals_.resize( size( 1 ) );
605 Impl::referenceIntegrationOuterNormals( topologyId, dim, &(integrationNormals_[ 0 ]) );
609 Hybrid::forEach( std::make_index_sequence< dim+1 >{}, [ & ](
auto i ){ CreateGeometries< i >::apply( *
this, geometries_ ); } );
612 template<
int... codim >
613 static std::tuple< std::vector< typename Codim< codim >::Geometry >... >
614 makeGeometryTable ( std::integer_sequence< int, codim... > );
617 typedef decltype( makeGeometryTable( std::make_integer_sequence< int, dim+1 >() ) ) GeometryTable;
622 std::vector< Coordinate > baryCenters_[ dim+1 ];
623 std::vector< Coordinate > integrationNormals_;
626 GeometryTable geometries_;
628 std::vector< SubEntityInfo > info_[ dim+1 ];
632 template<
class ctype,
int dim >
633 struct ReferenceElementImplementation< ctype, dim >::SubEntityInfo
638 static constexpr std::size_t maxSubEntityCount()
640 std::size_t maxCount=0;
641 for(std::size_t codim=0; codim<=dim; ++codim)
642 maxCount =
std::max(maxCount,
binomial(std::size_t(dim),codim)*(1 << codim));
646 using SubEntityFlags = std::bitset<maxSubEntityCount()>;
655 using iterator = Base::iterator;
656 using const_iterator = Base::const_iterator;
658 SubEntityRange(
const iterator& begin,
const iterator& end,
const SubEntityFlags& contains) :
660 containsPtr_(&contains),
666 containsPtr_(nullptr),
670 std::size_t size()
const
675 bool contains(std::size_t i)
const
677 return (*containsPtr_)[i];
681 const SubEntityFlags* containsPtr_;
689 : numbering_( nullptr )
691 std::fill( offset_.begin(), offset_.end(), 0 );
694 SubEntityInfo (
const SubEntityInfo &other )
695 : offset_( other.offset_ ),
696 type_( other.type_ ),
697 containsSubentity_( other.containsSubentity_ )
699 numbering_ = allocate();
700 std::copy( other.numbering_, other.numbering_ + capacity(), numbering_ );
703 ~SubEntityInfo () { deallocate( numbering_ ); }
705 const SubEntityInfo &operator= (
const SubEntityInfo &other )
708 offset_ = other.offset_;
710 deallocate( numbering_ );
711 numbering_ = allocate();
712 std::copy( other.numbering_, other.numbering_ + capacity(), numbering_ );
714 containsSubentity_ = other.containsSubentity_;
719 int size (
int cc )
const
721 assert( (cc >= 0) && (cc <= dim) );
722 return (offset_[ cc+1 ] - offset_[ cc ]);
725 int number (
int ii,
int cc )
const
727 assert( (ii >= 0) && (ii < size( cc )) );
728 return numbering_[ offset_[ cc ] + ii ];
731 auto numbers (
int cc )
const
733 return SubEntityRange( numbering_ + offset_[ cc ], numbering_ + offset_[ cc+1 ], containsSubentity_[cc]);
738 void initialize (
unsigned int topologyId,
int codim,
unsigned int i )
740 const unsigned int subId = Impl::subTopologyId( topologyId, dim, codim, i );
744 for(
int cc = 0; cc <= codim; ++cc )
746 for(
int cc = codim; cc <= dim; ++cc )
747 offset_[ cc+1 ] = offset_[ cc ] + Impl::size( subId, dim-codim, cc-codim );
750 deallocate( numbering_ );
751 numbering_ = allocate();
752 for(
int cc = codim; cc <= dim; ++cc )
753 Impl::subTopologyNumbering( topologyId, dim, codim, i, cc-codim, numbering_+offset_[ cc ], numbering_+offset_[ cc+1 ] );
756 for(std::size_t cc=0; cc<= dim; ++cc)
758 containsSubentity_[cc].reset();
759 for(std::size_t idx=0; idx<std::size_t(size(cc)); ++idx)
760 containsSubentity_[cc][number(idx,cc)] =
true;
765 int codim ()
const {
return dim - type().dim(); }
767 unsigned int *allocate () {
return (capacity() != 0 ?
new unsigned int[ capacity() ] :
nullptr); }
768 void deallocate (
unsigned int *ptr ) {
delete[] ptr; }
769 unsigned int capacity ()
const {
return offset_[ dim+1 ]; }
772 unsigned int *numbering_;
773 std::array< unsigned int, dim+2 > offset_;
775 std::array< SubEntityFlags, dim+1> containsSubentity_;
779 template<
class ctype,
int dim >
780 template<
int codim >
781 struct ReferenceElementImplementation< ctype, dim >::CreateGeometries
785 subRefElement(
const ReferenceElementImplementation< ctype, dim > &refElement,
int i, std::integral_constant< int, cc > )
791 subRefElement(
const ReferenceElementImplementation< ctype, dim > &refElement,
792 [[maybe_unused]]
int i, std::integral_constant<int, 0>)
797 static void apply (
const ReferenceElementImplementation< ctype, dim > &refElement, GeometryTable &geometries )
799 const int size = refElement.size( codim );
800 std::vector< FieldVector< ctype, dim > > origins( size );
801 std::vector< FieldMatrix< ctype, dim - codim, dim > > jacobianTransposeds( size );
802 Impl::referenceEmbeddings( refElement.type().id(), dim, codim, &(origins[ 0 ]), &(jacobianTransposeds[ 0 ]) );
804 std::get< codim >( geometries ).reserve( size );
805 for(
int i = 0; i < size; ++i )
807 typename Codim< codim >::Geometry geometry( subRefElement( refElement, i, std::integral_constant< int, codim >() ), origins[ i ], jacobianTransposeds[ i ] );
808 std::get< codim >( geometries ).push_back( geometry );
An implementation of the Geometry interface for affine geometries.
Simple range between a begin and an end iterator.
Definition: iteratorrange.hh:20
Various implementations of the power function for run-time and static arguments.
Traits for type conversions and type information.
GeometryType
Type representing VTK's entity geometry types.
Definition: common.hh:130
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.
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:266
auto max(ADLTag< 0 >, const V &v1, const V &v2)
implements binary Simd::max()
Definition: defaults.hh:79
Some useful basic math stuff.
Dune namespace.
Definition: alignedallocator.hh:11
static constexpr T binomial(const T &n, const T &k) noexcept
calculate the binomial coefficient n over k as a constexpr
Definition: math.hh:128
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.