dune-geometry  3.0-git
referenceelements.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_REFERENCEELEMENTS_HH
4 #define DUNE_GEOMETRY_REFERENCEELEMENTS_HH
5 
6 #include <algorithm>
7 #include <limits>
8 
9 #include <dune/common/forloop.hh>
10 #include <dune/common/typetraits.hh>
11 #include <dune/common/visibility.hh>
12 #include <dune/common/unused.hh>
13 
18 
19 namespace Dune
20 {
21 
22  // Internal Forward Declarations
23  // -----------------------------
24 
25  template< class ctype, int dim >
26  class ReferenceElementContainer;
27 
28  template< class ctype, int dim >
29  struct ReferenceElements;
30 
31 
32 
51  template< class ctype, int dim >
52  class ReferenceElement
53  {
55 
56  friend class ReferenceElementContainer< ctype, dim >;
57 
58  struct SubEntityInfo;
59 
60  // make copy constructor private
61  ReferenceElement ( const This & );
62 
63  ReferenceElement () {}
64 
65  template< int codim > struct CreateGeometries;
66 
67  public:
69  template< int codim >
70  struct Codim
71  {
73  typedef AffineGeometry< ctype, dim-codim, dim > Geometry;
74  };
75 
80  int size ( int c ) const
81  {
82  assert( (c >= 0) && (c <= dim) );
83  return info_[ c ].size();
84  }
85 
97  int size ( int i, int c, int cc ) const
98  {
99  assert( (i >= 0) && (i < size( c )) );
100  return info_[ c ][ i ].size( cc );
101  }
102 
116  int subEntity ( int i, int c, int ii, int cc ) const
117  {
118  assert( (i >= 0) && (i < size( c )) );
119  return info_[ c ][ i ].number( ii, cc );
120  }
121 
130  const GeometryType &type ( int i, int c ) const
131  {
132  assert( (i >= 0) && (i < size( c )) );
133  return info_[ c ][ i ].type();
134  }
135 
137  const GeometryType &type () const { return type( 0, 0 ); }
138 
148  const FieldVector< ctype, dim > &position( int i, int c ) const
149  {
150  assert( (c >= 0) && (c <= dim) );
151  return baryCenters_[ c ][ i ];
152  }
153 
161  bool checkInside ( const FieldVector< ctype, dim > &local ) const
162  {
163  const ctype tolerance = ctype( 64 ) * std::numeric_limits< ctype >::epsilon();
164  return GenericGeometry::template checkInside< ctype, dim >( type().id(), dim, local, tolerance );
165  }
166 
178  template< int codim >
179  typename Codim< codim >::Geometry geometry ( int i ) const
180  {
181  std::integral_constant< int, codim > codimVariable;
182  return geometries_[ codimVariable ][ i ];
183  }
184 
186  ctype volume () const
187  {
188  return volume_;
189  }
190 
198  const FieldVector< ctype, dim > &integrationOuterNormal ( int face ) const
199  {
200  assert( (face >= 0) && (face < int( integrationNormals_.size() )) );
201  return integrationNormals_[ face ];
202  }
203 
204  private:
205  void initialize ( unsigned int topologyId )
206  {
207  assert( topologyId < GenericGeometry::numTopologies( dim ) );
208 
209  // set up subentities
210  for( int codim = 0; codim <= dim; ++codim )
211  {
212  const unsigned int size = GenericGeometry::size( topologyId, dim, codim );
213  info_[ codim ].resize( size );
214  for( unsigned int i = 0; i < size; ++i )
215  info_[ codim ][ i ].initialize( topologyId, codim, i );
216  }
217 
218  // compute corners
219  const unsigned int numVertices = size( dim );
220  baryCenters_[ dim ].resize( numVertices );
221  GenericGeometry::referenceCorners( topologyId, dim, &(baryCenters_[ dim ][ 0 ]) );
222 
223  // compute barycenters
224  for( int codim = 0; codim < dim; ++codim )
225  {
226  baryCenters_[ codim ].resize( size(codim) );
227  for( int i = 0; i < size( codim ); ++i )
228  {
229  baryCenters_[ codim ][ i ] = FieldVector< ctype, dim >( ctype( 0 ) );
230  const unsigned int numCorners = size( i, codim, dim );
231  for( unsigned int j = 0; j < numCorners; ++j )
232  baryCenters_[ codim ][ i ] += baryCenters_[ dim ][ subEntity( i, codim, j, dim ) ];
233  baryCenters_[ codim ][ i ] *= ctype( 1 ) / ctype( numCorners );
234  }
235  }
236 
237  // compute reference element volume
238  volume_ = GenericGeometry::template referenceVolume< ctype >( topologyId, dim );
239 
240  // compute integration outer normals
241  if( dim > 0 )
242  {
243  integrationNormals_.resize( size( 1 ) );
244  GenericGeometry::referenceIntegrationOuterNormals( topologyId, dim, &(integrationNormals_[ 0 ]) );
245  }
246 
247  // set up geometries
248  Dune::ForLoop< CreateGeometries, 0, dim >::apply( *this, geometries_ );
249  }
250 
252  template< int codim >
253  struct GeometryArray
254  : public std::vector< typename Codim< codim >::Geometry >
255  {};
256 
258  typedef GenericGeometry::CodimTable< GeometryArray, dim > GeometryTable;
259 
261  ctype volume_;
262 
263  std::vector< FieldVector< ctype, dim > > baryCenters_[ dim+1 ];
264  std::vector< FieldVector< ctype, dim > > integrationNormals_;
265 
267  GeometryTable geometries_;
268 
269  std::vector< SubEntityInfo > info_[ dim+1 ];
270  };
271 
273  template< class ctype, int dim >
274  struct ReferenceElement< ctype, dim >::SubEntityInfo
275  {
277  : numbering_( nullptr )
278  {
279  std::fill( offset_.begin(), offset_.end(), 0 );
280  }
281 
282  SubEntityInfo ( const SubEntityInfo &other )
283  : offset_( other.offset_ ),
284  type_( other.type_ )
285  {
286  numbering_ = allocate();
287  std::copy( other.numbering_, other.numbering_ + capacity(), numbering_ );
288  }
289 
290  ~SubEntityInfo () { deallocate( numbering_ ); }
291 
292  const SubEntityInfo &operator= ( const SubEntityInfo &other )
293  {
294  type_ = other.type_;
295  offset_ = other.offset_;
296 
297  deallocate( numbering_ );
298  numbering_ = allocate();
299  std::copy( other.numbering_, other.numbering_ + capacity(), numbering_ );
300 
301  return *this;
302  }
303 
304  int size ( int cc ) const
305  {
306  assert( (cc >= codim()) && (cc <= dim) );
307  return (offset_[ cc+1 ] - offset_[ cc ]);
308  }
309 
310  int number ( int ii, int cc ) const
311  {
312  assert( (ii >= 0) && (ii < size( cc )) );
313  return numbering_[ offset_[ cc ] + ii ];
314  }
315 
316  const GeometryType &type () const { return type_; }
317 
318  void initialize ( unsigned int topologyId, int codim, unsigned int i )
319  {
320  const unsigned int subId = GenericGeometry::subTopologyId( topologyId, dim, codim, i );
321  type_ = GeometryType( subId, dim-codim );
322 
323  // compute offsets
324  for( int cc = 0; cc <= codim; ++cc )
325  offset_[ cc ] = 0;
326  for( int cc = codim; cc <= dim; ++cc )
327  offset_[ cc+1 ] = offset_[ cc ] + GenericGeometry::size( subId, dim-codim, cc-codim );
328 
329  // compute subnumbering
330  deallocate( numbering_ );
331  numbering_ = allocate();
332  for( int cc = codim; cc <= dim; ++cc )
333  GenericGeometry::subTopologyNumbering( topologyId, dim, codim, i, cc-codim, numbering_+offset_[ cc ], numbering_+offset_[ cc+1 ] );
334  }
335 
336  protected:
337  int codim () const { return dim - type().dim(); }
338 
339  unsigned int *allocate () { return (capacity() != 0 ? new unsigned int[ capacity() ] : nullptr); }
340  void deallocate ( unsigned int *ptr ) { delete[] ptr; }
341  unsigned int capacity () const { return offset_[ dim+1 ]; }
342 
343  private:
344  unsigned int *numbering_;
345  std::array< unsigned int, dim+2 > offset_;
346  GeometryType type_;
347  };
348 
349 
350  template< class ctype, int dim >
351  template< int codim >
352  struct ReferenceElement< ctype, dim >::CreateGeometries
353  {
354  template< int cc >
355  static const ReferenceElement< ctype, dim-cc > &
356  subRefElement( const ReferenceElement< ctype, dim > &refElement, int i, std::integral_constant< int, cc > )
357  {
358  return ReferenceElements< ctype, dim-cc >::general( refElement.type( i, cc ) );
359  }
360 
361  static const ReferenceElement< ctype, dim > &
362  subRefElement( const ReferenceElement< ctype, dim > &refElement, int i, std::integral_constant< int, 0 > )
363  {
364  DUNE_UNUSED_PARAMETER(i);
365  return refElement;
366  }
367 
368  static void apply ( const ReferenceElement< ctype, dim > &refElement, GeometryTable &geometries )
369  {
370  const int size = refElement.size( codim );
371  std::vector< FieldVector< ctype, dim > > origins( size );
372  std::vector< FieldMatrix< ctype, dim - codim, dim > > jacobianTransposeds( size );
373  GenericGeometry::referenceEmbeddings( refElement.type().id(), dim, codim, &(origins[ 0 ]), &(jacobianTransposeds[ 0 ]) );
374 
375  std::integral_constant< int, codim > codimVariable;
376  geometries[ codimVariable ].reserve( size );
377  for( int i = 0; i < size; ++i )
378  {
379  typename Codim< codim >::Geometry geometry( subRefElement( refElement, i, codimVariable ), origins[ i ], jacobianTransposeds[ i ] );
380  geometries[ codimVariable ].push_back( geometry );
381  }
382  }
383  };
384 
385 
386 
387  // ReferenceElementContainer
388  // -------------------------
389 
390  template< class ctype, int dim >
392  {
393  static const unsigned int numTopologies = (1u << dim);
394 
395  public:
397  typedef const value_type *const_iterator;
398 
400  {
401  for( unsigned int topologyId = 0; topologyId < numTopologies; ++topologyId )
402  values_[ topologyId ].initialize( topologyId );
403  }
404 
405  const value_type &operator() ( const GeometryType &type ) const
406  {
407  assert( type.dim() == dim );
408  return values_[ type.id() ];
409  }
410 
411  const value_type &simplex () const
412  {
414  }
415 
416  const value_type &cube () const
417  {
419  }
420 
421  const value_type &pyramid () const
422  {
424  }
425 
426  const value_type &prism () const
427  {
429  }
430 
431  const_iterator begin () const { return values_; }
432  const_iterator end () const { return values_ + numTopologies; }
433 
434  private:
435  value_type values_[ numTopologies ];
436  };
437 
438 
439 
440  // ReferenceElements
441  // ------------------------
442 
453  template< class ctype, int dim >
455  {
457 
459  static const ReferenceElement< ctype, dim > &
460  general ( const GeometryType &type )
461  {
462  return container() ( type );
463  }
464 
467  {
468  return container().simplex();
469  }
470 
473  {
474  return container().cube();
475  }
476 
477  static Iterator begin () { return container().begin(); }
478  static Iterator end () { return container().end(); }
479 
480  private:
481  DUNE_EXPORT static const ReferenceElementContainer< ctype, dim > &container ()
482  {
484  return container;
485  }
486  };
487 
488 } // namespace Dune
489 
490 #endif // #ifndef DUNE_GEOMETRY_REFERENCEELEMENTS_HH
An implementation of the Geometry interface for affine geometries.
Definition: affinegeometry.hh:19
void subTopologyNumbering(unsigned int topologyId, int dim, int codim, unsigned int i, int subcodim, unsigned int *beginOut, unsigned int *endOut)
Definition: subtopologies.cc:85
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 numTopologies(int dim)
obtain the number of topologies of a given dimension
Definition: topologytypes.hh:134
unsigned int referenceCorners(unsigned int topologyId, int dim, FieldVector< ct, cdim > *corners)
Definition: referencedomain.hh:308
unsigned int referenceIntegrationOuterNormals(unsigned int topologyId, int dim, const FieldVector< ct, cdim > *origins, FieldVector< ct, cdim > *normals)
Definition: referencedomain.hh:471
unsigned int subTopologyId(unsigned int topologyId, int dim, int codim, unsigned int i)
Compute the topology id of a given subentity.
Definition: subtopologies.cc:47
unsigned int referenceEmbeddings(unsigned int topologyId, int dim, int codim, FieldVector< ct, cdim > *origins, FieldMatrix< ct, mydim, cdim > *jacobianTransposeds)
Definition: referencedomain.hh:406
This class provides access to geometric and topological properties of a reference element.
Definition: referenceelements.hh:29
int size(int c) const
number of subentities of codimension c
Definition: referenceelements.hh:80
int size(int i, int c, int cc) const
number of subentities of codimension cc of subentity (i,c)
Definition: referenceelements.hh:97
ctype volume() const
obtain the volume of the reference element
Definition: referenceelements.hh:186
int subEntity(int i, int c, int ii, int cc) const
obtain number of ii-th subentity with codim cc of (i,c)
Definition: referenceelements.hh:116
bool checkInside(const FieldVector< ctype, dim > &local) const
check if a coordinate is in the reference element
Definition: referenceelements.hh:161
const FieldVector< ctype, dim > & position(int i, int c) const
position of the barycenter of entity (i,c)
Definition: referenceelements.hh:148
Codim< codim >::Geometry geometry(int i) const
obtain the embedding of subentity (i,codim) into the reference element
Definition: referenceelements.hh:179
const GeometryType & type(int i, int c) const
obtain the type of subentity (i,c)
Definition: referenceelements.hh:130
const GeometryType & type() const
obtain the type of this reference element
Definition: referenceelements.hh:137
const FieldVector< ctype, dim > & integrationOuterNormal(int face) const
obtain the integration outer normal of the reference element
Definition: referenceelements.hh:198
Class providing access to the singletons of the reference elements.
Definition: referenceelements.hh:455
static const ReferenceElement< ctype, dim > & simplex()
get simplex reference elements
Definition: referenceelements.hh:466
static Iterator end()
Definition: referenceelements.hh:478
static const ReferenceElement< ctype, dim > & cube()
get hypercube reference elements
Definition: referenceelements.hh:472
ReferenceElementContainer< ctype, dim >::const_iterator Iterator
Definition: referenceelements.hh:456
static const ReferenceElement< ctype, dim > & general(const GeometryType &type)
get general reference elements
Definition: referenceelements.hh:460
static Iterator begin()
Definition: referenceelements.hh:477
Implementation of the Geometry interface for affine geometries.
Definition: affinegeometry.hh:39
Definition: topologytypes.hh:215
Definition: topologytypes.hh:232
Definition: topologytypes.hh:249
Definition: topologytypes.hh:260
Definition: referenceelements.hh:392
const value_type & simplex() const
Definition: referenceelements.hh:411
const value_type & pyramid() const
Definition: referenceelements.hh:421
const value_type * const_iterator
Definition: referenceelements.hh:397
const value_type & cube() const
Definition: referenceelements.hh:416
ReferenceElementContainer()
Definition: referenceelements.hh:399
const value_type & prism() const
Definition: referenceelements.hh:426
ReferenceElement< ctype, dim > value_type
Definition: referenceelements.hh:396
const value_type & operator()(const GeometryType &type) const
Definition: referenceelements.hh:405
const_iterator end() const
Definition: referenceelements.hh:432
const_iterator begin() const
Definition: referenceelements.hh:431
Collection of types depending on the codimension.
Definition: referenceelements.hh:71
AffineGeometry< ctype, dim-codim, dim > Geometry
type of geometry embedding a subentity into the reference element
Definition: referenceelements.hh:73
topological information about the subentities of a reference element
Definition: referenceelements.hh:275
int number(int ii, int cc) const
Definition: referenceelements.hh:310
int size(int cc) const
Definition: referenceelements.hh:304
~SubEntityInfo()
Definition: referenceelements.hh:290
unsigned int * allocate()
Definition: referenceelements.hh:339
const GeometryType & type() const
Definition: referenceelements.hh:316
SubEntityInfo()
Definition: referenceelements.hh:276
SubEntityInfo(const SubEntityInfo &other)
Definition: referenceelements.hh:282
unsigned int capacity() const
Definition: referenceelements.hh:341
void deallocate(unsigned int *ptr)
Definition: referenceelements.hh:340
void initialize(unsigned int topologyId, int codim, unsigned int i)
Definition: referenceelements.hh:318
int codim() const
Definition: referenceelements.hh:337
Unique label for each type of entities that can occur in DUNE grids.
Definition: type.hh:25
unsigned int dim() const
Return dimension of the type.
Definition: type.hh:321
unsigned int id() const
Return the topology id the type.
Definition: type.hh:326