│ │ │ │
│ │ │ │
│ │ │ │
│ │ │ │ -
5#ifndef DUNE_GEOMETRY_REFERENCEELEMENTIMPLEMENTATION_HH
│ │ │ │ -
6#define DUNE_GEOMETRY_REFERENCEELEMENTIMPLEMENTATION_HH
│ │ │ │ +
5#ifndef DUNE_GEOMETRY_MULTILINEARGEOMETRY_HH
│ │ │ │ +
6#define DUNE_GEOMETRY_MULTILINEARGEOMETRY_HH
│ │ │ │
│ │ │ │
│ │ │ │ -
│ │ │ │ -
│ │ │ │ +
│ │ │ │ +
│ │ │ │
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
14#include <dune/common/fmatrix.hh>
│ │ │ │ +
15#include <dune/common/fvector.hh>
│ │ │ │ +
16#include <dune/common/typetraits.hh>
│ │ │ │
│ │ │ │ -
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>
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
42 template<
class ctype,
int dim >
│ │ │ │ -
43 class ReferenceElementContainer;
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
46 template<
class ctype,
int dim >
│ │ │ │ -
47 struct ReferenceElements;
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
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;
│ │ │ │ -
│ │ │ │ -
62 unsigned int size (
unsigned int topologyId,
int dim,
int codim );
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
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 );
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
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 ) );
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
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 );
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
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 ) );
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
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 ) );
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
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 ) );
│ │ │ │ -
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
61 static ct
tolerance () {
return ct( 16 ) * std::numeric_limits< ct >::epsilon(); }
│ │ │ │ +
│ │ │ │ +
127 template<
int mydim,
int cdim >
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
130 typedef std::vector< FieldVector< ct, cdim > >
Type;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
149 static const bool v =
false;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
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 ) );
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ +
179 template<
class ct,
int mydim,
int cdim,
class Traits = MultiLinearGeometryTraits< ct > >
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
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 ) );
│ │ │ │ +
207 typedef FieldMatrix< ctype, coorddimension, mydimension >
Jacobian;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │
│ │ │ │ -
221 if( (0 < codim) && (codim <= 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 ) );
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
248 else if( codim < dim )
│ │ │ │ -
│ │ │ │ -
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 );
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
261 else if( codim == 0 )
│ │ │ │ -
│ │ │ │ -
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 );
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
280 template<
class ct,
int cdim >
│ │ │ │ -
│ │ │ │ -
282 referenceIntegrationOuterNormals (
unsigned int topologyId,
int dim,
│ │ │ │ -
283 const FieldVector< ct, cdim > *origins,
│ │ │ │ -
284 FieldVector< ct, cdim > *normals )
│ │ │ │ -
│ │ │ │ -
286 assert( (dim > 0) && (dim <= cdim) );
│ │ │ │ -
287 assert( topologyId < numTopologies( dim ) );
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
291 const unsigned int baseId = baseTopologyId( topologyId, dim );
│ │ │ │ -
292 if( isPrism( topologyId, dim ) )
│ │ │ │ -
│ │ │ │ -
294 const unsigned int numBaseFaces
│ │ │ │ -
295 = referenceIntegrationOuterNormals( baseId, dim-1, origins, normals );
│ │ │ │ -
│ │ │ │ -
297 for(
unsigned int i = 0; i < 2; ++i )
│ │ │ │ -
│ │ │ │ -
299 normals[ numBaseFaces+i ] = FieldVector< ct, cdim >( ct( 0 ) );
│ │ │ │ -
300 normals[ numBaseFaces+i ][ dim-1 ] = ct( 2*
int( i )-1 );
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
303 return numBaseFaces+2;
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
307 normals[ 0 ] = FieldVector< ct, cdim >( ct( 0 ) );
│ │ │ │ -
308 normals[ 0 ][ dim-1 ] = ct( -1 );
│ │ │ │ -
│ │ │ │ -
310 const unsigned int numBaseFaces
│ │ │ │ -
311 = referenceIntegrationOuterNormals( baseId, dim-1, origins+1, normals+1 );
│ │ │ │ -
312 for(
unsigned int i = 1; i <= numBaseFaces; ++i )
│ │ │ │ -
313 normals[ i ][ dim-1 ] = normals[ i ]*origins[ i ];
│ │ │ │ -
│ │ │ │ -
315 return numBaseFaces+1;
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
320 for(
unsigned int i = 0; i < 2; ++i )
│ │ │ │ -
│ │ │ │ -
322 normals[ i ] = FieldVector< ct, cdim >( ct( 0 ) );
│ │ │ │ -
323 normals[ i ][ 0 ] = ct( 2*
int( i )-1 );
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ +
│ │ │ │ +
222 static const bool hasSingleGeometryType = Traits::template hasSingleGeometryType< mydimension >::v;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
226 typedef typename std::conditional< hasSingleGeometryType, std::integral_constant< unsigned int, Traits::template hasSingleGeometryType< mydimension >::topologyId >,
unsigned int >
::type TopologyId;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
238 template<
class Corners >
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
254 template<
class Corners >
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
277 assert( (i >= 0) && (i <
corners()) );
│ │ │ │ +
278 return std::cref(corners_).get()[ i ];
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
294 auto cit = begin(std::cref(corners_).get());
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
314 const ctype tolerance = Traits::tolerance();
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
317 const bool affineMapping = this->
affine();
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
321 const GlobalCoordinate dglobal = (*this).global( x ) - globalCoord;
│ │ │ │ +
322 const bool invertible =
│ │ │ │ +
323 MatrixHelper::template xTRightInvA< mydimension, coorddimension >(
jacobianTransposed( x ), dglobal, dx );
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │
│ │ │ │ -
330 template<
class ct,
int cdim >
│ │ │ │ -
│ │ │ │ -
332 referenceIntegrationOuterNormals (
unsigned int topologyId,
int dim,
│ │ │ │ -
333 FieldVector< ct, cdim > *normals )
│ │ │ │ -
│ │ │ │ -
335 assert( (dim > 0) && (dim <= cdim) );
│ │ │ │ -
│ │ │ │ -
337 FieldVector< ct, cdim > *origins
│ │ │ │ -
338 =
new FieldVector< ct, cdim >[
size( topologyId, dim, 1 ) ];
│ │ │ │ -
339 referenceOrigins( topologyId, dim, 1, origins );
│ │ │ │ -
│ │ │ │ -
341 const unsigned int numFaces
│ │ │ │ -
342 = referenceIntegrationOuterNormals( topologyId, dim, origins, normals );
│ │ │ │ -
343 assert( numFaces ==
size( topologyId, dim, 1 ) );
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
375 template<
class ctype_,
int dim >
│ │ │ │ -
376 class ReferenceElementImplementation
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ +
│ │ │ │ +
331 if ( affineMapping )
break;
│ │ │ │ +
332 }
while( dx.two_norm2() > tolerance );
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │
│ │ │ │ -
382 using ctype = ctype_;
│ │ │ │ -
│ │ │ │ -
385 using CoordinateField = ctype;
│ │ │ │ +
│ │ │ │ +
382 auto cit = begin(std::cref(corners_).get());
│ │ │ │ +
383 jacobianTransposed< false >(
topologyId(), std::integral_constant< int, mydimension >(), cit,
ctype( 1 ),
local,
ctype( 1 ), jt );
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │
│ │ │ │ -
388 using Coordinate = Dune::FieldVector<ctype,dim>;
│ │ │ │ -
│ │ │ │ -
391 static constexpr int dimension = dim;
│ │ │ │ -
│ │ │ │ -
394 typedef ctype Volume;
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
398 friend class Impl::ReferenceElementContainer< ctype, dim >;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │
│ │ │ │ -
400 struct SubEntityInfo;
│ │ │ │ -
│ │ │ │ -
402 template<
int codim >
struct CreateGeometries;
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
406 template<
int codim >
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
410 typedef AffineGeometry< ctype, dim-codim, dim > Geometry;
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
414 ReferenceElementImplementation (
const ReferenceElementImplementation& ) =
delete;
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
417 ReferenceElementImplementation& operator= (
const ReferenceElementImplementation& ) =
delete;
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
420 ReferenceElementImplementation () =
default;
│ │ │ │ -
│ │ │ │ -
426 int size (
int c )
const
│ │ │ │ -
│ │ │ │ -
428 assert( (c >= 0) && (c <= dim) );
│ │ │ │ -
429 return info_[ c ].size();
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
443 int size (
int i,
int c,
int cc )
const
│ │ │ │ -
│ │ │ │ -
445 assert( (i >= 0) && (i <
size( c )) );
│ │ │ │ -
446 return info_[ c ][ i ].size( cc );
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
462 int subEntity (
int i,
int c,
int ii,
int cc )
const
│ │ │ │ -
│ │ │ │ -
464 assert( (i >= 0) && (i <
size( c )) );
│ │ │ │ -
465 return info_[ c ][ i ].number( ii, cc );
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
483 auto subEntities (
int i,
int c,
int cc )
const
│ │ │ │ -
│ │ │ │ -
485 assert( (i >= 0) && (i <
size( c )) );
│ │ │ │ -
486 return info_[ c ][ i ].numbers( cc );
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
497 const GeometryType &type (
int i,
int c )
const
│ │ │ │ -
│ │ │ │ -
499 assert( (i >= 0) && (i <
size( c )) );
│ │ │ │ -
500 return info_[ c ][ i ].type();
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
504 const GeometryType &type ()
const {
return type( 0, 0 ); }
│ │ │ │ -
│ │ │ │ -
515 const Coordinate &position(
int i,
int c )
const
│ │ │ │ -
│ │ │ │ -
517 assert( (c >= 0) && (c <= dim) );
│ │ │ │ -
518 return baryCenters_[ c ][ i ];
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
528 bool checkInside (
const Coordinate &local )
const
│ │ │ │ -
│ │ │ │ -
530 const ctype tolerance = ctype( 64 ) * std::numeric_limits< ctype >::epsilon();
│ │ │ │ -
531 return Impl::template checkInside< ctype, dim >( type().
id(), dim, local, tolerance );
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
545 template<
int codim >
│ │ │ │ -
546 typename Codim< codim >::Geometry geometry (
int i )
const
│ │ │ │ -
│ │ │ │ -
548 return std::get< codim >( geometries_ )[ i ];
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
552 Volume volume ()
const
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
564 const Coordinate &integrationOuterNormal (
int face )
const
│ │ │ │ -
│ │ │ │ -
566 assert( (face >= 0) && (face <
int( integrationNormals_.size() )) );
│ │ │ │ -
567 return integrationNormals_[ face ];
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
571 void initialize (
unsigned int topologyId )
│ │ │ │ -
│ │ │ │ -
573 assert( topologyId < Impl::numTopologies( dim ) );
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
576 for(
int codim = 0; codim <= dim; ++codim )
│ │ │ │ -
│ │ │ │ -
578 const unsigned int size = Impl::size( topologyId, dim, codim );
│ │ │ │ -
579 info_[ codim ].resize( size );
│ │ │ │ -
580 for(
unsigned int i = 0; i <
size; ++i )
│ │ │ │ -
581 info_[ codim ][ i ].initialize( topologyId, codim, i );
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
585 const unsigned int numVertices =
size( dim );
│ │ │ │ -
586 baryCenters_[ dim ].resize( numVertices );
│ │ │ │ -
587 Impl::referenceCorners( topologyId, dim, &(baryCenters_[ dim ][ 0 ]) );
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
590 for(
int codim = 0; codim < dim; ++codim )
│ │ │ │ -
│ │ │ │ -
592 baryCenters_[ codim ].resize(
size(codim) );
│ │ │ │ -
593 for(
int i = 0; i <
size( codim ); ++i )
│ │ │ │ -
│ │ │ │ -
595 baryCenters_[ codim ][ i ] = Coordinate( ctype( 0 ) );
│ │ │ │ -
596 const unsigned int numCorners =
size( i, codim, dim );
│ │ │ │ -
597 for(
unsigned int j = 0; j < numCorners; ++j )
│ │ │ │ -
598 baryCenters_[ codim ][ i ] += baryCenters_[ dim ][ subEntity( i, codim, j, dim ) ];
│ │ │ │ -
599 baryCenters_[ codim ][ i ] *= ctype( 1 ) / ctype( numCorners );
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
604 volume_ = Impl::template referenceVolume< ctype >( topologyId, dim );
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
609 integrationNormals_.resize(
size( 1 ) );
│ │ │ │ -
610 Impl::referenceIntegrationOuterNormals( topologyId, dim, &(integrationNormals_[ 0 ]) );
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
614 Hybrid::forEach( std::make_index_sequence< dim+1 >{}, [ & ](
auto i ){ CreateGeometries< i >::apply( *
this, geometries_ ); } );
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
617 template<
int... codim >
│ │ │ │ -
618 static std::tuple< std::vector< typename Codim< codim >::Geometry >... >
│ │ │ │ -
619 makeGeometryTable ( std::integer_sequence< int, codim... > );
│ │ │ │ -
│ │ │ │ -
622 typedef decltype( makeGeometryTable( std::make_integer_sequence< int, dim+1 >() ) ) GeometryTable;
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
627 std::vector< Coordinate > baryCenters_[ dim+1 ];
│ │ │ │ -
628 std::vector< Coordinate > integrationNormals_;
│ │ │ │ -
│ │ │ │ -
631 GeometryTable geometries_;
│ │ │ │ -
│ │ │ │ -
633 std::vector< SubEntityInfo > info_[ dim+1 ];
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
637 template<
class ctype,
int dim >
│ │ │ │ -
638 struct ReferenceElementImplementation< ctype, dim >::SubEntityInfo
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
643 static constexpr std::size_t maxSubEntityCount()
│ │ │ │ -
│ │ │ │ -
645 std::size_t maxCount=0;
│ │ │ │ -
646 for(std::size_t codim=0; codim<=dim; ++codim)
│ │ │ │ -
647 maxCount = std::max(maxCount, binomial(std::size_t(dim),codim)*(1 << codim));
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
651 using SubEntityFlags = std::bitset<maxSubEntityCount()>;
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
654 :
public Dune::IteratorRange<const unsigned int*>
│ │ │ │ -
│ │ │ │ -
656 using Base =
typename Dune::IteratorRange<const unsigned int*>;
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
660 using iterator = Base::iterator;
│ │ │ │ -
661 using const_iterator = Base::const_iterator;
│ │ │ │ -
│ │ │ │ -
663 SubEntityRange(
const iterator& begin,
const iterator& end,
const SubEntityFlags& contains) :
│ │ │ │ -
│ │ │ │ -
665 containsPtr_(&contains),
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
671 containsPtr_(nullptr),
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
675 std::size_t
size()
const
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
680 bool contains(std::size_t i)
const
│ │ │ │ -
│ │ │ │ -
682 return (*containsPtr_)[i];
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
686 const SubEntityFlags* containsPtr_;
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
691 using NumberRange =
typename Dune::IteratorRange<const unsigned int*>;
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
694 : numbering_( nullptr )
│ │ │ │ -
│ │ │ │ -
696 std::fill( offset_.begin(), offset_.end(), 0 );
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
699 SubEntityInfo (
const SubEntityInfo &other )
│ │ │ │ -
700 : offset_( other.offset_ ),
│ │ │ │ -
701 type_( other.type_ ),
│ │ │ │ -
702 containsSubentity_( other.containsSubentity_ )
│ │ │ │ -
│ │ │ │ -
704 numbering_ = allocate();
│ │ │ │ -
705 std::copy( other.numbering_, other.numbering_ + capacity(), numbering_ );
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
708 ~SubEntityInfo () { deallocate( numbering_ ); }
│ │ │ │ -
│ │ │ │ -
710 const SubEntityInfo &operator= (
const SubEntityInfo &other )
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
713 offset_ = other.offset_;
│ │ │ │ -
│ │ │ │ -
715 deallocate( numbering_ );
│ │ │ │ -
716 numbering_ = allocate();
│ │ │ │ -
717 std::copy( other.numbering_, other.numbering_ + capacity(), numbering_ );
│ │ │ │ -
│ │ │ │ -
719 containsSubentity_ = other.containsSubentity_;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
432 return topologyId( std::integral_constant< bool, hasSingleGeometryType >() );
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
435 template<
bool add,
int dim,
class CornerIterator >
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
439 template<
bool add,
class CornerIterator >
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
444 template<
bool add,
int rows,
int dim,
class CornerIterator >
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
447 const ctype &rf, FieldMatrix< ctype, rows, cdim > &jt );
│ │ │ │ +
448 template<
bool add,
int rows,
class CornerIterator >
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
451 const ctype &rf, FieldMatrix< ctype, rows, cdim > &jt );
│ │ │ │ +
│ │ │ │ +
453 template<
int dim,
class CornerIterator >
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
455 template<
class CornerIterator >
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
462 auto cit = begin(std::cref(corners_).get());
│ │ │ │ +
463 return affine(
topologyId(), std::integral_constant< int, mydimension >(), cit, jacobianT );
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
470 static unsigned int toUnsignedInt(
unsigned int i) {
return i; }
│ │ │ │ +
471 template<
unsigned int v>
│ │ │ │ +
472 static unsigned int toUnsignedInt(std::integral_constant<unsigned int,v> ) {
return v; }
│ │ │ │ +
│ │ │ │ +
474 unsigned int topologyId ( std::integral_constant< bool, false > )
const {
return refElement().type().id(); }
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
477 typename Traits::template CornerStorage< mydimension, coorddimension >::Type corners_;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
485 template<
class ct,
int mydim,
int cdim,
class Traits >
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
487 :
public FieldMatrix< ctype, coorddimension, mydimension >
│ │ │ │ +
│ │ │ │ +
489 typedef FieldMatrix< ctype, coorddimension, mydimension > Base;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
494 detInv_ = MatrixHelper::template rightInvA< mydimension, coorddimension >( jt,
static_cast< Base &
>( *
this ) );
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
499 detInv_ = MatrixHelper::template sqrtDetAAT< mydimension, coorddimension >( jt );
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
523 template<
class ct,
int mydim,
int cdim,
class Traits = MultiLinearGeometryTraits< ct > >
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
550 template<
class CornerStorage >
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
553 affine_(
Base::
affine( jacobianTransposed_ ) ),
│ │ │ │ +
554 jacobianInverseTransposedComputed_( false ),
│ │ │ │ +
555 integrationElementComputed_( false )
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
558 template<
class CornerStorage >
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
560 :
Base( gt, cornerStorage ),
│ │ │ │ +
561 affine_(
Base::
affine( jacobianTransposed_ ) ),
│ │ │ │ +
562 jacobianInverseTransposedComputed_( false ),
│ │ │ │ +
563 integrationElementComputed_( false )
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
609 if( jacobianInverseTransposedComputed_ )
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
612 MatrixHelper::template xTRightInvA< mydimension, coorddimension >( jacobianTransposed_,
global -
corner( 0 ),
local );
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
637 if( !integrationElementComputed_ )
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
640 integrationElementComputed_ =
true;
│ │ │ │ +
│ │ │ │ +
642 return jacobianInverseTransposed_.
detInv();
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
669 return jacobianTransposed_;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
684 if( !jacobianInverseTransposedComputed_ )
│ │ │ │ +
│ │ │ │ +
686 jacobianInverseTransposed_.
setup( jacobianTransposed_ );
│ │ │ │ +
687 jacobianInverseTransposedComputed_ =
true;
│ │ │ │ +
688 integrationElementComputed_ =
true;
│ │ │ │ +
│ │ │ │ +
690 return jacobianInverseTransposed_;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
724 int size (
int cc )
const
│ │ │ │ -
│ │ │ │ -
726 assert( (cc >= 0) && (cc <= dim) );
│ │ │ │ -
727 return (offset_[ cc+1 ] - offset_[ cc ]);
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
730 int number (
int ii,
int cc )
const
│ │ │ │ -
│ │ │ │ -
732 assert( (ii >= 0) && (ii <
size( cc )) );
│ │ │ │ -
733 return numbering_[ offset_[ cc ] + ii ];
│ │ │ │ -
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
725 mutable bool affine_ : 1;
│ │ │ │ +
│ │ │ │ +
727 mutable bool jacobianInverseTransposedComputed_ : 1;
│ │ │ │ +
728 mutable bool integrationElementComputed_ : 1;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │
│ │ │ │ -
736 auto numbers (
int cc )
const
│ │ │ │ -
│ │ │ │ -
738 return SubEntityRange( numbering_ + offset_[ cc ], numbering_ + offset_[ cc+1 ], containsSubentity_[cc]);
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
741 const GeometryType &type ()
const {
return type_; }
│ │ │ │ -
│ │ │ │ -
743 void initialize (
unsigned int topologyId,
int codim,
unsigned int i )
│ │ │ │ -
│ │ │ │ -
745 const unsigned int subId = Impl::subTopologyId( topologyId, dim, codim, i );
│ │ │ │ -
746 type_ = GeometryType( subId, dim-codim );
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
749 for(
int cc = 0; cc <= codim; ++cc )
│ │ │ │ -
│ │ │ │ -
751 for(
int cc = codim; cc <= dim; ++cc )
│ │ │ │ -
752 offset_[ cc+1 ] = offset_[ cc ] + Impl::size( subId, dim-codim, cc-codim );
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
755 deallocate( numbering_ );
│ │ │ │ -
756 numbering_ = allocate();
│ │ │ │ -
757 for(
int cc = codim; cc <= dim; ++cc )
│ │ │ │ -
758 Impl::subTopologyNumbering( topologyId, dim, codim, i, cc-codim, numbering_+offset_[ cc ], numbering_+offset_[ cc+1 ] );
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
761 for(std::size_t cc=0; cc<= dim; ++cc)
│ │ │ │ -
│ │ │ │ -
763 containsSubentity_[cc].reset();
│ │ │ │ -
764 for(std::size_t idx=0; idx<std::size_t(
size(cc)); ++idx)
│ │ │ │ -
765 containsSubentity_[cc][number(idx,cc)] =
true;
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
770 int codim ()
const {
return dim - type().dim(); }
│ │ │ │ -
│ │ │ │ -
772 unsigned int *allocate () {
return (capacity() != 0 ?
new unsigned int[ capacity() ] : nullptr); }
│ │ │ │ -
773 void deallocate (
unsigned int *ptr ) {
delete[] ptr; }
│ │ │ │ -
774 unsigned int capacity ()
const {
return offset_[ dim+1 ]; }
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
777 unsigned int *numbering_;
│ │ │ │ -
778 std::array< unsigned int, dim+2 > offset_;
│ │ │ │ -
│ │ │ │ -
780 std::array< SubEntityFlags, dim+1> containsSubentity_;
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
784 template<
class ctype,
int dim >
│ │ │ │ -
785 template<
int codim >
│ │ │ │ -
786 struct ReferenceElementImplementation< ctype, dim >::CreateGeometries
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
789 static typename ReferenceElements< ctype, dim-cc >::ReferenceElement
│ │ │ │ -
790 subRefElement(
const ReferenceElementImplementation< ctype, dim > &refElement,
int i, std::integral_constant< int, cc > )
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
796 subRefElement(
const ReferenceElementImplementation< ctype, dim > &refElement,
│ │ │ │ -
797 [[maybe_unused]]
int i, std::integral_constant<int, 0>)
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
802 static void apply (
const ReferenceElementImplementation< ctype, dim > &refElement, GeometryTable &geometries )
│ │ │ │ -
│ │ │ │ -
804 const int size = refElement.size( codim );
│ │ │ │ -
805 std::vector< FieldVector< ctype, dim > > origins( size );
│ │ │ │ -
806 std::vector< FieldMatrix< ctype, dim - codim, dim > > jacobianTransposeds( size );
│ │ │ │ -
807 Impl::referenceEmbeddings( refElement.type().id(), dim, codim, &(origins[ 0 ]), &(jacobianTransposeds[ 0 ]) );
│ │ │ │ -
│ │ │ │ -
809 std::get< codim >( geometries ).reserve( size );
│ │ │ │ -
810 for(
int i = 0; i <
size; ++i )
│ │ │ │ -
│ │ │ │ -
812 typename Codim< codim >::Geometry geometry( subRefElement( refElement, i, std::integral_constant< int, codim >() ), origins[ i ], jacobianTransposeds[ i ] );
│ │ │ │ -
813 std::get< codim >( geometries ).push_back( geometry );
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
An implementation of the Geometry interface for affine geometries.
│ │ │ │ -
A unique label for each type of element that can occur in a grid.
│ │ │ │ -
│ │ │ │ +
736 template<
class ct,
int mydim,
int cdim,
class Traits >
│ │ │ │ +
737 inline typename MultiLinearGeometry< ct, mydim, cdim, Traits >::JacobianInverseTransposed
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
741 jit.
setup( jacobianTransposed( local ) );
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
746 template<
class ct,
int mydim,
int cdim,
class Traits >
│ │ │ │ +
747 template<
bool add,
int dim,
class CornerIterator >
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
753 const ctype xn = df*x[ dim-1 ];
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
756 if( Impl::isPrism( toUnsignedInt(topologyId), mydimension, mydimension-dim ) )
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
759 global< add >( topologyId, std::integral_constant< int, dim-1 >(), cit, df, x, rf*cxn, y );
│ │ │ │ +
│ │ │ │ +
761 global< true >( topologyId, std::integral_constant< int, dim-1 >(), cit, df, x, rf*xn, y );
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
765 assert( Impl::isPyramid( toUnsignedInt(topologyId), mydimension, mydimension-dim ) );
│ │ │ │ +
│ │ │ │ +
767 if( cxn > Traits::tolerance() || cxn < -Traits::tolerance() )
│ │ │ │ +
768 global< add >( topologyId, std::integral_constant< int, dim-1 >(), cit, df/cxn, x, rf*cxn, y );
│ │ │ │ +
│ │ │ │ +
770 global< add >( topologyId, std::integral_constant< int, dim-1 >(), cit, df, x,
ctype( 0 ), y );
│ │ │ │ +
│ │ │ │ +
772 y.axpy( rf*xn, *cit );
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
777 template<
class ct,
int mydim,
int cdim,
class Traits >
│ │ │ │ +
778 template<
bool add,
class CornerIterator >
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
786 for(
int i = 0; i < coorddimension; ++i )
│ │ │ │ +
787 y[ i ] = (add ? y[ i ] + rf*origin[ i ] : rf*origin[ i ]);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
791 template<
class ct,
int mydim,
int cdim,
class Traits >
│ │ │ │ +
792 template<
bool add,
int rows,
int dim,
class CornerIterator >
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
796 const ctype &rf, FieldMatrix< ctype, rows, cdim > &jt )
│ │ │ │ +
│ │ │ │ +
798 assert( rows >= dim );
│ │ │ │ +
│ │ │ │ +
800 const ctype xn = df*x[ dim-1 ];
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
804 if( Impl::isPrism( toUnsignedInt(topologyId), mydimension, mydimension-dim ) )
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
807 jacobianTransposed< add >( topologyId, std::integral_constant< int, dim-1 >(), cit2, df, x, rf*cxn, jt );
│ │ │ │ +
│ │ │ │ +
809 jacobianTransposed< true >( topologyId, std::integral_constant< int, dim-1 >(), cit2, df, x, rf*xn, jt );
│ │ │ │ +
│ │ │ │ +
811 global< add >( topologyId, std::integral_constant< int, dim-1 >(), cit, df, x, -rf, jt[ dim-1 ] );
│ │ │ │ +
812 global< true >( topologyId, std::integral_constant< int, dim-1 >(), cit, df, x, rf, jt[ dim-1 ] );
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
816 assert( Impl::isPyramid( toUnsignedInt(topologyId), mydimension, mydimension-dim ) );
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
858 ctype dfcxn = (cxn > Traits::tolerance() || cxn < -Traits::tolerance()) ?
ctype(df / cxn) :
ctype(0);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
863 global< add >( topologyId, std::integral_constant< int, dim-1 >(), cit, dfcxn, x, -rf, jt[ dim-1 ] );
│ │ │ │ +
│ │ │ │ +
865 jt[ dim-1 ].axpy( rf, *cit );
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
870 FieldMatrix<
ctype, dim-1, coorddimension > jt2;
│ │ │ │ +
│ │ │ │ +
872 jacobianTransposed< false >( topologyId, std::integral_constant< int, dim-1 >(), cit2, dfcxn, x, rf, jt2 );
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
876 for(
int j = 0; j < dim-1; ++j )
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
879 jt[ dim-1 ].axpy( dfcxn*x[ j ], jt2[ j ] );
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
885 jacobianTransposed< false >( topologyId, std::integral_constant< int, dim-1 >(), cit2, dfcxn, x, rf, jt );
│ │ │ │ +
│ │ │ │ +
887 for(
int j = 0; j < dim-1; ++j )
│ │ │ │ +
888 jt[ dim-1 ].axpy( dfcxn*x[ j ], jt[ j ] );
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
893 template<
class ct,
int mydim,
int cdim,
class Traits >
│ │ │ │ +
894 template<
bool add,
int rows,
class CornerIterator >
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
898 const ctype &, FieldMatrix< ctype, rows, cdim > & )
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
905 template<
class ct,
int mydim,
int cdim,
class Traits >
│ │ │ │ +
906 template<
int dim,
class CornerIterator >
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
911 if( !affine( topologyId, std::integral_constant< int, dim-1 >(), cit, jt ) )
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
915 if( Impl::isPrism( toUnsignedInt(topologyId), mydimension, mydimension-dim ) )
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
918 if( !affine( topologyId, std::integral_constant< int, dim-1 >(), cit, jtTop ) )
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
923 for(
int i = 0; i < dim-1; ++i )
│ │ │ │ +
924 norm += (jtTop[ i ] - jt[ i ]).two_norm2();
│ │ │ │ +
925 if( norm >= Traits::tolerance() )
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
930 jt[ dim-1 ] = orgTop - orgBottom;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
934 template<
class ct,
int mydim,
int cdim,
class Traits >
│ │ │ │ +
935 template<
class CornerIterator >
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
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:194
│ │ │ │ +
Class providing access to the singletons of the reference elements.
Definition referenceelements.hh:128
│ │ │ │
typename Container::ReferenceElement ReferenceElement
The reference element type.
Definition referenceelements.hh:146
│ │ │ │ -
static const ReferenceElement & general(const GeometryType &type)
get general reference elements
Definition referenceelements.hh:156
│ │ │ │ +
default traits class for MultiLinearGeometry
Definition multilineargeometry.hh:39
│ │ │ │ +
Impl::FieldMatrixHelper< ct > MatrixHelper
helper structure containing some matrix routines
Definition multilineargeometry.hh:58
│ │ │ │ +
static ct tolerance()
tolerance to numerical algorithms
Definition multilineargeometry.hh:61
│ │ │ │ +
template specifying the storage for the corners
Definition multilineargeometry.hh:129
│ │ │ │ +
std::vector< FieldVector< ct, cdim > > Type
Definition multilineargeometry.hh:130
│ │ │ │ +
will there be only one geometry type for a dimension?
Definition multilineargeometry.hh:148
│ │ │ │ +
static const unsigned int topologyId
Definition multilineargeometry.hh:150
│ │ │ │ +
static const bool v
Definition multilineargeometry.hh:149
│ │ │ │ +
generic geometry implementation based on corner coordinates
Definition multilineargeometry.hh:181
│ │ │ │ +
static void global(TopologyId topologyId, std::integral_constant< int, dim >, CornerIterator &cit, const ctype &df, const LocalCoordinate &x, const ctype &rf, GlobalCoordinate &y)
Definition multilineargeometry.hh:749
│ │ │ │ +
static const int mydimension
geometry dimension
Definition multilineargeometry.hh:189
│ │ │ │ +
Dune::GeometryType type() const
obtain the name of the reference element
Definition multilineargeometry.hh:269
│ │ │ │ +
Traits::MatrixHelper MatrixHelper
Definition multilineargeometry.hh:225
│ │ │ │ +
FieldVector< ctype, coorddimension > GlobalCoordinate
type of global coordinates
Definition multilineargeometry.hh:196
│ │ │ │ +
static void jacobianTransposed(TopologyId topologyId, std::integral_constant< int, 0 >, CornerIterator &cit, const ctype &df, const LocalCoordinate &x, const ctype &rf, FieldMatrix< ctype, rows, cdim > &jt)
Definition multilineargeometry.hh:896
│ │ │ │ +
JacobianTransposed jacobianTransposed(const LocalCoordinate &local) const
obtain the transposed of the Jacobian
Definition multilineargeometry.hh:377
│ │ │ │ +
ReferenceElement refElement() const
Definition multilineargeometry.hh:425
│ │ │ │ +
GlobalCoordinate global(const LocalCoordinate &local) const
evaluate the mapping
Definition multilineargeometry.hh:290
│ │ │ │ +
GlobalCoordinate center() const
obtain the centroid of the mapping's image
Definition multilineargeometry.hh:282
│ │ │ │ +
GlobalCoordinate corner(int i) const
obtain coordinates of the i-th corner
Definition multilineargeometry.hh:275
│ │ │ │ +
Dune::ReferenceElements< ctype, mydimension > ReferenceElements
Definition multilineargeometry.hh:214
│ │ │ │ +
ct ctype
coordinate type
Definition multilineargeometry.hh:186
│ │ │ │ +
static void jacobianTransposed(TopologyId topologyId, std::integral_constant< int, dim >, CornerIterator &cit, const ctype &df, const LocalCoordinate &x, const ctype &rf, FieldMatrix< ctype, rows, cdim > &jt)
Definition multilineargeometry.hh:794
│ │ │ │ +
static const int coorddimension
coordinate dimension
Definition multilineargeometry.hh:191
│ │ │ │ +
int corners() const
obtain number of corners of the corresponding reference element
Definition multilineargeometry.hh:272
│ │ │ │ +
TopologyId topologyId() const
Definition multilineargeometry.hh:430
│ │ │ │ +
friend ReferenceElement referenceElement(const MultiLinearGeometry &geometry)
Definition multilineargeometry.hh:395
│ │ │ │ +
LocalCoordinate local(const GlobalCoordinate &globalCoord) const
evaluate the inverse mapping
Definition multilineargeometry.hh:312
│ │ │ │ +
static void global(TopologyId topologyId, std::integral_constant< int, 0 >, CornerIterator &cit, const ctype &df, const LocalCoordinate &x, const ctype &rf, GlobalCoordinate &y)
Definition multilineargeometry.hh:780
│ │ │ │ +
Volume volume() const
obtain the volume of the mapping's image
Definition multilineargeometry.hh:363
│ │ │ │ +
FieldVector< ctype, mydimension > LocalCoordinate
type of local coordinates
Definition multilineargeometry.hh:194
│ │ │ │ +
MultiLinearGeometry(const ReferenceElement &refElement, const Corners &corners)
constructor
Definition multilineargeometry.hh:239
│ │ │ │ +
static bool affine(TopologyId topologyId, std::integral_constant< int, 0 >, CornerIterator &cit, JacobianTransposed &jt)
Definition multilineargeometry.hh:937
│ │ │ │ +
std::conditional< hasSingleGeometryType, std::integral_constant< unsignedint, Traits::templatehasSingleGeometryType< mydimension >::topologyId >, unsignedint >::type TopologyId
Definition multilineargeometry.hh:226
│ │ │ │ +
ctype Volume
type of volume
Definition multilineargeometry.hh:198
│ │ │ │ +
static bool affine(TopologyId topologyId, std::integral_constant< int, dim >, CornerIterator &cit, JacobianTransposed &jt)
Definition multilineargeometry.hh:908
│ │ │ │ +
JacobianInverse jacobianInverse(const LocalCoordinate &local) const
Obtain the Jacobian's inverse.
Definition multilineargeometry.hh:418
│ │ │ │ +
MultiLinearGeometry(Dune::GeometryType gt, const Corners &corners)
constructor
Definition multilineargeometry.hh:255
│ │ │ │ +
FieldMatrix< ctype, mydimension, coorddimension > JacobianTransposed
type of jacobian transposed
Definition multilineargeometry.hh:201
│ │ │ │ +
ReferenceElements::ReferenceElement ReferenceElement
type of reference element
Definition multilineargeometry.hh:219
│ │ │ │ +
JacobianInverseTransposed jacobianInverseTransposed(const LocalCoordinate &local) const
obtain the transposed of the Jacobian's inverse
Definition multilineargeometry.hh:738
│ │ │ │ +
FieldMatrix< ctype, coorddimension, mydimension > Jacobian
Type for the Jacobian matrix.
Definition multilineargeometry.hh:207
│ │ │ │ +
bool affine() const
is this mapping affine?
Definition multilineargeometry.hh:262
│ │ │ │ +
FieldMatrix< ctype, mydimension, coorddimension > JacobianInverse
Type for the inverse Jacobian matrix.
Definition multilineargeometry.hh:210
│ │ │ │ +
bool affine(JacobianTransposed &jacobianT) const
Definition multilineargeometry.hh:458
│ │ │ │ +
Volume integrationElement(const LocalCoordinate &local) const
obtain the integration element
Definition multilineargeometry.hh:350
│ │ │ │ +
Jacobian jacobian(const LocalCoordinate &local) const
Obtain the Jacobian.
Definition multilineargeometry.hh:407
│ │ │ │ +
Definition multilineargeometry.hh:488
│ │ │ │ +
void setup(const JacobianTransposed &jt)
Definition multilineargeometry.hh:492
│ │ │ │ +
ctype det() const
Definition multilineargeometry.hh:502
│ │ │ │ +
ctype detInv() const
Definition multilineargeometry.hh:503
│ │ │ │ +
void setupDeterminant(const JacobianTransposed &jt)
Definition multilineargeometry.hh:497
│ │ │ │ +
Implement a MultiLinearGeometry with additional caching.
Definition multilineargeometry.hh:526
│ │ │ │ +
GlobalCoordinate global(const LocalCoordinate &local) const
evaluate the mapping
Definition multilineargeometry.hh:580
│ │ │ │ +
Base::ReferenceElement ReferenceElement
Definition multilineargeometry.hh:534
│ │ │ │ +
bool affine() const
is this mapping affine?
Definition multilineargeometry.hh:567
│ │ │ │ +
CachedMultiLinearGeometry(const ReferenceElement &referenceElement, const CornerStorage &cornerStorage)
Definition multilineargeometry.hh:551
│ │ │ │ +
ReferenceElement refElement() const
Definition multilineargeometry.hh:425
│ │ │ │ +
LocalCoordinate local(const GlobalCoordinate &global) const
evaluate the inverse mapping
Definition multilineargeometry.hh:604
│ │ │ │ +
Base::MatrixHelper MatrixHelper
Definition multilineargeometry.hh:531
│ │ │ │ +
Base::LocalCoordinate LocalCoordinate
Definition multilineargeometry.hh:541
│ │ │ │ +
JacobianInverse jacobianInverse(const LocalCoordinate &local) const
Obtain the Jacobian's inverse.
Definition multilineargeometry.hh:713
│ │ │ │ +
JacobianTransposed jacobianTransposed(const LocalCoordinate &local) const
obtain the transposed of the Jacobian
Definition multilineargeometry.hh:666
│ │ │ │ +
GlobalCoordinate corner(int i) const
obtain coordinates of the i-th corner
Definition multilineargeometry.hh:275
│ │ │ │ +
Volume volume() const
obtain the volume of the mapping's image
Definition multilineargeometry.hh:649
│ │ │ │ +
CachedMultiLinearGeometry(Dune::GeometryType gt, const CornerStorage &cornerStorage)
Definition multilineargeometry.hh:559
│ │ │ │ +
ctype integrationElement(const LocalCoordinate &local) const
obtain the integration element
Definition multilineargeometry.hh:633
│ │ │ │ +
Base::ctype ctype
Definition multilineargeometry.hh:536
│ │ │ │ +
Jacobian jacobian(const LocalCoordinate &local) const
Obtain the Jacobian.
Definition multilineargeometry.hh:702
│ │ │ │ +
Base::JacobianInverseTransposed JacobianInverseTransposed
Definition multilineargeometry.hh:546
│ │ │ │ +
Base::JacobianTransposed JacobianTransposed
Definition multilineargeometry.hh:545
│ │ │ │ +
Base::JacobianInverse JacobianInverse
Definition multilineargeometry.hh:548
│ │ │ │ +
Base::Jacobian Jacobian
Definition multilineargeometry.hh:547
│ │ │ │ +
Base::Volume Volume
Definition multilineargeometry.hh:543
│ │ │ │ +
Base::GlobalCoordinate GlobalCoordinate
Definition multilineargeometry.hh:542
│ │ │ │ +
GlobalCoordinate center() const
obtain the centroid of the mapping's image
Definition multilineargeometry.hh:572
│ │ │ │ +
JacobianInverseTransposed jacobianInverseTransposed(const LocalCoordinate &local) const
obtain the transposed of the Jacobian's inverse
Definition multilineargeometry.hh:680
│ │ │ │ +
Unique label for each type of entities that can occur in DUNE grids.
Definition type.hh:114
│ │ │ │
