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
19namespace 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 >
53 {
55
56 friend class ReferenceElementContainer< ctype, dim >;
57
58 struct SubEntityInfo;
59
60 // make copy constructor private
61 ReferenceElement ( const This & );
62
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
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:
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:53
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
const FieldVector< ctype, dim > & position(int i, int c) const
position of the barycenter of entity (i,c)
Definition referenceelements.hh:148
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 GeometryType & type(int i, int c) const
obtain the type of subentity (i,c)
Definition referenceelements.hh:130
Codim< codim >::Geometry geometry(int i) const
obtain the embedding of subentity (i,codim) into the reference element
Definition referenceelements.hh:179
const FieldVector< ctype, dim > & integrationOuterNormal(int face) const
obtain the integration outer normal of the reference element
Definition referenceelements.hh:198
const GeometryType & type() const
obtain the type of this reference element
Definition referenceelements.hh:137
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 const ReferenceElement< ctype, dim > & cube()
get hypercube reference elements
Definition referenceelements.hh:472
static Iterator end()
Definition referenceelements.hh:478
static const ReferenceElement< ctype, dim > & general(const GeometryType &type)
get general reference elements
Definition referenceelements.hh:460
static Iterator begin()
Definition referenceelements.hh:477
ReferenceElementContainer< ctype, dim >::const_iterator Iterator
Definition referenceelements.hh:456
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 * const_iterator
Definition referenceelements.hh:397
const value_type & operator()(const GeometryType &type) const
Definition referenceelements.hh:405
ReferenceElementContainer()
Definition referenceelements.hh:399
const value_type & pyramid() const
Definition referenceelements.hh:421
ReferenceElement< ctype, dim > value_type
Definition referenceelements.hh:396
const_iterator end() const
Definition referenceelements.hh:432
const value_type & cube() const
Definition referenceelements.hh:416
const_iterator begin() const
Definition referenceelements.hh:431
const value_type & prism() const
Definition referenceelements.hh:426
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
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
const GeometryType & type() const
Definition referenceelements.hh:316
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