dune-geometry  2.8.0
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 #include <bitset>
15 
16 #include <dune/common/fmatrix.hh>
17 #include <dune/common/fvector.hh>
18 #include <dune/common/hybridutilities.hh>
19 #include <dune/common/typetraits.hh>
20 #include <dune/common/iteratorrange.hh>
21 #include <dune/common/power.hh>
22 #include <dune/common/math.hh>
23 
26 #include <dune/geometry/type.hh>
27 
28 namespace Dune
29 {
30 
31  namespace Geo
32  {
33 
34 #ifndef DOXYGEN
35 
36  // Internal Forward Declarations
37  // -----------------------------
38 
39  namespace Impl
40  {
41  template< class ctype, int dim >
42  class ReferenceElementContainer;
43  }
44 
45  template< class ctype, int dim >
46  struct ReferenceElements;
47 
48 
49 
50  namespace Impl
51  {
52 
53  using Dune::Impl::isPrism;
54  using Dune::Impl::isPyramid;
55  using Dune::Impl::baseTopologyId;
56  using Dune::Impl::prismConstruction;
57  using Dune::Impl::pyramidConstruction;
58  using Dune::Impl::numTopologies;
59 
61  unsigned int size ( unsigned int topologyId, int dim, int codim );
62 
63 
64 
72  unsigned int subTopologyId ( unsigned int topologyId, int dim, int codim, unsigned int i );
73 
74 
75 
76  // subTopologyNumbering
77  // --------------------
78 
79  void subTopologyNumbering ( unsigned int topologyId, int dim, int codim, unsigned int i, int subcodim,
80  unsigned int *beginOut, unsigned int *endOut );
81 
82 
83 
84 
85  // checkInside
86  // -----------
87 
88  template< class ct, int cdim >
89  inline bool
90  checkInside ( unsigned int topologyId, int dim, const FieldVector< ct, cdim > &x, ct tolerance, ct factor = ct( 1 ) )
91  {
92  assert( (dim >= 0) && (dim <= cdim) );
93  assert( topologyId < numTopologies( dim ) );
94 
95  if( dim > 0 )
96  {
97  const ct baseFactor = (isPrism( topologyId, dim ) ? factor : factor - x[ dim-1 ]);
98  if( (x[ dim-1 ] > -tolerance) && (factor - x[ dim-1 ] > -tolerance) )
99  return checkInside< ct, cdim >( baseTopologyId( topologyId, dim ), dim-1, x, tolerance, baseFactor );
100  else
101  return false;
102  }
103  else
104  return true;
105  }
106 
107 
108 
109  // referenceCorners
110  // ----------------
111 
112  template< class ct, int cdim >
113  inline unsigned int
114  referenceCorners ( unsigned int topologyId, int dim, FieldVector< ct, cdim > *corners )
115  {
116  assert( (dim >= 0) && (dim <= cdim) );
117  assert( topologyId < numTopologies( dim ) );
118 
119  if( dim > 0 )
120  {
121  const unsigned int nBaseCorners
122  = referenceCorners( baseTopologyId( topologyId, dim ), dim-1, corners );
123  assert( nBaseCorners == size( baseTopologyId( topologyId, dim ), dim-1, dim-1 ) );
124  if( isPrism( topologyId, dim ) )
125  {
126  std::copy( corners, corners + nBaseCorners, corners + nBaseCorners );
127  for( unsigned int i = 0; i < nBaseCorners; ++i )
128  corners[ i+nBaseCorners ][ dim-1 ] = ct( 1 );
129  return 2*nBaseCorners;
130  }
131  else
132  {
133  corners[ nBaseCorners ] = FieldVector< ct, cdim >( ct( 0 ) );
134  corners[ nBaseCorners ][ dim-1 ] = ct( 1 );
135  return nBaseCorners+1;
136  }
137  }
138  else
139  {
140  *corners = FieldVector< ct, cdim >( ct( 0 ) );
141  return 1;
142  }
143  }
144 
145 
146 
147  // referenceVolume
148  // ---------------
149 
150  unsigned long referenceVolumeInverse ( unsigned int topologyId, int dim );
151 
152  template< class ct >
153  inline ct referenceVolume ( unsigned int topologyId, int dim )
154  {
155  return ct( 1 ) / ct( referenceVolumeInverse( topologyId, dim ) );
156  }
157 
158 
159 
160  // referenceOrigins
161  // ----------------
162 
163  template< class ct, int cdim >
164  inline unsigned int
165  referenceOrigins ( unsigned int topologyId, int dim, int codim, FieldVector< ct, cdim > *origins )
166  {
167  assert( (dim >= 0) && (dim <= cdim) );
168  assert( topologyId < numTopologies( dim ) );
169  assert( (codim >= 0) && (codim <= dim) );
170 
171  if( codim > 0 )
172  {
173  const unsigned int baseId = baseTopologyId( topologyId, dim );
174  if( isPrism( topologyId, dim ) )
175  {
176  const unsigned int n = (codim < dim ? referenceOrigins( baseId, dim-1, codim, origins ) : 0);
177  const unsigned int m = referenceOrigins( baseId, dim-1, codim-1, origins+n );
178  for( unsigned int i = 0; i < m; ++i )
179  {
180  origins[ n+m+i ] = origins[ n+i ];
181  origins[ n+m+i ][ dim-1 ] = ct( 1 );
182  }
183  return n+2*m;
184  }
185  else
186  {
187  const unsigned int m = referenceOrigins( baseId, dim-1, codim-1, origins );
188  if( codim == dim )
189  {
190  origins[ m ] = FieldVector< ct, cdim >( ct( 0 ) );
191  origins[ m ][ dim-1 ] = ct( 1 );
192  return m+1;
193  }
194  else
195  return m+referenceOrigins( baseId, dim-1, codim, origins+m );
196  }
197  }
198  else
199  {
200  origins[ 0 ] = FieldVector< ct, cdim >( ct( 0 ) );
201  return 1;
202  }
203  }
204 
205 
206 
207  // referenceEmbeddings
208  // -------------------
209 
210  template< class ct, int cdim, int mydim >
211  inline unsigned int
212  referenceEmbeddings ( unsigned int topologyId, int dim, int codim,
213  FieldVector< ct, cdim > *origins,
214  FieldMatrix< ct, mydim, cdim > *jacobianTransposeds )
215  {
216  assert( (0 <= codim) && (codim <= dim) && (dim <= cdim) );
217  assert( (dim - codim <= mydim) && (mydim <= cdim) );
218  assert( topologyId < numTopologies( dim ) );
219 
220  if( codim > 0 )
221  {
222  const unsigned int baseId = baseTopologyId( topologyId, dim );
223  if( isPrism( topologyId, dim ) )
224  {
225  const unsigned int n = (codim < dim ? referenceEmbeddings( baseId, dim-1, codim, origins, jacobianTransposeds ) : 0);
226  for( unsigned int i = 0; i < n; ++i )
227  jacobianTransposeds[ i ][ dim-codim-1 ][ dim-1 ] = ct( 1 );
228 
229  const unsigned int m = referenceEmbeddings( baseId, dim-1, codim-1, origins+n, jacobianTransposeds+n );
230  std::copy( origins+n, origins+n+m, origins+n+m );
231  std::copy( jacobianTransposeds+n, jacobianTransposeds+n+m, jacobianTransposeds+n+m );
232  for( unsigned int i = 0; i < m; ++i )
233  origins[ n+m+i ][ dim-1 ] = ct( 1 );
234 
235  return n+2*m;
236  }
237  else
238  {
239  const unsigned int m = referenceEmbeddings( baseId, dim-1, codim-1, origins, jacobianTransposeds );
240  if( codim == dim )
241  {
242  origins[ m ] = FieldVector< ct, cdim >( ct( 0 ) );
243  origins[ m ][ dim-1 ] = ct( 1 );
244  jacobianTransposeds[ m ] = FieldMatrix< ct, mydim, cdim >( ct( 0 ) );
245  return m+1;
246  }
247  else
248  {
249  const unsigned int n = referenceEmbeddings( baseId, dim-1, codim, origins+m, jacobianTransposeds+m );
250  for( unsigned int i = 0; i < n; ++i )
251  {
252  for( int k = 0; k < dim-1; ++k )
253  jacobianTransposeds[ m+i ][ dim-codim-1 ][ k ] = -origins[ m+i ][ k ];
254  jacobianTransposeds[ m+i ][ dim-codim-1 ][ dim-1 ] = ct( 1 );
255  }
256  return m+n;
257  }
258  }
259  }
260  else
261  {
262  origins[ 0 ] = FieldVector< ct, cdim >( ct( 0 ) );
263  jacobianTransposeds[ 0 ] = FieldMatrix< ct, mydim, cdim >( ct( 0 ) );
264  for( int k = 0; k < dim; ++k )
265  jacobianTransposeds[ 0 ][ k ][ k ] = ct( 1 );
266  return 1;
267  }
268  }
269 
270 
271 
272  // referenceIntegrationOuterNormals
273  // --------------------------------
274 
275  template< class ct, int cdim >
276  inline unsigned int
277  referenceIntegrationOuterNormals ( unsigned int topologyId, int dim,
278  const FieldVector< ct, cdim > *origins,
279  FieldVector< ct, cdim > *normals )
280  {
281  assert( (dim > 0) && (dim <= cdim) );
282  assert( topologyId < numTopologies( dim ) );
283 
284  if( dim > 1 )
285  {
286  const unsigned int baseId = baseTopologyId( topologyId, dim );
287  if( isPrism( topologyId, dim ) )
288  {
289  const unsigned int numBaseFaces
290  = referenceIntegrationOuterNormals( baseId, dim-1, origins, normals );
291 
292  for( unsigned int i = 0; i < 2; ++i )
293  {
294  normals[ numBaseFaces+i ] = FieldVector< ct, cdim >( ct( 0 ) );
295  normals[ numBaseFaces+i ][ dim-1 ] = ct( 2*int( i )-1 );
296  }
297 
298  return numBaseFaces+2;
299  }
300  else
301  {
302  normals[ 0 ] = FieldVector< ct, cdim >( ct( 0 ) );
303  normals[ 0 ][ dim-1 ] = ct( -1 );
304 
305  const unsigned int numBaseFaces
306  = referenceIntegrationOuterNormals( baseId, dim-1, origins+1, normals+1 );
307  for( unsigned int i = 1; i <= numBaseFaces; ++i )
308  normals[ i ][ dim-1 ] = normals[ i ]*origins[ i ];
309 
310  return numBaseFaces+1;
311  }
312  }
313  else
314  {
315  for( unsigned int i = 0; i < 2; ++i )
316  {
317  normals[ i ] = FieldVector< ct, cdim >( ct( 0 ) );
318  normals[ i ][ 0 ] = ct( 2*int( i )-1 );
319  }
320 
321  return 2;
322  }
323  }
324 
325  template< class ct, int cdim >
326  inline unsigned int
327  referenceIntegrationOuterNormals ( unsigned int topologyId, int dim,
328  FieldVector< ct, cdim > *normals )
329  {
330  assert( (dim > 0) && (dim <= cdim) );
331 
332  FieldVector< ct, cdim > *origins
333  = new FieldVector< ct, cdim >[ size( topologyId, dim, 1 ) ];
334  referenceOrigins( topologyId, dim, 1, origins );
335 
336  const unsigned int numFaces
337  = referenceIntegrationOuterNormals( topologyId, dim, origins, normals );
338  assert( numFaces == size( topologyId, dim, 1 ) );
339 
340  delete[] origins;
341 
342  return numFaces;
343  }
344 
345  } // namespace Impl
346 
347 
348 
349  // ReferenceElement
350  // ----------------
351 
370  template< class ctype_, int dim >
371  class ReferenceElementImplementation
372  {
373 
374  public:
375 
377  using ctype = ctype_;
378 
380  using CoordinateField = ctype;
381 
383  using Coordinate = Dune::FieldVector<ctype,dim>;
384 
386  static constexpr int dimension = dim;
387 
389  typedef ctype Volume;
390 
391  private:
392 
393  friend class Impl::ReferenceElementContainer< ctype, dim >;
394 
395  struct SubEntityInfo;
396 
397  template< int codim > struct CreateGeometries;
398 
399  public:
401  template< int codim >
402  struct Codim
403  {
405  typedef AffineGeometry< ctype, dim-codim, dim > Geometry;
406  };
407 
408  // ReferenceElement cannot be copied.
409  ReferenceElementImplementation ( const ReferenceElementImplementation& ) = delete;
410 
411  // ReferenceElementImplementation cannot be copied.
412  ReferenceElementImplementation& operator= ( const ReferenceElementImplementation& ) = delete;
413 
414  // ReferenceElementImplementation is default-constructible (required for storage in std::array)
415  ReferenceElementImplementation () = default;
416 
421  int size ( int c ) const
422  {
423  assert( (c >= 0) && (c <= dim) );
424  return info_[ c ].size();
425  }
426 
438  int size ( int i, int c, int cc ) const
439  {
440  assert( (i >= 0) && (i < size( c )) );
441  return info_[ c ][ i ].size( cc );
442  }
443 
457  int subEntity ( int i, int c, int ii, int cc ) const
458  {
459  assert( (i >= 0) && (i < size( c )) );
460  return info_[ c ][ i ].number( ii, cc );
461  }
462 
478  auto subEntities ( int i, int c, int cc ) const
479  {
480  assert( (i >= 0) && (i < size( c )) );
481  return info_[ c ][ i ].numbers( cc );
482  }
483 
492  const GeometryType &type ( int i, int c ) const
493  {
494  assert( (i >= 0) && (i < size( c )) );
495  return info_[ c ][ i ].type();
496  }
497 
499  const GeometryType &type () const { return type( 0, 0 ); }
500 
510  const Coordinate &position( int i, int c ) const
511  {
512  assert( (c >= 0) && (c <= dim) );
513  return baryCenters_[ c ][ i ];
514  }
515 
523  bool checkInside ( const Coordinate &local ) const
524  {
525  const ctype tolerance = ctype( 64 ) * std::numeric_limits< ctype >::epsilon();
526  return Impl::template checkInside< ctype, dim >( type().id(), dim, local, tolerance );
527  }
528 
540  template< int codim >
541  typename Codim< codim >::Geometry geometry ( int i ) const
542  {
543  return std::get< codim >( geometries_ )[ i ];
544  }
545 
547  Volume volume () const
548  {
549  return volume_;
550  }
551 
559  const Coordinate &integrationOuterNormal ( int face ) const
560  {
561  assert( (face >= 0) && (face < int( integrationNormals_.size() )) );
562  return integrationNormals_[ face ];
563  }
564 
565  private:
566  void initialize ( unsigned int topologyId )
567  {
568  assert( topologyId < Impl::numTopologies( dim ) );
569 
570  // set up subentities
571  for( int codim = 0; codim <= dim; ++codim )
572  {
573  const unsigned int size = Impl::size( topologyId, dim, codim );
574  info_[ codim ].resize( size );
575  for( unsigned int i = 0; i < size; ++i )
576  info_[ codim ][ i ].initialize( topologyId, codim, i );
577  }
578 
579  // compute corners
580  const unsigned int numVertices = size( dim );
581  baryCenters_[ dim ].resize( numVertices );
582  Impl::referenceCorners( topologyId, dim, &(baryCenters_[ dim ][ 0 ]) );
583 
584  // compute barycenters
585  for( int codim = 0; codim < dim; ++codim )
586  {
587  baryCenters_[ codim ].resize( size(codim) );
588  for( int i = 0; i < size( codim ); ++i )
589  {
590  baryCenters_[ codim ][ i ] = Coordinate( ctype( 0 ) );
591  const unsigned int numCorners = size( i, codim, dim );
592  for( unsigned int j = 0; j < numCorners; ++j )
593  baryCenters_[ codim ][ i ] += baryCenters_[ dim ][ subEntity( i, codim, j, dim ) ];
594  baryCenters_[ codim ][ i ] *= ctype( 1 ) / ctype( numCorners );
595  }
596  }
597 
598  // compute reference element volume
599  volume_ = Impl::template referenceVolume< ctype >( topologyId, dim );
600 
601  // compute integration outer normals
602  if( dim > 0 )
603  {
604  integrationNormals_.resize( size( 1 ) );
605  Impl::referenceIntegrationOuterNormals( topologyId, dim, &(integrationNormals_[ 0 ]) );
606  }
607 
608  // set up geometries
609  Hybrid::forEach( std::make_index_sequence< dim+1 >{}, [ & ]( auto i ){ CreateGeometries< i >::apply( *this, geometries_ ); } );
610  }
611 
612  template< int... codim >
613  static std::tuple< std::vector< typename Codim< codim >::Geometry >... >
614  makeGeometryTable ( std::integer_sequence< int, codim... > );
615 
617  typedef decltype( makeGeometryTable( std::make_integer_sequence< int, dim+1 >() ) ) GeometryTable;
618 
620  ctype volume_;
621 
622  std::vector< Coordinate > baryCenters_[ dim+1 ];
623  std::vector< Coordinate > integrationNormals_;
624 
626  GeometryTable geometries_;
627 
628  std::vector< SubEntityInfo > info_[ dim+1 ];
629  };
630 
632  template< class ctype, int dim >
633  struct ReferenceElementImplementation< ctype, dim >::SubEntityInfo
634  {
635  // Compute upper bound for the number of subsentities.
636  // If someone knows an explicit formal feel free to
637  // implement it here.
638  static constexpr std::size_t maxSubEntityCount()
639  {
640  std::size_t maxCount=0;
641  for(std::size_t codim=0; codim<=dim; ++codim)
642  maxCount = std::max(maxCount, binomial(std::size_t(dim),codim)*(1 << codim));
643  return maxCount;
644  }
645 
646  using SubEntityFlags = std::bitset<maxSubEntityCount()>;
647 
648  class SubEntityRange
649  : public Dune::IteratorRange<const unsigned int*>
650  {
651  using Base = typename Dune::IteratorRange<const unsigned int*>;
652 
653  public:
654 
655  using iterator = Base::iterator;
656  using const_iterator = Base::const_iterator;
657 
658  SubEntityRange(const iterator& begin, const iterator& end, const SubEntityFlags& contains) :
659  Base(begin, end),
660  containsPtr_(&contains),
661  size_(end-begin)
662  {}
663 
664  SubEntityRange() :
665  Base(),
666  containsPtr_(nullptr),
667  size_(0)
668  {}
669 
670  std::size_t size() const
671  {
672  return size_;
673  }
674 
675  bool contains(std::size_t i) const
676  {
677  return (*containsPtr_)[i];
678  }
679 
680  private:
681  const SubEntityFlags* containsPtr_;
682  std::size_t size_;
683  std::size_t offset_;
684  };
685 
686  using NumberRange = typename Dune::IteratorRange<const unsigned int*>;
687 
688  SubEntityInfo ()
689  : numbering_( nullptr )
690  {
691  std::fill( offset_.begin(), offset_.end(), 0 );
692  }
693 
694  SubEntityInfo ( const SubEntityInfo &other )
695  : offset_( other.offset_ ),
696  type_( other.type_ ),
697  containsSubentity_( other.containsSubentity_ )
698  {
699  numbering_ = allocate();
700  std::copy( other.numbering_, other.numbering_ + capacity(), numbering_ );
701  }
702 
703  ~SubEntityInfo () { deallocate( numbering_ ); }
704 
705  const SubEntityInfo &operator= ( const SubEntityInfo &other )
706  {
707  type_ = other.type_;
708  offset_ = other.offset_;
709 
710  deallocate( numbering_ );
711  numbering_ = allocate();
712  std::copy( other.numbering_, other.numbering_ + capacity(), numbering_ );
713 
714  containsSubentity_ = other.containsSubentity_;
715 
716  return *this;
717  }
718 
719  int size ( int cc ) const
720  {
721  assert( (cc >= 0) && (cc <= dim) );
722  return (offset_[ cc+1 ] - offset_[ cc ]);
723  }
724 
725  int number ( int ii, int cc ) const
726  {
727  assert( (ii >= 0) && (ii < size( cc )) );
728  return numbering_[ offset_[ cc ] + ii ];
729  }
730 
731  auto numbers ( int cc ) const
732  {
733  return SubEntityRange( numbering_ + offset_[ cc ], numbering_ + offset_[ cc+1 ], containsSubentity_[cc]);
734  }
735 
736  const GeometryType &type () const { return type_; }
737 
738  void initialize ( unsigned int topologyId, int codim, unsigned int i )
739  {
740  const unsigned int subId = Impl::subTopologyId( topologyId, dim, codim, i );
741  type_ = GeometryType( subId, dim-codim );
742 
743  // compute offsets
744  for( int cc = 0; cc <= codim; ++cc )
745  offset_[ cc ] = 0;
746  for( int cc = codim; cc <= dim; ++cc )
747  offset_[ cc+1 ] = offset_[ cc ] + Impl::size( subId, dim-codim, cc-codim );
748 
749  // compute subnumbering
750  deallocate( numbering_ );
751  numbering_ = allocate();
752  for( int cc = codim; cc <= dim; ++cc )
753  Impl::subTopologyNumbering( topologyId, dim, codim, i, cc-codim, numbering_+offset_[ cc ], numbering_+offset_[ cc+1 ] );
754 
755  // initialize containsSubentity lookup-table
756  for(std::size_t cc=0; cc<= dim; ++cc)
757  {
758  containsSubentity_[cc].reset();
759  for(std::size_t idx=0; idx<std::size_t(size(cc)); ++idx)
760  containsSubentity_[cc][number(idx,cc)] = true;
761  }
762  }
763 
764  protected:
765  int codim () const { return dim - type().dim(); }
766 
767  unsigned int *allocate () { return (capacity() != 0 ? new unsigned int[ capacity() ] : nullptr); }
768  void deallocate ( unsigned int *ptr ) { delete[] ptr; }
769  unsigned int capacity () const { return offset_[ dim+1 ]; }
770 
771  private:
772  unsigned int *numbering_;
773  std::array< unsigned int, dim+2 > offset_;
774  GeometryType type_;
775  std::array< SubEntityFlags, dim+1> containsSubentity_;
776  };
777 
778 
779  template< class ctype, int dim >
780  template< int codim >
781  struct ReferenceElementImplementation< ctype, dim >::CreateGeometries
782  {
783  template< int cc >
784  static typename ReferenceElements< ctype, dim-cc >::ReferenceElement
785  subRefElement( const ReferenceElementImplementation< ctype, dim > &refElement, int i, std::integral_constant< int, cc > )
786  {
787  return ReferenceElements< ctype, dim-cc >::general( refElement.type( i, cc ) );
788  }
789 
791  subRefElement(const ReferenceElementImplementation< ctype, dim > &refElement,
792  [[maybe_unused]] int i, std::integral_constant<int, 0>)
793  {
794  return refElement;
795  }
796 
797  static void apply ( const ReferenceElementImplementation< ctype, dim > &refElement, GeometryTable &geometries )
798  {
799  const int size = refElement.size( codim );
800  std::vector< FieldVector< ctype, dim > > origins( size );
801  std::vector< FieldMatrix< ctype, dim - codim, dim > > jacobianTransposeds( size );
802  Impl::referenceEmbeddings( refElement.type().id(), dim, codim, &(origins[ 0 ]), &(jacobianTransposeds[ 0 ]) );
803 
804  std::get< codim >( geometries ).reserve( size );
805  for( int i = 0; i < size; ++i )
806  {
807  typename Codim< codim >::Geometry geometry( subRefElement( refElement, i, std::integral_constant< int, codim >() ), origins[ i ], jacobianTransposeds[ i ] );
808  std::get< codim >( geometries ).push_back( geometry );
809  }
810  }
811  };
812 
813 #endif // DOXYGEN
814 
815  } // namespace Geo
816 
817 } // namespace Dune
818 
819 #endif // #ifndef DUNE_GEOMETRY_REFERENCEELEMENTIMPLEMENTATION_HH
An implementation of the Geometry interface for affine geometries.
A unique label for each type of element that can occur in a grid.
unspecified-type ReferenceElement
Returns the type of reference element for the argument type T.
Definition: referenceelements.hh:495
Definition: affinegeometry.hh:19
@ size
Definition: quadraturerules.hh:143
int binomial(int upper, int lower)
calculate
Definition: simplex.cc:295
typename Container::ReferenceElement ReferenceElement
The reference element type.
Definition: referenceelements.hh:186
static const ReferenceElement & general(const GeometryType &type)
get general reference elements
Definition: referenceelements.hh:196