3#ifndef DUNE_GEOMETRY_GENERICGEOMETRY_REFERENCEDOMAIN_HH
4#define DUNE_GEOMETRY_GENERICGEOMETRY_REFERENCEDOMAIN_HH
14#include <dune/geometry/genericgeometry/topologytypes.hh>
15#include <dune/geometry/genericgeometry/subtopologies.hh>
20 namespace GenericGeometry
26 template<
class Topology >
27 struct ReferenceDomain;
34 template<
class Topology >
35 class ReferenceDomainBase;
39 class ReferenceDomainBase< Point >
41 typedef Point Topology;
43 friend struct ReferenceDomain< Topology >;
44 friend class ReferenceDomainBase< Prism< Topology > >;
45 friend class ReferenceDomainBase< Pyramid< Topology > >;
47 static const unsigned int numNormals = 0;
49 template<
class ctype,
int dim >
54 assert( i < Topology::numCorners );
57 template<
class ctype,
int dim >
66 template<
class ctype,
int dim >
72 assert( i < numNormals );
75 template<
class ctype >
76 static ctype volume ()
83 template<
class BaseTopology >
84 class ReferenceDomainBase< Prism< BaseTopology > >
86 typedef Prism< BaseTopology > Topology;
88 friend struct ReferenceDomain< Topology >;
89 friend class ReferenceDomainBase< Prism< Topology > >;
90 friend class ReferenceDomainBase< Pyramid< Topology > >;
92 static const unsigned int numNormals = Size< Topology, 1 >::value;
94 static const unsigned int dimension = Topology::dimension;
95 static const unsigned int myindex = dimension - 1;
97 template<
class ctype,
int dim >
100 assert( i < Topology::numCorners );
101 const unsigned int j = i % BaseTopology::numCorners;
102 ReferenceDomainBase< BaseTopology >::corner( j, x );
103 if( i >= BaseTopology::numCorners )
104 x[ myindex ] = ctype( 1 );
107 template<
class ctype,
int dim >
111 const ctype xn = x[ myindex ];
112 const ctype cxn = factor - xn;
113 return (xn > -1e-12) && (cxn > -1e-12)
114 && ReferenceDomainBase< BaseTopology >::checkInside( x, factor );
117 template<
class ctype,
int dim >
121 typedef ReferenceDomainBase< BaseTopology > BaseReferenceDomain;
123 if( i >= BaseReferenceDomain::numNormals )
125 const unsigned int j = i - BaseReferenceDomain::numNormals;
126 n[ myindex ] = (j == 0 ? ctype( -1 ) : ctype( 1 ));
129 BaseReferenceDomain::integrationOuterNormal( i, n );
132 template<
class ctype >
133 static ctype volume ()
135 typedef ReferenceDomainBase< BaseTopology > BaseReferenceDomain;
136 return BaseReferenceDomain::template volume< ctype >();
141 template<
class BaseTopology >
142 class ReferenceDomainBase< Pyramid< BaseTopology > >
144 typedef Pyramid< BaseTopology > Topology;
146 friend struct ReferenceDomain< Topology >;
147 friend class ReferenceDomainBase< Prism< Topology > >;
148 friend class ReferenceDomainBase< Pyramid< Topology > >;
150 static const unsigned int numNormals = Size< Topology, 1 >::value;
152 static const unsigned int dimension = Topology::dimension;
153 static const unsigned int myindex = dimension - 1;
156 struct MultiDimensional
158 template<
class ctype,
int dim >
162 multiDimensionalIntegrationOuterNormal( i, n );
167 struct OneDimensional
169 template<
class ctype,
int dim >
173 n[ myindex ] = (i > 0) ? ctype( 1 ) : ctype( -1 );
177 template<
class ctype,
int dim >
180 assert( i < Topology::numCorners );
181 if( i < BaseTopology::numCorners )
182 ReferenceDomainBase< BaseTopology >::corner( i, x );
184 x[ myindex ] = ctype( 1 );
187 template<
class ctype,
int dim >
191 const ctype xn = x[ myindex ];
192 const ctype cxn = factor - xn;
193 return (xn > -1e-12) && (cxn > -1e-12)
194 && ReferenceDomainBase< BaseTopology >::checkInside( x, cxn );
197 template<
class ctype,
int dim >
201 typedef ReferenceDomainBase< BaseTopology > BaseReferenceDomain;
202 typedef SubTopologyNumbering< BaseTopology, 1, dimension-2 > Numbering;
206 const unsigned int j = Numbering::number( i-1, 0 );
208 BaseReferenceDomain::corner( j, x );
210 BaseReferenceDomain::integrationOuterNormal ( i-1, n );
211 n[ myindex ] = (x * n);
214 n[ myindex ] = ctype( -1 );
217 template<
class ctype,
int dim >
221 conditional< (dimension > 1), MultiDimensional<true>, OneDimensional<false> > :: type
222 ::integrationOuterNormal( i, n );
225 template<
class ctype >
226 static ctype volume ()
228 typedef ReferenceDomainBase< BaseTopology > BaseReferenceDomain;
229 const ctype baseVolume = BaseReferenceDomain::template volume< ctype >();
230 return baseVolume / ctype( (
unsigned int)(dimension) );
240 template<
class Topology >
241 struct ReferenceDomain
243 static const unsigned int numCorners = Topology::numCorners;
244 static const unsigned int dimension = Topology::dimension;
246 static const unsigned int numNormals
247 = ReferenceDomainBase< Topology >::numNormals;
249 template<
class ctype >
253 ReferenceDomainBase< Topology >::corner( i, x );
256 template<
class ctype >
259 return ReferenceDomainBase< Topology >::checkInside( x, ctype( 1 ) );
262 template<
class ctype >
267 return ReferenceDomainBase< Topology >::integrationOuterNormal( i, n );
270 template<
class ctype >
271 static ctype volume ()
273 return ReferenceDomainBase< Topology >::template volume< ctype >();
282 template<
class ct,
int cdim >
284 checkInside (
unsigned int topologyId,
int dim,
const FieldVector< ct, cdim > &x, ct tolerance, ct factor = ct( 1 ) )
286 assert( (dim >= 0) && (dim <= cdim) );
287 assert( topologyId < numTopologies( dim ) );
291 const ct baseFactor = (isPrism( topologyId, dim ) ? factor : factor - x[ dim-1 ]);
292 if( (x[ dim-1 ] > -tolerance) && (factor - x[ dim-1 ] > -tolerance) )
293 return checkInside< ct, cdim >( baseTopologyId( topologyId, dim ), dim-1, x, tolerance, baseFactor );
306 template<
class ct,
int cdim >
308 referenceCorners (
unsigned int topologyId,
int dim, FieldVector< ct, cdim > *corners )
310 assert( (dim >= 0) && (dim <= cdim) );
311 assert( topologyId < numTopologies( dim ) );
315 const unsigned int nBaseCorners
316 = referenceCorners( baseTopologyId( topologyId, dim ), dim-1, corners );
317 assert( nBaseCorners == size( baseTopologyId( topologyId, dim ), dim-1, dim-1 ) );
318 if( isPrism( topologyId, dim ) )
320 std::copy( corners, corners + nBaseCorners, corners + nBaseCorners );
321 for(
unsigned int i = 0; i < nBaseCorners; ++i )
322 corners[ i+nBaseCorners ][ dim-1 ] = ct( 1 );
323 return 2*nBaseCorners;
327 corners[ nBaseCorners ] = FieldVector< ct, cdim >( ct( 0 ) );
328 corners[ nBaseCorners ][ dim-1 ] = ct( 1 );
329 return nBaseCorners+1;
334 *corners = FieldVector< ct, cdim >( ct( 0 ) );
344 unsigned long referenceVolumeInverse (
unsigned int topologyId,
int dim );
347 inline ct referenceVolume (
unsigned int topologyId,
int dim )
349 return ct( 1 ) / ct( referenceVolumeInverse( topologyId, dim ) );
357 template<
class ct,
int cdim >
359 referenceOrigins (
unsigned int topologyId,
int dim,
int codim, FieldVector< ct, cdim > *origins )
361 assert( (dim >= 0) && (dim <= cdim) );
362 assert( topologyId < numTopologies( dim ) );
363 assert( (codim >= 0) && (codim <= dim) );
367 const unsigned int baseId = baseTopologyId( topologyId, dim );
368 if( isPrism( topologyId, dim ) )
370 const unsigned int n = (codim < dim ? referenceOrigins( baseId, dim-1, codim, origins ) : 0);
371 const unsigned int m = referenceOrigins( baseId, dim-1, codim-1, origins+n );
372 for(
unsigned int i = 0; i < m; ++i )
374 origins[ n+m+i ] = origins[ n+i ];
375 origins[ n+m+i ][ dim-1 ] = ct( 1 );
381 const unsigned int m = referenceOrigins( baseId, dim-1, codim-1, origins );
384 origins[ m ] = FieldVector< ct, cdim >( ct( 0 ) );
385 origins[ m ][ dim-1 ] = ct( 1 );
389 return m+referenceOrigins( baseId, dim-1, codim, origins+m );
394 origins[ 0 ] = FieldVector< ct, cdim >( ct( 0 ) );
404 template<
class ct,
int cdim,
int mydim >
406 referenceEmbeddings (
unsigned int topologyId,
int dim,
int codim,
407 FieldVector< ct, cdim > *origins,
408 FieldMatrix< ct, mydim, cdim > *jacobianTransposeds )
410 assert( (0 <= codim) && (codim <= dim) && (dim <= cdim) );
411 assert( (dim - codim <= mydim) && (mydim <= cdim) );
412 assert( topologyId < numTopologies( dim ) );
416 const unsigned int baseId = baseTopologyId( topologyId, dim );
417 if( isPrism( topologyId, dim ) )
419 const unsigned int n = (codim < dim ? referenceEmbeddings( baseId, dim-1, codim, origins, jacobianTransposeds ) : 0);
420 for(
unsigned int i = 0; i < n; ++i )
421 jacobianTransposeds[ i ][ dim-codim-1 ][ dim-1 ] = ct( 1 );
423 const unsigned int m = referenceEmbeddings( baseId, dim-1, codim-1, origins+n, jacobianTransposeds+n );
424 std::copy( origins+n, origins+n+m, origins+n+m );
425 std::copy( jacobianTransposeds+n, jacobianTransposeds+n+m, jacobianTransposeds+n+m );
426 for(
unsigned int i = 0; i < m; ++i )
427 origins[ n+m+i ][ dim-1 ] = ct( 1 );
433 const unsigned int m = referenceEmbeddings( baseId, dim-1, codim-1, origins, jacobianTransposeds );
436 origins[ m ] = FieldVector< ct, cdim >( ct( 0 ) );
437 origins[ m ][ dim-1 ] = ct( 1 );
438 jacobianTransposeds[ m ] = FieldMatrix< ct, mydim, cdim >( ct( 0 ) );
443 const unsigned int n = referenceEmbeddings( baseId, dim-1, codim, origins+m, jacobianTransposeds+m );
444 for(
unsigned int i = 0; i < n; ++i )
446 for(
int k = 0; k < dim-1; ++k )
447 jacobianTransposeds[ m+i ][ dim-codim-1 ][ k ] = -origins[ m+i ][ k ];
448 jacobianTransposeds[ m+i ][ dim-codim-1 ][ dim-1 ] = ct( 1 );
456 origins[ 0 ] = FieldVector< ct, cdim >( ct( 0 ) );
457 jacobianTransposeds[ 0 ] = FieldMatrix< ct, mydim, cdim >( ct( 0 ) );
458 for(
int k = 0; k < dim; ++k )
459 jacobianTransposeds[ 0 ][ k ][ k ] = ct( 1 );
469 template<
class ct,
int cdim >
471 referenceIntegrationOuterNormals (
unsigned int topologyId,
int dim,
472 const FieldVector< ct, cdim > *origins,
473 FieldVector< ct, cdim > *normals )
475 assert( (dim > 0) && (dim <= cdim) );
476 assert( topologyId < numTopologies( dim ) );
480 const unsigned int baseId = baseTopologyId( topologyId, dim );
481 if( isPrism( topologyId, dim ) )
483 const unsigned int numBaseFaces
484 = referenceIntegrationOuterNormals( baseId, dim-1, origins, normals );
486 for(
unsigned int i = 0; i < 2; ++i )
488 normals[ numBaseFaces+i ] = FieldVector< ct, cdim >( ct( 0 ) );
489 normals[ numBaseFaces+i ][ dim-1 ] = ct( 2*
int( i )-1 );
492 return numBaseFaces+2;
496 normals[ 0 ] = FieldVector< ct, cdim >( ct( 0 ) );
497 normals[ 0 ][ dim-1 ] = ct( -1 );
499 const unsigned int numBaseFaces
500 = referenceIntegrationOuterNormals( baseId, dim-1, origins+1, normals+1 );
501 for(
unsigned int i = 1; i <= numBaseFaces; ++i )
502 normals[ i ][ dim-1 ] = normals[ i ]*origins[ i ];
504 return numBaseFaces+1;
509 for(
unsigned int i = 0; i < 2; ++i )
511 normals[ i ] = FieldVector< ct, cdim >( ct( 0 ) );
512 normals[ i ][ 0 ] = ct( 2*
int( i )-1 );
519 template<
class ct,
int cdim >
521 referenceIntegrationOuterNormals (
unsigned int topologyId,
int dim,
522 FieldVector< ct, cdim > *normals )
524 assert( (dim > 0) && (dim <= cdim) );
526 FieldVector< ct, cdim > *origins
527 =
new FieldVector< ct, cdim >[ size( topologyId, dim, 1 ) ];
528 referenceOrigins( topologyId, dim, 1, origins );
530 const unsigned int numFaces
531 = referenceIntegrationOuterNormals( topologyId, dim, origins, normals );
532 assert( numFaces == size( topologyId, dim, 1 ) );
Fallback implementation of the std::array class (a static array)
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.
Dune namespace.
Definition: alignment.hh:10
Traits for type conversions and type information.
Definition of the DUNE_UNUSED macro for the case that config.h is not available.
#define DUNE_UNUSED_PARAMETER(parm)
A macro to mark intentional unused function parameters with.
Definition: unused.hh:18