3 #ifndef DUNE_GEOMETRY_GENERICGEOMETRY_REFERENCEDOMAIN_HH
4 #define DUNE_GEOMETRY_GENERICGEOMETRY_REFERENCEDOMAIN_HH
8 #include <dune/common/array.hh>
9 #include <dune/common/fmatrix.hh>
10 #include <dune/common/fvector.hh>
11 #include <dune/common/typetraits.hh>
12 #include <dune/common/unused.hh>
20 namespace GenericGeometry
26 template<
class Topology >
34 template<
class Topology >
47 static const unsigned int numNormals = 0;
49 template<
class ctype,
int dim >
50 static void corner (
unsigned int i, FieldVector< ctype, dim > &n )
52 DUNE_UNUSED_PARAMETER(i);
53 DUNE_UNUSED_PARAMETER(n);
54 assert( i < Topology::numCorners );
57 template<
class ctype,
int dim >
59 checkInside (
const FieldVector< ctype, dim > &x, ctype factor )
61 DUNE_UNUSED_PARAMETER(x);
62 DUNE_UNUSED_PARAMETER(factor);
66 template<
class ctype,
int dim >
68 integrationOuterNormal (
unsigned int i, FieldVector< ctype, dim > &n )
70 DUNE_UNUSED_PARAMETER(i);
71 DUNE_UNUSED_PARAMETER(n);
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 >
98 static void corner (
unsigned int i, FieldVector< ctype, dim > &x )
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 >
109 checkInside (
const FieldVector< ctype, dim > &x, ctype factor )
111 const ctype xn = x[ myindex ];
112 const ctype cxn = factor - xn;
113 return (xn > -1e-12) && (cxn > -1e-12)
117 template<
class ctype,
int dim >
119 integrationOuterNormal (
unsigned int i, FieldVector< ctype, dim > &n )
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 >
160 integrationOuterNormal (
unsigned int i, FieldVector< ctype, dim > &n )
162 multiDimensionalIntegrationOuterNormal( i, n );
167 struct OneDimensional
169 template<
class ctype,
int dim >
171 integrationOuterNormal (
unsigned int i, FieldVector< ctype, dim > &n )
173 n[ myindex ] = (i > 0) ? ctype( 1 ) : ctype( -1 );
177 template<
class ctype,
int dim >
178 static void corner (
unsigned int i, FieldVector< ctype, dim > &x )
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 >
189 checkInside (
const FieldVector< ctype, dim > &x, ctype factor )
191 const ctype xn = x[ myindex ];
192 const ctype cxn = factor - xn;
193 return (xn > -1e-12) && (cxn > -1e-12)
197 template<
class ctype,
int dim >
199 multiDimensionalIntegrationOuterNormal (
unsigned int i, FieldVector< ctype, dim > &n )
201 typedef ReferenceDomainBase< BaseTopology > BaseReferenceDomain;
202 typedef SubTopologyNumbering< BaseTopology, 1, dimension-2 > Numbering;
206 const unsigned int j = Numbering::number( i-1, 0 );
207 FieldVector< ctype, dim > x( ctype( 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 >
219 integrationOuterNormal (
unsigned int i, FieldVector< ctype, dim > &n )
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
244 static const unsigned int dimension = Topology::dimension;
246 static const unsigned int numNormals
249 template<
class ctype >
250 static void corner (
unsigned int i, FieldVector< ctype, dimension > &x )
256 template<
class ctype >
257 static bool checkInside (
const FieldVector< ctype, dimension > &x )
262 template<
class ctype >
270 template<
class 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) );
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 >
310 assert( (dim >= 0) && (dim <= cdim) );
315 const unsigned int nBaseCorners
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 ) );
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) );
363 assert( (codim >= 0) && (codim <= 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 );
384 origins[ m ] = FieldVector< ct, cdim >( ct( 0 ) );
385 origins[ m ][ dim-1 ] = ct( 1 );
394 origins[ 0 ] = FieldVector< ct, cdim >( ct( 0 ) );
404 template<
class ct,
int cdim,
int mydim >
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) );
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 >
472 const FieldVector< ct, cdim > *origins,
473 FieldVector< ct, cdim > *normals )
475 assert( (dim > 0) && (dim <= cdim) );
481 if(
isPrism( topologyId, dim ) )
483 const unsigned int numBaseFaces
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
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 >
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 ) ];
530 const unsigned int numFaces
532 assert( numFaces ==
size( topologyId, dim, 1 ) );
543 #endif // DUNE_GEOMETRY_GENERICGEOMETRY_REFERENCEDOMAIN_HH
unsigned int referenceOrigins(unsigned int topologyId, int dim, int codim, FieldVector< ct, cdim > *origins)
Definition: referencedomain.hh:359
static bool checkInside(const FieldVector< ctype, dimension > &x)
Definition: referencedomain.hh:257
Definition: referencedomain.hh:27
Definition: topologytypes.hh:40
Definition: affinegeometry.hh:18
unsigned int referenceCorners(unsigned int topologyId, int dim, FieldVector< ct, cdim > *corners)
Definition: referencedomain.hh:308
unsigned int size(unsigned int topologyId, int dim, int codim)
Compute the number of subentities of a given codimension.
Definition: subtopologies.cc:16
unsigned int baseTopologyId(unsigned int topologyId, int dim, int codim=1)
obtain the base topology of a given codimension
Definition: topologytypes.hh:201
Definition: topologytypes.hh:56
Definition: topologytypes.hh:25
unsigned long referenceVolumeInverse(unsigned int topologyId, int dim)
Definition: referencedomain.cc:16
static ctype volume()
Definition: referencedomain.hh:271
unsigned int referenceEmbeddings(unsigned int topologyId, int dim, int codim, FieldVector< ct, cdim > *origins, FieldMatrix< ct, mydim, cdim > *jacobianTransposeds)
Definition: referencedomain.hh:406
bool checkInside(unsigned int topologyId, int dim, const FieldVector< ct, cdim > &x, ct tolerance, ct factor=ct(1))
Definition: referencedomain.hh:284
bool isPrism(unsigned int topologyId, int dim, int codim=0)
check whether a prism construction was used to create a given codimension
Definition: topologytypes.hh:168
static const unsigned int numCorners
Definition: referencedomain.hh:243
static void integrationOuterNormal(unsigned int i, FieldVector< ctype, dimension > &n)
Definition: referencedomain.hh:264
Definition: topologytypes.hh:270
Definition: referencedomain.hh:35
unsigned int numTopologies(int dim)
obtain the number of topologies of a given dimension
Definition: topologytypes.hh:134
ct referenceVolume(unsigned int topologyId, int dim)
Definition: referencedomain.hh:347
unsigned int referenceIntegrationOuterNormals(unsigned int topologyId, int dim, const FieldVector< ct, cdim > *origins, FieldVector< ct, cdim > *normals)
Definition: referencedomain.hh:471
static void corner(unsigned int i, FieldVector< ctype, dimension > &x)
Definition: referencedomain.hh:250