dune-geometry  2.4
referencedomain.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=4 sw=2 sts=2:
3 #ifndef DUNE_GEOMETRY_GENERICGEOMETRY_REFERENCEDOMAIN_HH
4 #define DUNE_GEOMETRY_GENERICGEOMETRY_REFERENCEDOMAIN_HH
5 
6 #include <algorithm>
7 
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>
13 
16 
17 namespace Dune
18 {
19 
20  namespace GenericGeometry
21  {
22 
23  // Internal Forward Declarations
24  // -----------------------------
25 
26  template< class Topology >
28 
29 
30 
31  // ReferenceDomain
32  // ---------------
33 
34  template< class Topology >
36 
38  template<>
39  class ReferenceDomainBase< Point >
40  {
41  typedef Point Topology;
42 
43  friend struct ReferenceDomain< Topology >;
44  friend class ReferenceDomainBase< Prism< Topology > >;
45  friend class ReferenceDomainBase< Pyramid< Topology > >;
46 
47  static const unsigned int numNormals = 0;
48 
49  template< class ctype, int dim >
50  static void corner ( unsigned int i, FieldVector< ctype, dim > &n )
51  {
52  DUNE_UNUSED_PARAMETER(i);
53  DUNE_UNUSED_PARAMETER(n);
54  assert( i < Topology::numCorners );
55  }
56 
57  template< class ctype, int dim >
58  static bool
59  checkInside ( const FieldVector< ctype, dim > &x, ctype factor )
60  {
61  DUNE_UNUSED_PARAMETER(x);
62  DUNE_UNUSED_PARAMETER(factor);
63  return true;
64  }
65 
66  template< class ctype, int dim >
67  static void
68  integrationOuterNormal ( unsigned int i, FieldVector< ctype, dim > &n )
69  {
70  DUNE_UNUSED_PARAMETER(i);
71  DUNE_UNUSED_PARAMETER(n);
72  assert( i < numNormals );
73  }
74 
75  template< class ctype >
76  static ctype volume ()
77  {
78  return ctype( 1 );
79  }
80  };
81 
82 
83  template< class BaseTopology >
84  class ReferenceDomainBase< Prism< BaseTopology > >
85  {
86  typedef Prism< BaseTopology > Topology;
87 
88  friend struct ReferenceDomain< Topology >;
89  friend class ReferenceDomainBase< Prism< Topology > >;
90  friend class ReferenceDomainBase< Pyramid< Topology > >;
91 
92  static const unsigned int numNormals = Size< Topology, 1 >::value;
93 
94  static const unsigned int dimension = Topology::dimension;
95  static const unsigned int myindex = dimension - 1;
96 
97  template< class ctype, int dim >
98  static void corner ( unsigned int i, FieldVector< ctype, dim > &x )
99  {
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 );
105  }
106 
107  template< class ctype, int dim >
108  static bool
109  checkInside ( const FieldVector< ctype, dim > &x, ctype factor )
110  {
111  const ctype xn = x[ myindex ];
112  const ctype cxn = factor - xn;
113  return (xn > -1e-12) && (cxn > -1e-12)
115  }
116 
117  template< class ctype, int dim >
118  static void
119  integrationOuterNormal ( unsigned int i, FieldVector< ctype, dim > &n )
120  {
121  typedef ReferenceDomainBase< BaseTopology > BaseReferenceDomain;
122 
123  if( i >= BaseReferenceDomain::numNormals )
124  {
125  const unsigned int j = i - BaseReferenceDomain::numNormals;
126  n[ myindex ] = (j == 0 ? ctype( -1 ) : ctype( 1 ));
127  }
128  else
129  BaseReferenceDomain::integrationOuterNormal( i, n );
130  }
131 
132  template< class ctype >
133  static ctype volume ()
134  {
135  typedef ReferenceDomainBase< BaseTopology > BaseReferenceDomain;
136  return BaseReferenceDomain::template volume< ctype >();
137  }
138  };
139 
140 
141  template< class BaseTopology >
142  class ReferenceDomainBase< Pyramid< BaseTopology > >
143  {
144  typedef Pyramid< BaseTopology > Topology;
145 
146  friend struct ReferenceDomain< Topology >;
147  friend class ReferenceDomainBase< Prism< Topology > >;
148  friend class ReferenceDomainBase< Pyramid< Topology > >;
149 
150  static const unsigned int numNormals = Size< Topology, 1 >::value;
151 
152  static const unsigned int dimension = Topology::dimension;
153  static const unsigned int myindex = dimension - 1;
154 
155  template< bool >
156  struct MultiDimensional
157  {
158  template< class ctype, int dim >
159  static void
160  integrationOuterNormal ( unsigned int i, FieldVector< ctype, dim > &n )
161  {
162  multiDimensionalIntegrationOuterNormal( i, n );
163  }
164  };
165 
166  template< bool >
167  struct OneDimensional
168  {
169  template< class ctype, int dim >
170  static void
171  integrationOuterNormal ( unsigned int i, FieldVector< ctype, dim > &n )
172  {
173  n[ myindex ] = (i > 0) ? ctype( 1 ) : ctype( -1 );
174  }
175  };
176 
177  template< class ctype, int dim >
178  static void corner ( unsigned int i, FieldVector< ctype, dim > &x )
179  {
180  assert( i < Topology::numCorners );
181  if( i < BaseTopology::numCorners )
182  ReferenceDomainBase< BaseTopology >::corner( i, x );
183  else
184  x[ myindex ] = ctype( 1 );
185  }
186 
187  template< class ctype, int dim >
188  static bool
189  checkInside ( const FieldVector< ctype, dim > &x, ctype factor )
190  {
191  const ctype xn = x[ myindex ];
192  const ctype cxn = factor - xn;
193  return (xn > -1e-12) && (cxn > -1e-12)
195  }
196 
197  template< class ctype, int dim >
198  static void
199  multiDimensionalIntegrationOuterNormal ( unsigned int i, FieldVector< ctype, dim > &n )
200  {
201  typedef ReferenceDomainBase< BaseTopology > BaseReferenceDomain;
202  typedef SubTopologyNumbering< BaseTopology, 1, dimension-2 > Numbering;
203 
204  if( i > 0 )
205  {
206  const unsigned int j = Numbering::number( i-1, 0 );
207  FieldVector< ctype, dim > x( ctype( 0 ) );
208  BaseReferenceDomain::corner( j, x );
209 
210  BaseReferenceDomain::integrationOuterNormal ( i-1, n );
211  n[ myindex ] = (x * n);
212  }
213  else
214  n[ myindex ] = ctype( -1 );
215  }
216 
217  template< class ctype, int dim >
218  static void
219  integrationOuterNormal ( unsigned int i, FieldVector< ctype, dim > &n )
220  {
221  conditional< (dimension > 1), MultiDimensional<true>, OneDimensional<false> > :: type
222  ::integrationOuterNormal( i, n );
223  }
224 
225  template< class ctype >
226  static ctype volume ()
227  {
228  typedef ReferenceDomainBase< BaseTopology > BaseReferenceDomain;
229  const ctype baseVolume = BaseReferenceDomain::template volume< ctype >();
230  return baseVolume / ctype( (unsigned int)(dimension) ); // linker problem when using dimension directly
231  }
232  };
237  // ReferenceDomain
238  // ---------------
239 
240  template< class Topology >
241  struct ReferenceDomain
242  {
243  static const unsigned int numCorners = Topology::numCorners;
244  static const unsigned int dimension = Topology::dimension;
245 
246  static const unsigned int numNormals
248 
249  template< class ctype >
250  static void corner ( unsigned int i, FieldVector< ctype, dimension > &x )
251  {
252  x = ctype( 0 );
254  }
255 
256  template< class ctype >
257  static bool checkInside ( const FieldVector< ctype, dimension > &x )
258  {
259  return ReferenceDomainBase< Topology >::checkInside( x, ctype( 1 ) );
260  }
261 
262  template< class ctype >
263  static void
264  integrationOuterNormal ( unsigned int i, FieldVector< ctype, dimension > &n )
265  {
266  n = ctype( 0 );
268  }
269 
270  template< class ctype >
271  static ctype volume ()
272  {
273  return ReferenceDomainBase< Topology >::template volume< ctype >();
274  }
275  };
276 
277 
278 
279  // checkInside
280  // -----------
281 
282  template< class ct, int cdim >
283  inline bool
284  checkInside ( unsigned int topologyId, int dim, const FieldVector< ct, cdim > &x, ct tolerance, ct factor = ct( 1 ) )
285  {
286  assert( (dim >= 0) && (dim <= cdim) );
287  assert( topologyId < numTopologies( dim ) );
288 
289  if( dim > 0 )
290  {
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 );
294  else
295  return false;
296  }
297  else
298  return true;
299  }
300 
301 
302 
303  // referenceCorners
304  // ----------------
305 
306  template< class ct, int cdim >
307  inline unsigned int
308  referenceCorners ( unsigned int topologyId, int dim, FieldVector< ct, cdim > *corners )
309  {
310  assert( (dim >= 0) && (dim <= cdim) );
311  assert( topologyId < numTopologies( dim ) );
312 
313  if( dim > 0 )
314  {
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 ) )
319  {
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;
324  }
325  else
326  {
327  corners[ nBaseCorners ] = FieldVector< ct, cdim >( ct( 0 ) );
328  corners[ nBaseCorners ][ dim-1 ] = ct( 1 );
329  return nBaseCorners+1;
330  }
331  }
332  else
333  {
334  *corners = FieldVector< ct, cdim >( ct( 0 ) );
335  return 1;
336  }
337  }
338 
339 
340 
341  // referenceVolume
342  // ---------------
343 
344  unsigned long referenceVolumeInverse ( unsigned int topologyId, int dim );
345 
346  template< class ct >
347  inline ct referenceVolume ( unsigned int topologyId, int dim )
348  {
349  return ct( 1 ) / ct( referenceVolumeInverse( topologyId, dim ) );
350  }
351 
352 
353 
354  // referenceOrigins
355  // ----------------
356 
357  template< class ct, int cdim >
358  inline unsigned int
359  referenceOrigins ( unsigned int topologyId, int dim, int codim, FieldVector< ct, cdim > *origins )
360  {
361  assert( (dim >= 0) && (dim <= cdim) );
362  assert( topologyId < numTopologies( dim ) );
363  assert( (codim >= 0) && (codim <= dim) );
364 
365  if( codim > 0 )
366  {
367  const unsigned int baseId = baseTopologyId( topologyId, dim );
368  if( isPrism( topologyId, dim ) )
369  {
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 )
373  {
374  origins[ n+m+i ] = origins[ n+i ];
375  origins[ n+m+i ][ dim-1 ] = ct( 1 );
376  }
377  return n+2*m;
378  }
379  else
380  {
381  const unsigned int m = referenceOrigins( baseId, dim-1, codim-1, origins );
382  if( codim == dim )
383  {
384  origins[ m ] = FieldVector< ct, cdim >( ct( 0 ) );
385  origins[ m ][ dim-1 ] = ct( 1 );
386  return m+1;
387  }
388  else
389  return m+referenceOrigins( baseId, dim-1, codim, origins+m );
390  }
391  }
392  else
393  {
394  origins[ 0 ] = FieldVector< ct, cdim >( ct( 0 ) );
395  return 1;
396  }
397  }
398 
399 
400 
401  // referenceEmbeddings
402  // -------------------
403 
404  template< class ct, int cdim, int mydim >
405  inline unsigned int
406  referenceEmbeddings ( unsigned int topologyId, int dim, int codim,
407  FieldVector< ct, cdim > *origins,
408  FieldMatrix< ct, mydim, cdim > *jacobianTransposeds )
409  {
410  assert( (0 <= codim) && (codim <= dim) && (dim <= cdim) );
411  assert( (dim - codim <= mydim) && (mydim <= cdim) );
412  assert( topologyId < numTopologies( dim ) );
413 
414  if( codim > 0 )
415  {
416  const unsigned int baseId = baseTopologyId( topologyId, dim );
417  if( isPrism( topologyId, dim ) )
418  {
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 );
422 
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 );
428 
429  return n+2*m;
430  }
431  else
432  {
433  const unsigned int m = referenceEmbeddings( baseId, dim-1, codim-1, origins, jacobianTransposeds );
434  if( codim == dim )
435  {
436  origins[ m ] = FieldVector< ct, cdim >( ct( 0 ) );
437  origins[ m ][ dim-1 ] = ct( 1 );
438  jacobianTransposeds[ m ] = FieldMatrix< ct, mydim, cdim >( ct( 0 ) );
439  return m+1;
440  }
441  else
442  {
443  const unsigned int n = referenceEmbeddings( baseId, dim-1, codim, origins+m, jacobianTransposeds+m );
444  for( unsigned int i = 0; i < n; ++i )
445  {
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 );
449  }
450  return m+n;
451  }
452  }
453  }
454  else
455  {
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 );
460  return 1;
461  }
462  }
463 
464 
465 
466  // referenceIntegrationOuterNormals
467  // --------------------------------
468 
469  template< class ct, int cdim >
470  inline unsigned int
471  referenceIntegrationOuterNormals ( unsigned int topologyId, int dim,
472  const FieldVector< ct, cdim > *origins,
473  FieldVector< ct, cdim > *normals )
474  {
475  assert( (dim > 0) && (dim <= cdim) );
476  assert( topologyId < numTopologies( dim ) );
477 
478  if( dim > 1 )
479  {
480  const unsigned int baseId = baseTopologyId( topologyId, dim );
481  if( isPrism( topologyId, dim ) )
482  {
483  const unsigned int numBaseFaces
484  = referenceIntegrationOuterNormals( baseId, dim-1, origins, normals );
485 
486  for( unsigned int i = 0; i < 2; ++i )
487  {
488  normals[ numBaseFaces+i ] = FieldVector< ct, cdim >( ct( 0 ) );
489  normals[ numBaseFaces+i ][ dim-1 ] = ct( 2*int( i )-1 );
490  }
491 
492  return numBaseFaces+2;
493  }
494  else
495  {
496  normals[ 0 ] = FieldVector< ct, cdim >( ct( 0 ) );
497  normals[ 0 ][ dim-1 ] = ct( -1 );
498 
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 ];
503 
504  return numBaseFaces+1;
505  }
506  }
507  else
508  {
509  for( unsigned int i = 0; i < 2; ++i )
510  {
511  normals[ i ] = FieldVector< ct, cdim >( ct( 0 ) );
512  normals[ i ][ 0 ] = ct( 2*int( i )-1 );
513  }
514 
515  return 2;
516  }
517  }
518 
519  template< class ct, int cdim >
520  inline unsigned int
521  referenceIntegrationOuterNormals ( unsigned int topologyId, int dim,
522  FieldVector< ct, cdim > *normals )
523  {
524  assert( (dim > 0) && (dim <= cdim) );
525 
526  FieldVector< ct, cdim > *origins
527  = new FieldVector< ct, cdim >[ size( topologyId, dim, 1 ) ];
528  referenceOrigins( topologyId, dim, 1, origins );
529 
530  const unsigned int numFaces
531  = referenceIntegrationOuterNormals( topologyId, dim, origins, normals );
532  assert( numFaces == size( topologyId, dim, 1 ) );
533 
534  delete[] origins;
535 
536  return numFaces;
537  }
538 
539  } // namespace GenericGeometry
540 
541 } // namespace Dune
542 
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