│ │ │ │
│ │ │ │
│ │ │ │
│ │ │ │ -
5#ifndef DUNE_DEPRECATED_TOPOLOGY_HH
│ │ │ │ -
6#define DUNE_DEPRECATED_TOPOLOGY_HH
│ │ │ │ +
5#ifndef DUNE_GEOMETRY_REFERENCEELEMENTIMPLEMENTATION_HH
│ │ │ │ +
6#define DUNE_GEOMETRY_REFERENCEELEMENTIMPLEMENTATION_HH
│ │ │ │
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
15 struct PointDeprecationHelper
│ │ │ │ -
│ │ │ │ -
17 static const unsigned int dimension = 0;
│ │ │ │ -
18 static const unsigned int numCorners = 1;
│ │ │ │ -
│ │ │ │ -
20 static const unsigned int id = 0;
│ │ │ │ -
│ │ │ │ -
22 static std::string name () {
return "p"; }
│ │ │ │ -
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
18#include <dune/common/fmatrix.hh>
│ │ │ │ +
19#include <dune/common/fvector.hh>
│ │ │ │ +
20#include <dune/common/hybridutilities.hh>
│ │ │ │ +
21#include <dune/common/typetraits.hh>
│ │ │ │ +
22#include <dune/common/iteratorrange.hh>
│ │ │ │ +
23#include <dune/common/math.hh>
│ │ │ │
│ │ │ │ -
25 using Point [[deprecated(
"Use GeometryTypes::vertex instead.")]] = PointDeprecationHelper;
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
28 template<
class BaseTopology >
│ │ │ │ -
29 struct [[deprecated(
"Use GeometryTypes::prismaticExtension(GeometryType gt) instead.")]] Prism
│ │ │ │ -
│ │ │ │ -
31 static const unsigned int dimension = BaseTopology::dimension + 1;
│ │ │ │ -
32 static const unsigned int numCorners = 2 * BaseTopology::numCorners;
│ │ │ │ -
│ │ │ │ -
34 static const unsigned int id = BaseTopology::id | ((
unsigned int)prismConstruction << (dimension-1));
│ │ │ │ -
│ │ │ │ -
36 static std::string name () {
return BaseTopology::name() +
"l"; }
│ │ │ │ -
│ │ │ │ -
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │
│ │ │ │ -
40 template<
class BaseTopology >
│ │ │ │ -
41 struct [[deprecated(
"Use GeometryTypes::conicalExtension(GeometryType gt) instead.")]] Pyramid
│ │ │ │ -
│ │ │ │ -
43 static const unsigned int dimension = BaseTopology::dimension + 1;
│ │ │ │ -
44 static const unsigned int numCorners = BaseTopology::numCorners + 1;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
42 template<
class ctype,
int dim >
│ │ │ │ +
43 class ReferenceElementContainer;
│ │ │ │ +
│ │ │ │
│ │ │ │ -
46 static const unsigned int id = BaseTopology::id | ((
unsigned int)pyramidConstruction << (dimension-1));
│ │ │ │ -
│ │ │ │ -
48 static std::string name () {
return BaseTopology::name() +
"o"; }
│ │ │ │ -
│ │ │ │ +
46 template<
class ctype,
int dim >
│ │ │ │ +
47 struct ReferenceElements;
│ │ │ │ +
│ │ │ │ +
│ │ │ │
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
56 template<
class Topology >
│ │ │ │ -
57 struct [[deprecated(
"Use GeometryType::isSimplex() instead.")]] IsSimplex
│ │ │ │ -
58 :
public std::integral_constant< bool, (Topology::id >> 1) == 0 >
│ │ │ │ -
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
54 using Dune::Impl::isPrism;
│ │ │ │ +
55 using Dune::Impl::isPyramid;
│ │ │ │ +
56 using Dune::Impl::baseTopologyId;
│ │ │ │ +
57 using Dune::Impl::prismConstruction;
│ │ │ │ +
58 using Dune::Impl::pyramidConstruction;
│ │ │ │ +
59 using Dune::Impl::numTopologies;
│ │ │ │
│ │ │ │ -
61 template<
class Topology >
│ │ │ │ -
62 struct [[deprecated(
"Use GeometryType::isCube() instead.")]] IsCube
│ │ │ │ -
63 :
public std::integral_constant< bool, (Topology::id | 1) == (1 << Topology::dimension) - 1 >
│ │ │ │ -
│ │ │ │ +
62 unsigned int size (
unsigned int topologyId,
int dim,
int codim );
│ │ │ │ +
│ │ │ │ +
│ │ │ │
│ │ │ │ -
78 [[deprecated("Use GeometryType::isPrismatic() or GeometryType::isConical() instead.")]]
│ │ │ │ -
79 inline static bool isTopology ( TopologyConstruction construction, unsigned int topologyId, int dim, int codim = 0 ) noexcept
│ │ │ │ -
│ │ │ │ -
81 assert( (dim > 0) && (topologyId < numTopologies( dim )) );
│ │ │ │ -
82 assert( (0 <= codim) && (codim <= dim) );
│ │ │ │ -
83 return (codim >= (dim-1)) || (((topologyId >> (dim-codim-1)) & 1) == (unsigned int)construction);
│ │ │ │ -
│ │ │ │ +
73 unsigned int subTopologyId (
unsigned int topologyId,
int dim,
int codim,
unsigned int i );
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
80 void subTopologyNumbering (
unsigned int topologyId,
int dim,
int codim,
unsigned int i,
int subcodim,
│ │ │ │ +
81 unsigned int *beginOut,
unsigned int *endOut );
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
90 template< unsigned int dim >
│ │ │ │ -
91 struct [[deprecated("Use GeometryTypes::simplex(dim) instead.")]] SimplexTopology
│ │ │ │ -
│ │ │ │ -
93 typedef Pyramid< typename SimplexTopology< dim-1 >::type > type;
│ │ │ │ -
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
89 template<
class ct,
int cdim >
│ │ │ │ +
│ │ │ │ +
91 checkInside (
unsigned int topologyId,
int dim,
const FieldVector< ct, cdim > &x, ct tolerance, ct factor = ct( 1 ) )
│ │ │ │ +
│ │ │ │ +
93 assert( (dim >= 0) && (dim <= cdim) );
│ │ │ │ +
94 assert( topologyId < numTopologies( dim ) );
│ │ │ │
│ │ │ │ -
│ │ │ │ -
97 struct [[deprecated("Use GeometryTypes::simplex(dim) instead.")]] SimplexTopology< 0 >
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
107 template< unsigned int dim >
│ │ │ │ -
108 struct [[deprecated("Use GeometryTypes::cube(dim) instead.")]] CubeTopology
│ │ │ │ -
│ │ │ │ -
110 typedef Prism< typename CubeTopology< dim-1 >::type > type;
│ │ │ │ -
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
98 const ct baseFactor = (isPrism( topologyId, dim ) ? factor : factor - x[ dim-1 ]);
│ │ │ │ +
99 if( (x[ dim-1 ] > -tolerance) && (factor - x[ dim-1 ] > -tolerance) )
│ │ │ │ +
100 return checkInside< ct, cdim >( baseTopologyId( topologyId, dim ), dim-1, x, tolerance, baseFactor );
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │
│ │ │ │ -
│ │ │ │ -
114 struct [[deprecated("Use GeometryTypes::simplex(dim) instead.")]] CubeTopology< 0 >
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ +
113 template<
class ct,
int cdim >
│ │ │ │ +
│ │ │ │ +
115 referenceCorners (
unsigned int topologyId,
int dim, FieldVector< ct, cdim > *corners )
│ │ │ │ +
│ │ │ │ +
117 assert( (dim >= 0) && (dim <= cdim) );
│ │ │ │ +
118 assert( topologyId < numTopologies( dim ) );
│ │ │ │
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
124 template< unsigned int dim >
│ │ │ │ -
125 struct [[deprecated]] PyramidTopology
│ │ │ │ -
│ │ │ │ -
127 typedef Pyramid< typename CubeTopology< dim-1 >::type > type;
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
135 template< unsigned int dim >
│ │ │ │ -
136 struct [[deprecated]] PrismTopology
│ │ │ │ -
│ │ │ │ -
138 typedef Prism< typename SimplexTopology< dim-1 >::type > type;
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
122 const unsigned int nBaseCorners
│ │ │ │ +
123 = referenceCorners( baseTopologyId( topologyId, dim ), dim-1, corners );
│ │ │ │ +
124 assert( nBaseCorners ==
size( baseTopologyId( topologyId, dim ), dim-1, dim-1 ) );
│ │ │ │ +
125 if( isPrism( topologyId, dim ) )
│ │ │ │ +
│ │ │ │ +
127 std::copy( corners, corners + nBaseCorners, corners + nBaseCorners );
│ │ │ │ +
128 for(
unsigned int i = 0; i < nBaseCorners; ++i )
│ │ │ │ +
129 corners[ i+nBaseCorners ][ dim-1 ] = ct( 1 );
│ │ │ │ +
130 return 2*nBaseCorners;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
134 corners[ nBaseCorners ] = FieldVector< ct, cdim >( ct( 0 ) );
│ │ │ │ +
135 corners[ nBaseCorners ][ dim-1 ] = ct( 1 );
│ │ │ │ +
136 return nBaseCorners+1;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
141 *corners = FieldVector< ct, cdim >( ct( 0 ) );
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │
│ │ │ │ -
147 template< template< class > class Operation, int dim, class Topology = PointDeprecationHelper >
│ │ │ │ -
148 struct [[deprecated("Use IfGeometryType instead.")]] IfTopology
│ │ │ │ -
│ │ │ │ -
150 template< class... Args >
│ │ │ │ -
151 static auto apply ( unsigned int topologyId, Args &&... args )
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
154 return IfTopology< Operation, dim-1, Prism< Topology > >::apply( topologyId >> 1, std::forward< Args >( args )... );
│ │ │ │ -
│ │ │ │ -
156 return IfTopology< Operation, dim-1, Pyramid< Topology > >::apply( topologyId >> 1, std::forward< Args >( args )... );
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
151 unsigned long referenceVolumeInverse (
unsigned int topologyId,
int dim );
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
154 inline ct referenceVolume (
unsigned int topologyId,
int dim )
│ │ │ │ +
│ │ │ │ +
156 return ct( 1 ) / ct( referenceVolumeInverse( topologyId, dim ) );
│ │ │ │
│ │ │ │ -
│ │ │ │ +
│ │ │ │
│ │ │ │ -
160 template< template< class > class Operation, class Topology >
│ │ │ │ -
161 struct [[deprecated("Use IfGeometryType instead.")]] IfTopology< Operation, 0, Topology >
│ │ │ │ -
│ │ │ │ -
163 template<
class... Args >
│ │ │ │ -
164 static auto apply ([[maybe_unused]]
unsigned int topologyId, Args &&... args)
│ │ │ │ -
│ │ │ │ -
166 return Operation< Topology >::apply( std::forward< Args >( args )... );
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
164 template<
class ct,
int cdim >
│ │ │ │ +
│ │ │ │ +
166 referenceOrigins (
unsigned int topologyId,
int dim,
int codim, FieldVector< ct, cdim > *origins )
│ │ │ │ +
│ │ │ │ +
168 assert( (dim >= 0) && (dim <= cdim) );
│ │ │ │ +
169 assert( topologyId < numTopologies( dim ) );
│ │ │ │ +
170 assert( (codim >= 0) && (codim <= dim) );
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
174 const unsigned int baseId = baseTopologyId( topologyId, dim );
│ │ │ │ +
175 if( isPrism( topologyId, dim ) )
│ │ │ │ +
│ │ │ │ +
177 const unsigned int n = (codim < dim ? referenceOrigins( baseId, dim-1, codim, origins ) : 0);
│ │ │ │ +
178 const unsigned int m = referenceOrigins( baseId, dim-1, codim-1, origins+n );
│ │ │ │ +
179 for(
unsigned int i = 0; i < m; ++i )
│ │ │ │ +
│ │ │ │ +
181 origins[ n+m+i ] = origins[ n+i ];
│ │ │ │ +
182 origins[ n+m+i ][ dim-1 ] = ct( 1 );
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
188 const unsigned int m = referenceOrigins( baseId, dim-1, codim-1, origins );
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
191 origins[ m ] = FieldVector< ct, cdim >( ct( 0 ) );
│ │ │ │ +
192 origins[ m ][ dim-1 ] = ct( 1 );
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
196 return m+referenceOrigins( baseId, dim-1, codim, origins+m );
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
201 origins[ 0 ] = FieldVector< ct, cdim >( ct( 0 ) );
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
211 template<
class ct,
int cdim,
int mydim >
│ │ │ │ +
│ │ │ │ +
213 referenceEmbeddings (
unsigned int topologyId,
int dim,
int codim,
│ │ │ │ +
214 FieldVector< ct, cdim > *origins,
│ │ │ │ +
215 FieldMatrix< ct, mydim, cdim > *jacobianTransposeds )
│ │ │ │ +
│ │ │ │ +
217 assert( (0 <= codim) && (codim <= dim) && (dim <= cdim) );
│ │ │ │ +
218 assert( (dim - codim <= mydim) && (mydim <= cdim) );
│ │ │ │ +
219 assert( topologyId < numTopologies( dim ) );
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
223 const unsigned int baseId = baseTopologyId( topologyId, dim );
│ │ │ │ +
224 if( isPrism( topologyId, dim ) )
│ │ │ │ +
│ │ │ │ +
226 const unsigned int n = (codim < dim ? referenceEmbeddings( baseId, dim-1, codim, origins, jacobianTransposeds ) : 0);
│ │ │ │ +
227 for(
unsigned int i = 0; i < n; ++i )
│ │ │ │ +
228 jacobianTransposeds[ i ][ dim-codim-1 ][ dim-1 ] = ct( 1 );
│ │ │ │ +
│ │ │ │ +
230 const unsigned int m = referenceEmbeddings( baseId, dim-1, codim-1, origins+n, jacobianTransposeds+n );
│ │ │ │ +
231 std::copy( origins+n, origins+n+m, origins+n+m );
│ │ │ │ +
232 std::copy( jacobianTransposeds+n, jacobianTransposeds+n+m, jacobianTransposeds+n+m );
│ │ │ │ +
233 for(
unsigned int i = 0; i < m; ++i )
│ │ │ │ +
234 origins[ n+m+i ][ dim-1 ] = ct( 1 );
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
240 const unsigned int m = referenceEmbeddings( baseId, dim-1, codim-1, origins, jacobianTransposeds );
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
243 origins[ m ] = FieldVector< ct, cdim >( ct( 0 ) );
│ │ │ │ +
244 origins[ m ][ dim-1 ] = ct( 1 );
│ │ │ │ +
245 jacobianTransposeds[ m ] = FieldMatrix< ct, mydim, cdim >( ct( 0 ) );
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
250 const unsigned int n = referenceEmbeddings( baseId, dim-1, codim, origins+m, jacobianTransposeds+m );
│ │ │ │ +
251 for(
unsigned int i = 0; i < n; ++i )
│ │ │ │ +
│ │ │ │ +
253 for(
int k = 0; k < dim-1; ++k )
│ │ │ │ +
254 jacobianTransposeds[ m+i ][ dim-codim-1 ][ k ] = -origins[ m+i ][ k ];
│ │ │ │ +
255 jacobianTransposeds[ m+i ][ dim-codim-1 ][ dim-1 ] = ct( 1 );
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
263 origins[ 0 ] = FieldVector< ct, cdim >( ct( 0 ) );
│ │ │ │ +
264 jacobianTransposeds[ 0 ] = FieldMatrix< ct, mydim, cdim >( ct( 0 ) );
│ │ │ │ +
265 for(
int k = 0; k < dim; ++k )
│ │ │ │ +
266 jacobianTransposeds[ 0 ][ k ][ k ] = ct( 1 );
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
276 template<
class ct,
int cdim >
│ │ │ │ +
│ │ │ │ +
278 referenceIntegrationOuterNormals (
unsigned int topologyId,
int dim,
│ │ │ │ +
279 const FieldVector< ct, cdim > *origins,
│ │ │ │ +
280 FieldVector< ct, cdim > *normals )
│ │ │ │ +
│ │ │ │ +
282 assert( (dim > 0) && (dim <= cdim) );
│ │ │ │ +
283 assert( topologyId < numTopologies( dim ) );
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
287 const unsigned int baseId = baseTopologyId( topologyId, dim );
│ │ │ │ +
288 if( isPrism( topologyId, dim ) )
│ │ │ │ +
│ │ │ │ +
290 const unsigned int numBaseFaces
│ │ │ │ +
291 = referenceIntegrationOuterNormals( baseId, dim-1, origins, normals );
│ │ │ │ +
│ │ │ │ +
293 for(
unsigned int i = 0; i < 2; ++i )
│ │ │ │ +
│ │ │ │ +
295 normals[ numBaseFaces+i ] = FieldVector< ct, cdim >( ct( 0 ) );
│ │ │ │ +
296 normals[ numBaseFaces+i ][ dim-1 ] = ct( 2*
int( i )-1 );
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
299 return numBaseFaces+2;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
303 normals[ 0 ] = FieldVector< ct, cdim >( ct( 0 ) );
│ │ │ │ +
304 normals[ 0 ][ dim-1 ] = ct( -1 );
│ │ │ │ +
│ │ │ │ +
306 const unsigned int numBaseFaces
│ │ │ │ +
307 = referenceIntegrationOuterNormals( baseId, dim-1, origins+1, normals+1 );
│ │ │ │ +
308 for(
unsigned int i = 1; i <= numBaseFaces; ++i )
│ │ │ │ +
309 normals[ i ][ dim-1 ] = normals[ i ]*origins[ i ];
│ │ │ │ +
│ │ │ │ +
311 return numBaseFaces+1;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
316 for(
unsigned int i = 0; i < 2; ++i )
│ │ │ │ +
│ │ │ │ +
318 normals[ i ] = FieldVector< ct, cdim >( ct( 0 ) );
│ │ │ │ +
319 normals[ i ][ 0 ] = ct( 2*
int( i )-1 );
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
326 template<
class ct,
int cdim >
│ │ │ │ +
│ │ │ │ +
328 referenceIntegrationOuterNormals (
unsigned int topologyId,
int dim,
│ │ │ │ +
329 FieldVector< ct, cdim > *normals )
│ │ │ │ +
│ │ │ │ +
331 assert( (dim > 0) && (dim <= cdim) );
│ │ │ │ +
│ │ │ │ +
333 FieldVector< ct, cdim > *origins
│ │ │ │ +
334 =
new FieldVector< ct, cdim >[
size( topologyId, dim, 1 ) ];
│ │ │ │ +
335 referenceOrigins( topologyId, dim, 1, origins );
│ │ │ │ +
│ │ │ │ +
337 const unsigned int numFaces
│ │ │ │ +
338 = referenceIntegrationOuterNormals( topologyId, dim, origins, normals );
│ │ │ │ +
339 assert( numFaces ==
size( topologyId, dim, 1 ) );
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
371 template<
class ctype_,
int dim >
│ │ │ │ +
372 class ReferenceElementImplementation
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
378 using ctype = ctype_;
│ │ │ │ +
│ │ │ │ +
381 using CoordinateField = ctype;
│ │ │ │ +
│ │ │ │ +
384 using Coordinate = Dune::FieldVector<ctype,dim>;
│ │ │ │ +
│ │ │ │ +
387 static constexpr int dimension = dim;
│ │ │ │ +
│ │ │ │ +
390 typedef ctype Volume;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
394 friend class Impl::ReferenceElementContainer< ctype, dim >;
│ │ │ │ +
│ │ │ │ +
396 struct SubEntityInfo;
│ │ │ │ +
│ │ │ │ +
398 template<
int codim >
struct CreateGeometries;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
402 template<
int codim >
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
406 typedef AffineGeometry< ctype, dim-codim, dim > Geometry;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
410 ReferenceElementImplementation (
const ReferenceElementImplementation& ) =
delete;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
413 ReferenceElementImplementation& operator= (
const ReferenceElementImplementation& ) =
delete;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
416 ReferenceElementImplementation () =
default;
│ │ │ │ +
│ │ │ │ +
422 int size (
int c )
const
│ │ │ │ +
│ │ │ │ +
424 assert( (c >= 0) && (c <= dim) );
│ │ │ │ +
425 return info_[ c ].size();
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
439 int size (
int i,
int c,
int cc )
const
│ │ │ │ +
│ │ │ │ +
441 assert( (i >= 0) && (i <
size( c )) );
│ │ │ │ +
442 return info_[ c ][ i ].size( cc );
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
458 int subEntity (
int i,
int c,
int ii,
int cc )
const
│ │ │ │ +
│ │ │ │ +
460 assert( (i >= 0) && (i <
size( c )) );
│ │ │ │ +
461 return info_[ c ][ i ].number( ii, cc );
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
479 auto subEntities (
int i,
int c,
int cc )
const
│ │ │ │ +
│ │ │ │ +
481 assert( (i >= 0) && (i <
size( c )) );
│ │ │ │ +
482 return info_[ c ][ i ].numbers( cc );
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
493 const GeometryType &type (
int i,
int c )
const
│ │ │ │ +
│ │ │ │ +
495 assert( (i >= 0) && (i <
size( c )) );
│ │ │ │ +
496 return info_[ c ][ i ].type();
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
500 const GeometryType &type ()
const {
return type( 0, 0 ); }
│ │ │ │ +
│ │ │ │ +
511 const Coordinate &position(
int i,
int c )
const
│ │ │ │ +
│ │ │ │ +
513 assert( (c >= 0) && (c <= dim) );
│ │ │ │ +
514 return baryCenters_[ c ][ i ];
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
524 bool checkInside (
const Coordinate &local )
const
│ │ │ │ +
│ │ │ │ +
526 const ctype tolerance = ctype( 64 ) * std::numeric_limits< ctype >::epsilon();
│ │ │ │ +
527 return Impl::template checkInside< ctype, dim >( type().
id(), dim, local, tolerance );
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
541 template<
int codim >
│ │ │ │ +
542 typename Codim< codim >::Geometry geometry (
int i )
const
│ │ │ │ +
│ │ │ │ +
544 return std::get< codim >( geometries_ )[ i ];
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
548 Volume volume ()
const
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
560 const Coordinate &integrationOuterNormal (
int face )
const
│ │ │ │ +
│ │ │ │ +
562 assert( (face >= 0) && (face <
int( integrationNormals_.size() )) );
│ │ │ │ +
563 return integrationNormals_[ face ];
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
567 void initialize (
unsigned int topologyId )
│ │ │ │ +
│ │ │ │ +
569 assert( topologyId < Impl::numTopologies( dim ) );
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
572 for(
int codim = 0; codim <= dim; ++codim )
│ │ │ │ +
│ │ │ │ +
574 const unsigned int size = Impl::size( topologyId, dim, codim );
│ │ │ │ +
575 info_[ codim ].resize( size );
│ │ │ │ +
576 for(
unsigned int i = 0; i <
size; ++i )
│ │ │ │ +
577 info_[ codim ][ i ].initialize( topologyId, codim, i );
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
581 const unsigned int numVertices =
size( dim );
│ │ │ │ +
582 baryCenters_[ dim ].resize( numVertices );
│ │ │ │ +
583 Impl::referenceCorners( topologyId, dim, &(baryCenters_[ dim ][ 0 ]) );
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
586 for(
int codim = 0; codim < dim; ++codim )
│ │ │ │ +
│ │ │ │ +
588 baryCenters_[ codim ].resize(
size(codim) );
│ │ │ │ +
589 for(
int i = 0; i <
size( codim ); ++i )
│ │ │ │ +
│ │ │ │ +
591 baryCenters_[ codim ][ i ] = Coordinate( ctype( 0 ) );
│ │ │ │ +
592 const unsigned int numCorners =
size( i, codim, dim );
│ │ │ │ +
593 for(
unsigned int j = 0; j < numCorners; ++j )
│ │ │ │ +
594 baryCenters_[ codim ][ i ] += baryCenters_[ dim ][ subEntity( i, codim, j, dim ) ];
│ │ │ │ +
595 baryCenters_[ codim ][ i ] *= ctype( 1 ) / ctype( numCorners );
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
600 volume_ = Impl::template referenceVolume< ctype >( topologyId, dim );
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
605 integrationNormals_.resize(
size( 1 ) );
│ │ │ │ +
606 Impl::referenceIntegrationOuterNormals( topologyId, dim, &(integrationNormals_[ 0 ]) );
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
610 Hybrid::forEach( std::make_index_sequence< dim+1 >{}, [ & ](
auto i ){ CreateGeometries< i >::apply( *
this, geometries_ ); } );
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
613 template<
int... codim >
│ │ │ │ +
614 static std::tuple< std::vector< typename Codim< codim >::Geometry >... >
│ │ │ │ +
615 makeGeometryTable ( std::integer_sequence< int, codim... > );
│ │ │ │ +
│ │ │ │ +
618 typedef decltype( makeGeometryTable( std::make_integer_sequence< int, dim+1 >() ) ) GeometryTable;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
623 std::vector< Coordinate > baryCenters_[ dim+1 ];
│ │ │ │ +
624 std::vector< Coordinate > integrationNormals_;
│ │ │ │ +
│ │ │ │ +
627 GeometryTable geometries_;
│ │ │ │ +
│ │ │ │ +
629 std::vector< SubEntityInfo > info_[ dim+1 ];
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
633 template<
class ctype,
int dim >
│ │ │ │ +
634 struct ReferenceElementImplementation< ctype, dim >::SubEntityInfo
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
639 static constexpr std::size_t maxSubEntityCount()
│ │ │ │ +
│ │ │ │ +
641 std::size_t maxCount=0;
│ │ │ │ +
642 for(std::size_t codim=0; codim<=dim; ++codim)
│ │ │ │ +
643 maxCount = std::max(maxCount,
binomial(std::size_t(dim),codim)*(1 << codim));
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
647 using SubEntityFlags = std::bitset<maxSubEntityCount()>;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
650 :
public Dune::IteratorRange<const unsigned int*>
│ │ │ │ +
│ │ │ │ +
652 using Base =
typename Dune::IteratorRange<const unsigned int*>;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
656 using iterator = Base::iterator;
│ │ │ │ +
657 using const_iterator = Base::const_iterator;
│ │ │ │ +
│ │ │ │ +
659 SubEntityRange(
const iterator& begin,
const iterator& end,
const SubEntityFlags& contains) :
│ │ │ │ +
│ │ │ │ +
661 containsPtr_(&contains),
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
667 containsPtr_(nullptr),
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
671 std::size_t
size()
const
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
676 bool contains(std::size_t i)
const
│ │ │ │ +
│ │ │ │ +
678 return (*containsPtr_)[i];
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
682 const SubEntityFlags* containsPtr_;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
687 using NumberRange =
typename Dune::IteratorRange<const unsigned int*>;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
690 : numbering_( nullptr )
│ │ │ │ +
│ │ │ │ +
692 std::fill( offset_.begin(), offset_.end(), 0 );
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
695 SubEntityInfo (
const SubEntityInfo &other )
│ │ │ │ +
696 : offset_( other.offset_ ),
│ │ │ │ +
697 type_( other.type_ ),
│ │ │ │ +
698 containsSubentity_( other.containsSubentity_ )
│ │ │ │ +
│ │ │ │ +
700 numbering_ = allocate();
│ │ │ │ +
701 std::copy( other.numbering_, other.numbering_ + capacity(), numbering_ );
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
704 ~SubEntityInfo () { deallocate( numbering_ ); }
│ │ │ │ +
│ │ │ │ +
706 const SubEntityInfo &operator= (
const SubEntityInfo &other )
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
709 offset_ = other.offset_;
│ │ │ │ +
│ │ │ │ +
711 deallocate( numbering_ );
│ │ │ │ +
712 numbering_ = allocate();
│ │ │ │ +
713 std::copy( other.numbering_, other.numbering_ + capacity(), numbering_ );
│ │ │ │ +
│ │ │ │ +
715 containsSubentity_ = other.containsSubentity_;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
720 int size (
int cc )
const
│ │ │ │ +
│ │ │ │ +
722 assert( (cc >= 0) && (cc <= dim) );
│ │ │ │ +
723 return (offset_[ cc+1 ] - offset_[ cc ]);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
726 int number (
int ii,
int cc )
const
│ │ │ │ +
│ │ │ │ +
728 assert( (ii >= 0) && (ii <
size( cc )) );
│ │ │ │ +
729 return numbering_[ offset_[ cc ] + ii ];
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
732 auto numbers (
int cc )
const
│ │ │ │ +
│ │ │ │ +
734 return SubEntityRange( numbering_ + offset_[ cc ], numbering_ + offset_[ cc+1 ], containsSubentity_[cc]);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
737 const GeometryType &type ()
const {
return type_; }
│ │ │ │ +
│ │ │ │ +
739 void initialize (
unsigned int topologyId,
int codim,
unsigned int i )
│ │ │ │ +
│ │ │ │ +
741 const unsigned int subId = Impl::subTopologyId( topologyId, dim, codim, i );
│ │ │ │ +
742 type_ = GeometryType( subId, dim-codim );
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
745 for(
int cc = 0; cc <= codim; ++cc )
│ │ │ │ +
│ │ │ │ +
747 for(
int cc = codim; cc <= dim; ++cc )
│ │ │ │ +
748 offset_[ cc+1 ] = offset_[ cc ] + Impl::size( subId, dim-codim, cc-codim );
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
751 deallocate( numbering_ );
│ │ │ │ +
752 numbering_ = allocate();
│ │ │ │ +
753 for(
int cc = codim; cc <= dim; ++cc )
│ │ │ │ +
754 Impl::subTopologyNumbering( topologyId, dim, codim, i, cc-codim, numbering_+offset_[ cc ], numbering_+offset_[ cc+1 ] );
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
757 for(std::size_t cc=0; cc<= dim; ++cc)
│ │ │ │ +
│ │ │ │ +
759 containsSubentity_[cc].reset();
│ │ │ │ +
760 for(std::size_t idx=0; idx<std::size_t(
size(cc)); ++idx)
│ │ │ │ +
761 containsSubentity_[cc][number(idx,cc)] =
true;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
766 int codim ()
const {
return dim - type().dim(); }
│ │ │ │ +
│ │ │ │ +
768 unsigned int *allocate () {
return (capacity() != 0 ?
new unsigned int[ capacity() ] : nullptr); }
│ │ │ │ +
769 void deallocate (
unsigned int *ptr ) {
delete[] ptr; }
│ │ │ │ +
770 unsigned int capacity ()
const {
return offset_[ dim+1 ]; }
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
773 unsigned int *numbering_;
│ │ │ │ +
774 std::array< unsigned int, dim+2 > offset_;
│ │ │ │ +
│ │ │ │ +
776 std::array< SubEntityFlags, dim+1> containsSubentity_;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
780 template<
class ctype,
int dim >
│ │ │ │ +
781 template<
int codim >
│ │ │ │ +
782 struct ReferenceElementImplementation< ctype, dim >::CreateGeometries
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
785 static typename ReferenceElements< ctype, dim-cc >::ReferenceElement
│ │ │ │ +
786 subRefElement(
const ReferenceElementImplementation< ctype, dim > &refElement,
int i, std::integral_constant< int, cc > )
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
792 subRefElement(
const ReferenceElementImplementation< ctype, dim > &refElement,
│ │ │ │ +
793 [[maybe_unused]]
int i, std::integral_constant<int, 0>)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
798 static void apply (
const ReferenceElementImplementation< ctype, dim > &refElement, GeometryTable &geometries )
│ │ │ │ +
│ │ │ │ +
800 const int size = refElement.size( codim );
│ │ │ │ +
801 std::vector< FieldVector< ctype, dim > > origins( size );
│ │ │ │ +
802 std::vector< FieldMatrix< ctype, dim - codim, dim > > jacobianTransposeds( size );
│ │ │ │ +
803 Impl::referenceEmbeddings( refElement.type().id(), dim, codim, &(origins[ 0 ]), &(jacobianTransposeds[ 0 ]) );
│ │ │ │ +
│ │ │ │ +
805 std::get< codim >( geometries ).reserve( size );
│ │ │ │ +
806 for(
int i = 0; i <
size; ++i )
│ │ │ │ +
│ │ │ │ +
808 typename Codim< codim >::Geometry geometry( subRefElement( refElement, i, std::integral_constant< int, codim >() ), origins[ i ], jacobianTransposeds[ i ] );
│ │ │ │ +
809 std::get< codim >( geometries ).push_back( geometry );
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
A unique label for each type of element that can occur in a grid.
│ │ │ │ +
│ │ │ │ +
An implementation of the Geometry interface for affine geometries.
│ │ │ │ +
Definition affinegeometry.hh:21
│ │ │ │ +
@ size
Definition quadraturerules.hh:145
│ │ │ │ +
int binomial(int upper, int lower)
calculate
Definition simplex.cc:305
│ │ │ │ +
typename Container::ReferenceElement ReferenceElement
The reference element type.
Definition referenceelements.hh:188
│ │ │ │ +
static const ReferenceElement & general(const GeometryType &type)
get general reference elements
Definition referenceelements.hh:198
│ │ │ │
