dune-geometry  2.2.0
referenceelements.hh
Go to the documentation of this file.
1 // -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=8 sw=2 sts=2:
3 #ifndef DUNE_GEOMETRY_REFERENCEELEMENTS_HH
4 #define DUNE_GEOMETRY_REFERENCEELEMENTS_HH
5 
6 #include <dune/common/forloop.hh>
7 #include <dune/common/typetraits.hh>
8 
14 
15 namespace Dune
16 {
17 
18  // Internal Forward Declarations
19  // -----------------------------
20 
21  template< class ctype, int dim >
22  class GenericReferenceElementContainer;
23 
24 
25 
26  // GenericReferenceElement
27  // -----------------------
28 
45  template< class ctype, int dim >
47  {
49 
50  friend class GenericReferenceElementContainer< ctype, dim >;
51 
52  // make copy constructor private
54 
56 
58  {
59  ForLoop< Destroy, 0, dim >::apply( mappings_ );
60  integral_constant< int, 0 > codim0Variable;
61  if(mappings_[ codim0Variable ].size())
62  delete mappings_[ codim0Variable ][ 0 ];
63  }
64 
65  class SubEntityInfo;
66  template< class Topology > class CornerStorage;
67  template< class Topology > struct Initialize;
68  template< int codim > struct Destroy;
69 
70  struct GeometryTraits
71  : public GenericGeometry::DefaultGeometryTraits< ctype, dim, dim >
72  {
73  typedef GenericGeometry::DefaultGeometryTraits< ctype, dim, dim > Base;
74 
75  typedef typename Base::CoordTraits CoordTraits;
76 
77  template< class Topology >
78  struct Mapping
79  {
81  };
82 
83  struct Caching
84  {
89  };
90 
91  };
92 
93  public:
95  template< int codim >
96  struct Codim
97  {
99  typedef GenericGeometry::HybridMapping< dim-codim, GeometryTraits > Mapping;
100  };
101 
102  private:
104  template< int codim >
105  struct MappingArray
106  : public std::vector< typename Codim< codim >::Mapping * >
107  {};
108 
111 
112  std::vector< SubEntityInfo > info_[ dim+1 ];
113 
115  ctype volume_;
116  std::vector< FieldVector< ctype, dim > > volumeNormals_;
117 
119  MappingsTable mappings_;
120 
121  public:
126  int size ( int c ) const
127  {
128  assert( (c >= 0) && (c <= dim) );
129  return info_[ c ].size();
130  }
131 
143  int size ( int i, int c, int cc ) const
144  {
145  assert( (c >= 0) && (c <= dim) );
146  return info_[ c ][ i ].size( cc );
147  }
148 
162  int subEntity ( int i, int c, int ii, int cc ) const
163  {
164  assert( (c >= 0) && (c <= dim) );
165  return info_[ c ][ i ].number( ii, cc );
166  }
167 
177  const FieldVector< ctype, dim > &position( int i, int c ) const
178  {
179  assert( (c >= 0) && (c <= dim) );
180  return info_[ c ][ i ].position();
181  }
182 
190  bool checkInside ( const FieldVector< ctype, dim > &local ) const
191  {
192  return checkInside< 0 >( local, 0 );
193  }
194 
209  template< int codim >
210  bool checkInside ( const FieldVector< ctype, dim-codim > &local, int i ) const
211  {
212  return mapping< codim >( i ).checkInside( local );
213  }
214 
236  template< int codim >
237  FieldVector< ctype, dim >
238  global( const FieldVector< ctype, dim-codim > &local, int i, int c ) const
239  {
240  if( c != codim )
241  DUNE_THROW( Exception, "Local Coordinate Type does not correspond to codimension c." );
242  assert( c == codim );
243  return mapping< codim >( i ).global( local );
244  }
245 
264  template< int codim >
265  FieldVector< ctype, dim >
266  global( const FieldVector< ctype, dim-codim > &local, int i ) const
267  {
268  return mapping< codim >( i ).global( local );
269  }
270 
286  template< int codim >
287  typename Codim< codim >::Mapping &mapping( int i ) const
288  {
289  integral_constant< int, codim > codimVariable;
290  return *(mappings_[ codimVariable ][ i ]);
291  }
292 
301  const GeometryType &type ( int i, int c ) const
302  {
303  assert( (c >= 0) && (c <= dim) );
304  return info_[ c ][ i ].type();
305  }
306 
308  const GeometryType &type () const { return type( 0, 0 ); }
309 
310  unsigned int topologyId ( int i, int c ) const DUNE_DEPRECATED
311  {
312  assert( (c >= 0) && (c <= dim) );
313  return info_[ c ][ i ].topologyId();
314  }
315 
317  ctype volume () const
318  {
319  return volume_;
320  }
321 
329  const FieldVector< ctype, dim > &volumeOuterNormal ( int face ) const
330  {
331  assert( (face >= 0) && (face < int( volumeNormals_.size())) );
332  return volumeNormals_[ face ];
333  }
334 
341  template< class Topology >
343  {
344  dune_static_assert( (Topology::dimension == dim),
345  "Cannot initialize reference element for different dimension." );
346  typedef Initialize< Topology > Init;
348 
349  // set up subentities
350  integral_constant< int, 0 > codim0Variable;
351  mappings_[ codim0Variable ].resize( 1 );
352  mappings_[ codim0Variable ][ 0 ] = new VirtualMapping( codim0Variable );
353 
354  Dune::ForLoop< Init::template Codim, 0, dim >::apply( info_, mappings_ );
355 
356  // compute reference element volume
357  typedef GenericGeometry::ReferenceDomain< Topology > ReferenceDomain;
358  volume_ = ReferenceDomain::template volume< ctype >();
359 
360  // compute normals
361  volumeNormals_.resize( ReferenceDomain::numNormals );
362  for( unsigned int i = 0; i < ReferenceDomain::numNormals; ++i )
363  ReferenceDomain::integrationOuterNormal( i ,volumeNormals_[ i ] );
364  }
365  };
366 
367 
371  template< class ctype, int dim >
373  {
374  template< class Topology, int codim > struct Initialize
375  {
376  template< int subcodim > struct SubCodim;
377  };
378 
379  int codim_;
380  std::vector< int > numbering_[ dim+1 ];
381  FieldVector< ctype, dim > baryCenter_;
382  GeometryType type_;
383 
384  public:
385  int size ( int cc ) const
386  {
387  assert( (cc >= codim_) && (cc <= dim) );
388  return numbering_[ cc ].size();
389  }
390 
391  int number ( int ii, int cc ) const
392  {
393  assert( (cc >= codim_) && (cc <= dim) );
394  return numbering_[ cc ][ ii ];
395  }
396 
397  const FieldVector< ctype, dim > &position () const
398  {
399  return baryCenter_;
400  }
401 
402  const GeometryType &type () const
403  {
404  return type_;
405  }
406 
407  unsigned int topologyId () const DUNE_DEPRECATED
408  {
409  return type_.id();
410  }
411 
412  template< class Topology, unsigned int codim, unsigned int i >
413  void initialize ()
414  {
415  typedef Initialize< Topology, codim > Init;
417 
418  codim_ = codim;
419 
420  const unsigned int iVariable = i;
421  Dune::ForLoop< Init::template SubCodim, 0, dim-codim >::apply( iVariable, numbering_ );
422 
423  baryCenter_ = ctype( 0 );
424  static const unsigned int numCorners = size( dim );
425  for( unsigned int j = 0; j < numCorners; ++j )
426  {
427  FieldVector< ctype, dim > corner;
428  RefDomain::corner( number( j, dim ), corner );
429  baryCenter_ += corner;
430  }
431  baryCenter_ *= ctype( 1 ) / ctype( numCorners );
432 
434  type_ = GeometryType( SubTopology::id, SubTopology::dimension );
435  // type_ = GenericGeometry::DuneGeometryType< SubTopology, GeometryType::simplex >::type();
436  }
437  };
438 
439 
440  template< class ctype, int dim >
441  template< class Topology >
443  {
445 
446  public:
447  static const unsigned int size = Topology::numCorners;
448 
449  template< class SubTopology >
450  struct SubStorage
451  {
453  };
454 
455  explicit CornerStorage ( const integral_constant< int, 0 > & )
456  {
457  for( unsigned int i = 0; i < size; ++i )
458  RefDomain::corner( i, coords_[ i ] );
459  }
460 
461  template< class Mapping, unsigned int codim >
462  explicit
464  {
465  for( unsigned int i = 0; i < size; ++i )
466  coords_[ i ] = coords[ i ];
467  }
468 
469  const FieldVector< ctype, dim > &operator[] ( unsigned int i ) const
470  {
471  return coords_[ i ];
472  }
473 
474  private:
475  FieldVector< ctype, dim > coords_[ size ];
476  };
477 
478 
479  template< class ctype, int dim >
480  template< class Topology, int codim >
481  template< int subcodim >
482  struct GenericReferenceElement< ctype, dim >::SubEntityInfo::Initialize< Topology, codim >::SubCodim
483  {
486 
487  static void apply ( unsigned int i, std::vector< int > (&numbering)[ dim+1 ] )
488  {
489  const unsigned int size = SubSize::size( i );
490  numbering[ codim+subcodim ].resize( size );
491  for( unsigned int j = 0; j < size; ++j )
492  numbering[ codim+subcodim ][ j ] = SubNumbering::number( i, j );
493  }
494  };
495 
496 
497  template< class ctype, int dim >
498  template< class Topology >
499  struct GenericReferenceElement< ctype, dim >::Initialize
500  {
502 
503  typedef typename GenericReferenceElement::template Codim< 0 >::Mapping ReferenceMapping;
504 
505  template< int codim >
506  struct Codim
507  {
508  template< int i >
509  struct SubTopology
510  {
511  static void apply ( std::vector< SubEntityInfo > &info )
512  {
513  info[ i ].template initialize< Topology, codim, i >();
514  }
515  };
516 
517  static void
518  apply ( std::vector< SubEntityInfo > (&info)[ dim+1 ],
519  MappingsTable &mappings )
520  {
522  info[ codim ].resize( size );
523  Dune::ForLoop< SubTopology, 0, size-1 >::apply( info[ codim ] );
524 
525  if( codim > 0 )
526  {
527  integral_constant< int, 0 > codim0Variable;
528  const ReferenceMapping &refMapping = *(mappings[ codim0Variable ][ 0 ]);
529 
530  typedef typename GenericGeometry::MappingProvider< ReferenceMapping, codim > MappingProvider;
531 
532  integral_constant< int, codim > codimVariable;
533  mappings[ codimVariable ].resize( size );
534  for( unsigned int i = 0; i < size; ++i ) {
535  char* storage = new char[MappingProvider::maxMappingSize];
536  mappings[ codimVariable ][ i ] = refMapping.template trace< codim >( i, storage );
537  }
538  }
539  }
540  };
541  };
542 
543 
544 
545  template< class ctype, int dim >
546  template< int codim >
547  struct GenericReferenceElement< ctype, dim >::Destroy
548  {
549  static void apply ( MappingsTable &mappings )
550  {
551  if (codim > 0 )
552  {
553  integral_constant< int, codim > codimVariable;
554  for( size_t i = 0; i < mappings[ codimVariable ].size(); ++i ) {
555  typedef typename Codim<codim>::Mapping Mapping;
556  mappings[ codimVariable ][ i ]->~Mapping();
557  char* storage = (char*)mappings[ codimVariable ][ i ];
558  delete[](storage);
559  }
560  }
561  }
562  };
563 
564 
565 
566  // GenericReferenceElementContainer
567  // --------------------------------
568 
569  template< class ctype, int dim >
571  {
572  static const unsigned int numTopologies = (1u << dim);
573 
574  public:
576  typedef const value_type *const_iterator;
577 
579  {
580  ForLoop< Builder, 0, numTopologies-1 >::apply( values_ );
581  }
582 
583  const value_type &operator() ( const unsigned int topologyId ) const DUNE_DEPRECATED
584  {
585  return values_[ topologyId ];
586  }
587 
588  const value_type &operator() ( const GeometryType &type ) const
589  {
590  assert( type.dim() == dim );
591  return values_[ type.id() ];
592  }
593 
594  const value_type &simplex () const
595  {
597  }
598 
599  const value_type &cube () const
600  {
602  }
603 
604  const value_type &pyramid () const
605  {
607  }
608 
609  const value_type &prism () const
610  {
612  }
613 
614  const_iterator begin () const { return values_; }
615  const_iterator end () const { return values_ + numTopologies; }
616 
617  static const GenericReferenceElementContainer &instance () DUNE_DEPRECATED
618  {
620  return inst;
621  }
622 
623  private:
624  template< int topologyId >
625  struct Builder
626  {
627  static void apply ( value_type (&values)[ numTopologies ] )
628  {
629  typedef typename GenericGeometry::Topology< topologyId, dim >::type Topology;
630  values[ topologyId ].template initializeTopology< Topology >();
631  }
632  };
633 
634  value_type values_[ numTopologies ];
635  };
636 
637 
638 
639  // GenericReferenceElements
640  // ------------------------
641 
650  template< class ctype, int dim >
652  {
654 
657  general ( const GeometryType &type )
658  {
659  return container()( type );
660  }
661 
664  {
665  return container().simplex();
666  }
667 
670  {
671  return container().cube();
672  }
673 
674  static Iterator begin () { return container().begin(); }
675  static Iterator end () { return container().end(); }
676 
677  private:
678  static const GenericReferenceElementContainer< ctype, dim > &container ()
679  {
681  return container;
682  }
683  };
684 
685 } // namespace Dune
686 
687 #endif // #ifndef DUNE_GEOMETRY_REFERENCEELEMENTS_HH