dune-geometry  2.6-git
referenceelementimplementation.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_REFERENCEELEMENTIMPLEMENTATION_HH
4 #define DUNE_GEOMETRY_REFERENCEELEMENTIMPLEMENTATION_HH
5 
6 #include <cassert>
7 
8 #include <algorithm>
9 #include <limits>
10 #include <tuple>
11 #include <utility>
12 #include <vector>
13 #include <array>
14 
15 #include <dune/common/fmatrix.hh>
16 #include <dune/common/fvector.hh>
17 #include <dune/common/hybridutilities.hh>
18 #include <dune/common/typetraits.hh>
19 #include <dune/common/unused.hh>
20 
23 #include <dune/geometry/type.hh>
24 
25 namespace Dune
26 {
27 
28  namespace Geo
29  {
30 
31 #ifndef DOXYGEN
32 
33  // Internal Forward Declarations
34  // -----------------------------
35 
36  namespace Impl
37  {
38  template< class ctype, int dim >
39  class ReferenceElementContainer;
40  }
41 
42  template< class ctype, int dim >
43  struct ReferenceElements;
44 
45 
46 
47  namespace Impl
48  {
49 
50  using Dune::Impl::isPrism;
51  using Dune::Impl::isPyramid;
52  using Dune::Impl::baseTopologyId;
53  using Dune::Impl::prismConstruction;
54  using Dune::Impl::pyramidConstruction;
55  using Dune::Impl::numTopologies;
56 
58  unsigned int size ( unsigned int topologyId, int dim, int codim );
59 
60 
61 
69  unsigned int subTopologyId ( unsigned int topologyId, int dim, int codim, unsigned int i );
70 
71 
72 
73  // subTopologyNumbering
74  // --------------------
75 
76  void subTopologyNumbering ( unsigned int topologyId, int dim, int codim, unsigned int i, int subcodim,
77  unsigned int *beginOut, unsigned int *endOut );
78 
79 
80 
81 
82  // checkInside
83  // -----------
84 
85  template< class ct, int cdim >
86  inline bool
87  checkInside ( unsigned int topologyId, int dim, const FieldVector< ct, cdim > &x, ct tolerance, ct factor = ct( 1 ) )
88  {
89  assert( (dim >= 0) && (dim <= cdim) );
90  assert( topologyId < numTopologies( dim ) );
91 
92  if( dim > 0 )
93  {
94  const ct baseFactor = (isPrism( topologyId, dim ) ? factor : factor - x[ dim-1 ]);
95  if( (x[ dim-1 ] > -tolerance) && (factor - x[ dim-1 ] > -tolerance) )
96  return checkInside< ct, cdim >( baseTopologyId( topologyId, dim ), dim-1, x, tolerance, baseFactor );
97  else
98  return false;
99  }
100  else
101  return true;
102  }
103 
104 
105 
106  // referenceCorners
107  // ----------------
108 
109  template< class ct, int cdim >
110  inline unsigned int
111  referenceCorners ( unsigned int topologyId, int dim, FieldVector< ct, cdim > *corners )
112  {
113  assert( (dim >= 0) && (dim <= cdim) );
114  assert( topologyId < numTopologies( dim ) );
115 
116  if( dim > 0 )
117  {
118  const unsigned int nBaseCorners
119  = referenceCorners( baseTopologyId( topologyId, dim ), dim-1, corners );
120  assert( nBaseCorners == size( baseTopologyId( topologyId, dim ), dim-1, dim-1 ) );
121  if( isPrism( topologyId, dim ) )
122  {
123  std::copy( corners, corners + nBaseCorners, corners + nBaseCorners );
124  for( unsigned int i = 0; i < nBaseCorners; ++i )
125  corners[ i+nBaseCorners ][ dim-1 ] = ct( 1 );
126  return 2*nBaseCorners;
127  }
128  else
129  {
130  corners[ nBaseCorners ] = FieldVector< ct, cdim >( ct( 0 ) );
131  corners[ nBaseCorners ][ dim-1 ] = ct( 1 );
132  return nBaseCorners+1;
133  }
134  }
135  else
136  {
137  *corners = FieldVector< ct, cdim >( ct( 0 ) );
138  return 1;
139  }
140  }
141 
142 
143 
144  // referenceVolume
145  // ---------------
146 
147  unsigned long referenceVolumeInverse ( unsigned int topologyId, int dim );
148 
149  template< class ct >
150  inline ct referenceVolume ( unsigned int topologyId, int dim )
151  {
152  return ct( 1 ) / ct( referenceVolumeInverse( topologyId, dim ) );
153  }
154 
155 
156 
157  // referenceOrigins
158  // ----------------
159 
160  template< class ct, int cdim >
161  inline unsigned int
162  referenceOrigins ( unsigned int topologyId, int dim, int codim, FieldVector< ct, cdim > *origins )
163  {
164  assert( (dim >= 0) && (dim <= cdim) );
165  assert( topologyId < numTopologies( dim ) );
166  assert( (codim >= 0) && (codim <= dim) );
167 
168  if( codim > 0 )
169  {
170  const unsigned int baseId = baseTopologyId( topologyId, dim );
171  if( isPrism( topologyId, dim ) )
172  {
173  const unsigned int n = (codim < dim ? referenceOrigins( baseId, dim-1, codim, origins ) : 0);
174  const unsigned int m = referenceOrigins( baseId, dim-1, codim-1, origins+n );
175  for( unsigned int i = 0; i < m; ++i )
176  {
177  origins[ n+m+i ] = origins[ n+i ];
178  origins[ n+m+i ][ dim-1 ] = ct( 1 );
179  }
180  return n+2*m;
181  }
182  else
183  {
184  const unsigned int m = referenceOrigins( baseId, dim-1, codim-1, origins );
185  if( codim == dim )
186  {
187  origins[ m ] = FieldVector< ct, cdim >( ct( 0 ) );
188  origins[ m ][ dim-1 ] = ct( 1 );
189  return m+1;
190  }
191  else
192  return m+referenceOrigins( baseId, dim-1, codim, origins+m );
193  }
194  }
195  else
196  {
197  origins[ 0 ] = FieldVector< ct, cdim >( ct( 0 ) );
198  return 1;
199  }
200  }
201 
202 
203 
204  // referenceEmbeddings
205  // -------------------
206 
207  template< class ct, int cdim, int mydim >
208  inline unsigned int
209  referenceEmbeddings ( unsigned int topologyId, int dim, int codim,
210  FieldVector< ct, cdim > *origins,
211  FieldMatrix< ct, mydim, cdim > *jacobianTransposeds )
212  {
213  assert( (0 <= codim) && (codim <= dim) && (dim <= cdim) );
214  assert( (dim - codim <= mydim) && (mydim <= cdim) );
215  assert( topologyId < numTopologies( dim ) );
216 
217  if( codim > 0 )
218  {
219  const unsigned int baseId = baseTopologyId( topologyId, dim );
220  if( isPrism( topologyId, dim ) )
221  {
222  const unsigned int n = (codim < dim ? referenceEmbeddings( baseId, dim-1, codim, origins, jacobianTransposeds ) : 0);
223  for( unsigned int i = 0; i < n; ++i )
224  jacobianTransposeds[ i ][ dim-codim-1 ][ dim-1 ] = ct( 1 );
225 
226  const unsigned int m = referenceEmbeddings( baseId, dim-1, codim-1, origins+n, jacobianTransposeds+n );
227  std::copy( origins+n, origins+n+m, origins+n+m );
228  std::copy( jacobianTransposeds+n, jacobianTransposeds+n+m, jacobianTransposeds+n+m );
229  for( unsigned int i = 0; i < m; ++i )
230  origins[ n+m+i ][ dim-1 ] = ct( 1 );
231 
232  return n+2*m;
233  }
234  else
235  {
236  const unsigned int m = referenceEmbeddings( baseId, dim-1, codim-1, origins, jacobianTransposeds );
237  if( codim == dim )
238  {
239  origins[ m ] = FieldVector< ct, cdim >( ct( 0 ) );
240  origins[ m ][ dim-1 ] = ct( 1 );
241  jacobianTransposeds[ m ] = FieldMatrix< ct, mydim, cdim >( ct( 0 ) );
242  return m+1;
243  }
244  else
245  {
246  const unsigned int n = referenceEmbeddings( baseId, dim-1, codim, origins+m, jacobianTransposeds+m );
247  for( unsigned int i = 0; i < n; ++i )
248  {
249  for( int k = 0; k < dim-1; ++k )
250  jacobianTransposeds[ m+i ][ dim-codim-1 ][ k ] = -origins[ m+i ][ k ];
251  jacobianTransposeds[ m+i ][ dim-codim-1 ][ dim-1 ] = ct( 1 );
252  }
253  return m+n;
254  }
255  }
256  }
257  else
258  {
259  origins[ 0 ] = FieldVector< ct, cdim >( ct( 0 ) );
260  jacobianTransposeds[ 0 ] = FieldMatrix< ct, mydim, cdim >( ct( 0 ) );
261  for( int k = 0; k < dim; ++k )
262  jacobianTransposeds[ 0 ][ k ][ k ] = ct( 1 );
263  return 1;
264  }
265  }
266 
267 
268 
269  // referenceIntegrationOuterNormals
270  // --------------------------------
271 
272  template< class ct, int cdim >
273  inline unsigned int
274  referenceIntegrationOuterNormals ( unsigned int topologyId, int dim,
275  const FieldVector< ct, cdim > *origins,
276  FieldVector< ct, cdim > *normals )
277  {
278  assert( (dim > 0) && (dim <= cdim) );
279  assert( topologyId < numTopologies( dim ) );
280 
281  if( dim > 1 )
282  {
283  const unsigned int baseId = baseTopologyId( topologyId, dim );
284  if( isPrism( topologyId, dim ) )
285  {
286  const unsigned int numBaseFaces
287  = referenceIntegrationOuterNormals( baseId, dim-1, origins, normals );
288 
289  for( unsigned int i = 0; i < 2; ++i )
290  {
291  normals[ numBaseFaces+i ] = FieldVector< ct, cdim >( ct( 0 ) );
292  normals[ numBaseFaces+i ][ dim-1 ] = ct( 2*int( i )-1 );
293  }
294 
295  return numBaseFaces+2;
296  }
297  else
298  {
299  normals[ 0 ] = FieldVector< ct, cdim >( ct( 0 ) );
300  normals[ 0 ][ dim-1 ] = ct( -1 );
301 
302  const unsigned int numBaseFaces
303  = referenceIntegrationOuterNormals( baseId, dim-1, origins+1, normals+1 );
304  for( unsigned int i = 1; i <= numBaseFaces; ++i )
305  normals[ i ][ dim-1 ] = normals[ i ]*origins[ i ];
306 
307  return numBaseFaces+1;
308  }
309  }
310  else
311  {
312  for( unsigned int i = 0; i < 2; ++i )
313  {
314  normals[ i ] = FieldVector< ct, cdim >( ct( 0 ) );
315  normals[ i ][ 0 ] = ct( 2*int( i )-1 );
316  }
317 
318  return 2;
319  }
320  }
321 
322  template< class ct, int cdim >
323  inline unsigned int
324  referenceIntegrationOuterNormals ( unsigned int topologyId, int dim,
325  FieldVector< ct, cdim > *normals )
326  {
327  assert( (dim > 0) && (dim <= cdim) );
328 
329  FieldVector< ct, cdim > *origins
330  = new FieldVector< ct, cdim >[ size( topologyId, dim, 1 ) ];
331  referenceOrigins( topologyId, dim, 1, origins );
332 
333  const unsigned int numFaces
334  = referenceIntegrationOuterNormals( topologyId, dim, origins, normals );
335  assert( numFaces == size( topologyId, dim, 1 ) );
336 
337  delete[] origins;
338 
339  return numFaces;
340  }
341 
342  } // namespace Impl
343 
344 
345 
346  // ReferenceElement
347  // ----------------
348 
367  template< class ctype_, int dim >
368  class ReferenceElementImplementation
369  {
370 
371  public:
372 
374  using ctype = ctype_;
375 
377  using CoordinateField = ctype;
378 
380  using Coordinate = Dune::FieldVector<ctype,dim>;
381 
383  static constexpr int dimension = dim;
384 
385  private:
386 
387  friend class Impl::ReferenceElementContainer< ctype, dim >;
388 
389  struct SubEntityInfo;
390 
391  template< int codim > struct CreateGeometries;
392 
393  public:
395  template< int codim >
396  struct Codim
397  {
399  typedef AffineGeometry< ctype, dim-codim, dim > Geometry;
400  };
401 
402  // ReferenceElement cannot be copied.
403  ReferenceElementImplementation ( const ReferenceElementImplementation& ) = delete;
404 
405  // ReferenceElementImplementation cannot be copied.
406  ReferenceElementImplementation& operator= ( const ReferenceElementImplementation& ) = delete;
407 
408  // ReferenceElementImplementation is default-constructible (required for storage in std::array)
409  ReferenceElementImplementation () = default;
410 
415  int size ( int c ) const
416  {
417  assert( (c >= 0) && (c <= dim) );
418  return info_[ c ].size();
419  }
420 
432  int size ( int i, int c, int cc ) const
433  {
434  assert( (i >= 0) && (i < size( c )) );
435  return info_[ c ][ i ].size( cc );
436  }
437 
451  int subEntity ( int i, int c, int ii, int cc ) const
452  {
453  assert( (i >= 0) && (i < size( c )) );
454  return info_[ c ][ i ].number( ii, cc );
455  }
456 
465  const GeometryType &type ( int i, int c ) const
466  {
467  assert( (i >= 0) && (i < size( c )) );
468  return info_[ c ][ i ].type();
469  }
470 
472  const GeometryType &type () const { return type( 0, 0 ); }
473 
483  const Coordinate &position( int i, int c ) const
484  {
485  assert( (c >= 0) && (c <= dim) );
486  return baryCenters_[ c ][ i ];
487  }
488 
496  bool checkInside ( const Coordinate &local ) const
497  {
498  const ctype tolerance = ctype( 64 ) * std::numeric_limits< ctype >::epsilon();
499  return Impl::template checkInside< ctype, dim >( type().id(), dim, local, tolerance );
500  }
501 
513  template< int codim >
514  typename Codim< codim >::Geometry geometry ( int i ) const
515  {
516  return std::get< codim >( geometries_ )[ i ];
517  }
518 
520  ctype volume () const
521  {
522  return volume_;
523  }
524 
532  const Coordinate &integrationOuterNormal ( int face ) const
533  {
534  assert( (face >= 0) && (face < int( integrationNormals_.size() )) );
535  return integrationNormals_[ face ];
536  }
537 
538  private:
539  void initialize ( unsigned int topologyId )
540  {
541  assert( topologyId < Impl::numTopologies( dim ) );
542 
543  // set up subentities
544  for( int codim = 0; codim <= dim; ++codim )
545  {
546  const unsigned int size = Impl::size( topologyId, dim, codim );
547  info_[ codim ].resize( size );
548  for( unsigned int i = 0; i < size; ++i )
549  info_[ codim ][ i ].initialize( topologyId, codim, i );
550  }
551 
552  // compute corners
553  const unsigned int numVertices = size( dim );
554  baryCenters_[ dim ].resize( numVertices );
555  Impl::referenceCorners( topologyId, dim, &(baryCenters_[ dim ][ 0 ]) );
556 
557  // compute barycenters
558  for( int codim = 0; codim < dim; ++codim )
559  {
560  baryCenters_[ codim ].resize( size(codim) );
561  for( int i = 0; i < size( codim ); ++i )
562  {
563  baryCenters_[ codim ][ i ] = Coordinate( ctype( 0 ) );
564  const unsigned int numCorners = size( i, codim, dim );
565  for( unsigned int j = 0; j < numCorners; ++j )
566  baryCenters_[ codim ][ i ] += baryCenters_[ dim ][ subEntity( i, codim, j, dim ) ];
567  baryCenters_[ codim ][ i ] *= ctype( 1 ) / ctype( numCorners );
568  }
569  }
570 
571  // compute reference element volume
572  volume_ = Impl::template referenceVolume< ctype >( topologyId, dim );
573 
574  // compute integration outer normals
575  if( dim > 0 )
576  {
577  integrationNormals_.resize( size( 1 ) );
578  Impl::referenceIntegrationOuterNormals( topologyId, dim, &(integrationNormals_[ 0 ]) );
579  }
580 
581  // set up geometries
582  Hybrid::forEach( std::make_index_sequence< dim+1 >{}, [ & ]( auto i ){ CreateGeometries< i >::apply( *this, geometries_ ); } );
583  }
584 
585  template< int... codim >
586  static std::tuple< std::vector< typename Codim< codim >::Geometry >... >
587  makeGeometryTable ( std::integer_sequence< int, codim... > );
588 
590  typedef decltype( makeGeometryTable( std::make_integer_sequence< int, dim+1 >() ) ) GeometryTable;
591 
593  ctype volume_;
594 
595  std::vector< Coordinate > baryCenters_[ dim+1 ];
596  std::vector< Coordinate > integrationNormals_;
597 
599  GeometryTable geometries_;
600 
601  std::vector< SubEntityInfo > info_[ dim+1 ];
602  };
603 
605  template< class ctype, int dim >
606  struct ReferenceElementImplementation< ctype, dim >::SubEntityInfo
607  {
608  SubEntityInfo ()
609  : numbering_( nullptr )
610  {
611  std::fill( offset_.begin(), offset_.end(), 0 );
612  }
613 
614  SubEntityInfo ( const SubEntityInfo &other )
615  : offset_( other.offset_ ),
616  type_( other.type_ )
617  {
618  numbering_ = allocate();
619  std::copy( other.numbering_, other.numbering_ + capacity(), numbering_ );
620  }
621 
622  ~SubEntityInfo () { deallocate( numbering_ ); }
623 
624  const SubEntityInfo &operator= ( const SubEntityInfo &other )
625  {
626  type_ = other.type_;
627  offset_ = other.offset_;
628 
629  deallocate( numbering_ );
630  numbering_ = allocate();
631  std::copy( other.numbering_, other.numbering_ + capacity(), numbering_ );
632 
633  return *this;
634  }
635 
636  int size ( int cc ) const
637  {
638  assert( (cc >= codim()) && (cc <= dim) );
639  return (offset_[ cc+1 ] - offset_[ cc ]);
640  }
641 
642  int number ( int ii, int cc ) const
643  {
644  assert( (ii >= 0) && (ii < size( cc )) );
645  return numbering_[ offset_[ cc ] + ii ];
646  }
647 
648  const GeometryType &type () const { return type_; }
649 
650  void initialize ( unsigned int topologyId, int codim, unsigned int i )
651  {
652  const unsigned int subId = Impl::subTopologyId( topologyId, dim, codim, i );
653  type_ = GeometryType( subId, dim-codim );
654 
655  // compute offsets
656  for( int cc = 0; cc <= codim; ++cc )
657  offset_[ cc ] = 0;
658  for( int cc = codim; cc <= dim; ++cc )
659  offset_[ cc+1 ] = offset_[ cc ] + Impl::size( subId, dim-codim, cc-codim );
660 
661  // compute subnumbering
662  deallocate( numbering_ );
663  numbering_ = allocate();
664  for( int cc = codim; cc <= dim; ++cc )
665  Impl::subTopologyNumbering( topologyId, dim, codim, i, cc-codim, numbering_+offset_[ cc ], numbering_+offset_[ cc+1 ] );
666  }
667 
668  protected:
669  int codim () const { return dim - type().dim(); }
670 
671  unsigned int *allocate () { return (capacity() != 0 ? new unsigned int[ capacity() ] : nullptr); }
672  void deallocate ( unsigned int *ptr ) { delete[] ptr; }
673  unsigned int capacity () const { return offset_[ dim+1 ]; }
674 
675  private:
676  unsigned int *numbering_;
677  std::array< unsigned int, dim+2 > offset_;
678  GeometryType type_;
679  };
680 
681 
682  template< class ctype, int dim >
683  template< int codim >
684  struct ReferenceElementImplementation< ctype, dim >::CreateGeometries
685  {
686  template< int cc >
687  static typename ReferenceElements< ctype, dim-cc >::ReferenceElement
688  subRefElement( const ReferenceElementImplementation< ctype, dim > &refElement, int i, std::integral_constant< int, cc > )
689  {
690  return ReferenceElements< ctype, dim-cc >::general( refElement.type( i, cc ) );
691  }
692 
694  subRefElement( const ReferenceElementImplementation< ctype, dim > &refElement, int i, std::integral_constant< int, 0 > )
695  {
696  DUNE_UNUSED_PARAMETER(i);
697  return refElement;
698  }
699 
700  static void apply ( const ReferenceElementImplementation< ctype, dim > &refElement, GeometryTable &geometries )
701  {
702  const int size = refElement.size( codim );
703  std::vector< FieldVector< ctype, dim > > origins( size );
704  std::vector< FieldMatrix< ctype, dim - codim, dim > > jacobianTransposeds( size );
705  Impl::referenceEmbeddings( refElement.type().id(), dim, codim, &(origins[ 0 ]), &(jacobianTransposeds[ 0 ]) );
706 
707  std::get< codim >( geometries ).reserve( size );
708  for( int i = 0; i < size; ++i )
709  {
710  typename Codim< codim >::Geometry geometry( subRefElement( refElement, i, std::integral_constant< int, codim >() ), origins[ i ], jacobianTransposeds[ i ] );
711  std::get< codim >( geometries ).push_back( geometry );
712  }
713  }
714  };
715 
716 #endif // DOXYGEN
717 
718  } // namespace Geo
719 
720 } // namespace Dune
721 
722 #endif // #ifndef DUNE_GEOMETRY_REFERENCEELEMENTIMPLEMENTATION_HH
type.hh
A unique label for each type of element that can occur in a grid.
affinegeometry.hh
An implementation of the Geometry interface for affine geometries.
Dune::QuadratureType::size
@ size
Definition: quadraturerules.hh:88
Dune::Geo::ReferenceElements::general
static const ReferenceElement & general(const GeometryType &type)
get general reference elements
Definition: referenceelements.hh:196
referenceelement.hh
Dune::Geo::ReferenceElements::ReferenceElement
typename Container::ReferenceElement ReferenceElement
The reference element type.
Definition: referenceelements.hh:186
Dune
Definition: affinegeometry.hh:18
Dune::ReferenceElement
unspecified-type ReferenceElement
Returns the type of reference element for the argument type T.
Definition: referenceelements.hh:495