dune-grid  2.2.1
alugrid/3d/geometry.hh
Go to the documentation of this file.
1 #ifndef DUNE_ALU3DGRIDGEOMETRY_HH
2 #define DUNE_ALU3DGRIDGEOMETRY_HH
3 
4 // System includes
5 
6 // Dune includes
7 #include <dune/common/misc.hh>
9 
10 // Local includes
11 #include "alu3dinclude.hh"
12 #include "topology.hh"
13 #include "mappings.hh"
14 
15 namespace Dune
16 {
17 
18  // Forward declarations
19  template<int cd, int dim, class GridImp>
20  class ALU3dGridEntity;
21  template<int cd, class GridImp >
22  class ALU3dGridEntityPointer;
23  template<int mydim, int coorddim, class GridImp>
24  class ALU3dGridGeometry;
25  template< ALU3dGridElementType, class >
26  class ALU3dGrid;
27  class BilinearSurfaceMapping;
28  class TrilinearMapping;
29 
30  template< class GridImp >
31  class ALU3dGridIntersectionIterator;
32 
33  template <int cdim>
35  {
36  public:
37  typedef FieldVector<alu3d_ctype, cdim> CoordinateVectorType;
38 
39  static const signed char invalid = -1; // means geometry is not meaningfull
40  static const signed char updated = 0; // means the point values have been set
41  static const signed char buildmapping = 1; // means updated and mapping was build
42 
43  struct CoordVecCopy
44  {
45 
46  // copy coordinate vector from field vector or alu3d_ctype[cdim]
47  template <class CoordPtrType>
48  static inline void copy(const CoordPtrType& p,
50  {
51  assert( cdim == 3 );
52  c[0] = p[0];
53  c[1] = p[1];
54  c[2] = p[2];
55  }
56 
57  template <class CoordPtrType>
58  void update(const CoordPtrType&,
59  const CoordPtrType&,
60  const CoordPtrType&,
61  const CoordPtrType&,
62  const CoordPtrType&,
63  const CoordPtrType&,
64  const CoordPtrType&,
65  const CoordPtrType& ) const
66  {
67  DUNE_THROW(InvalidStateException,"This method should not be called!");
68  }
69 
70  template <class CoordPtrType>
71  void update(const CoordPtrType&,
72  const CoordPtrType&,
73  const CoordPtrType&,
74  const CoordPtrType& ) const
75  {
76  DUNE_THROW(InvalidStateException,"This method should not be called!");
77  }
78 
79  template <class CoordPtrType>
80  void update(const CoordPtrType&,
81  const CoordPtrType&,
82  const CoordPtrType& ) const
83  {
84  DUNE_THROW(InvalidStateException,"This method should not be called!");
85  }
86  };
87 
89  template <int dummy, int dim,
91  public:
92  // geometry implementation for edges and vertices
93  template <int dummy, int dim, ALU3dGridElementType eltype>
94  class GeometryImpl : public CoordVecCopy
95  {
96  using CoordVecCopy :: copy ;
97 
98  enum { corners_ = dim+1 };
100  typedef FieldMatrix<alu3d_ctype, corners_ , cdim> CoordinateMatrixType;
101 
102  typedef LinearMapping<cdim, dim> MappingType;
103 
104  CoordinateMatrixType coord_;
105  MappingType liMap_;
106  signed char status_;
107 
108  public:
109  using CoordVecCopy :: update ;
110 
111  GeometryImpl() : coord_() , liMap_() , status_( invalid ) {}
112  GeometryImpl(const GeometryImpl& other)
113  : coord_(other.coord_),
114  liMap_(other.liMap_),
115  status_(other.status_)
116  {}
117 
118  // return coordinate vector
119  inline const CoordinateVectorType& operator [] (const int i) const
120  {
121  assert( valid() );
122  assert( i>=0 && i<corners_ );
123  return coord_[i];
124  }
125 
126  inline MappingType& mapping()
127  {
128  assert( valid() );
129  if( status_ == buildmapping ) return liMap_;
130 
131  liMap_.buildMapping( coord_[0] );
132  status_ = buildmapping ;
133  return liMap_;
134  }
135 
136  // update vertex
137  template <class CoordPtrType>
138  inline void update(const CoordPtrType& p0)
139  {
140  assert( corners_ == 1 );
141  copy( p0, coord_[0] );
142  // we need to update the mapping
143  status_ = updated ;
144  }
145 
146  // set status to invalid
147  void invalidate () { status_ = invalid ; }
148 
149  // return true if geometry is valid
150  bool valid () const { return status_ != invalid ; }
151  };
152 
153  // geometry implementation for edges and vertices
154  template <int dummy, ALU3dGridElementType eltype>
155  class GeometryImpl<dummy,1,eltype> : public CoordVecCopy
156  {
157  using CoordVecCopy :: copy ;
158 
159  enum { dim = 1 };
160  enum { corners_ = dim+1 };
162  typedef FieldMatrix<alu3d_ctype, corners_ , cdim> CoordinateMatrixType;
163 
165 
166  // for edges use LinearMapping<cdim, 1> here that has all features
167  // implemented
168 
169  CoordinateMatrixType coord_;
170  MappingType liMap_;
171  signed char status_;
172 
173  public:
174  using CoordVecCopy :: update ;
175 
176  GeometryImpl() : coord_() , liMap_() , status_( invalid ) {}
177  GeometryImpl(const GeometryImpl& other)
178  : coord_(other.coord_),
179  liMap_(other.liMap_),
180  status_(other.status_)
181  {}
182 
183  // return coordinate vector
184  inline const CoordinateVectorType& operator [] (const int i) const
185  {
186  assert( valid() );
187  assert( i>=0 && i<corners_ );
188  return coord_[i];
189  }
190 
191  inline MappingType& mapping()
192  {
193  assert( valid() );
194  if( status_ == buildmapping ) return liMap_;
195 
196  liMap_.buildMapping( coord_[0], coord_[1] );
197  status_ = buildmapping ;
198  return liMap_;
199  }
200 
201  // update edge
202  template <class CoordPtrType>
203  inline void update(const CoordPtrType& p0,
204  const CoordPtrType& p1)
205  {
206  assert( corners_ == 2 );
207  copy( p0, coord_[0] );
208  copy( p1, coord_[1] );
209  status_ = updated;
210  }
211 
212  // set status to invalid
213  void invalidate () { status_ = invalid ; }
214 
215  // return true if geometry is valid
216  bool valid () const { return status_ != invalid ; }
217  };
218 
219  // geom impl for simplex faces (triangles)
220  template <int dummy>
221  class GeometryImpl<dummy,2, tetra> : public CoordVecCopy
222  {
223  using CoordVecCopy :: copy ;
224 
225  enum { corners_ = 3 };
226 
228  typedef FieldMatrix<alu3d_ctype, corners_ , cdim> CoordinateMatrixType;
230 
232  CoordinateMatrixType coord_;
233 
234  MappingType liMap_;
235  signed char status_;
236 
237  public:
238  using CoordVecCopy :: update ;
239 
240  // constructor
241  GeometryImpl() : coord_(), liMap_(), status_( invalid ) {}
242  // copy constructor
243  GeometryImpl(const GeometryImpl& other) :
244  coord_(other.coord_),
245  liMap_(other.liMap_),
246  status_(other.status_)
247  {}
248 
249  // return coordinate vector
250  inline const CoordinateVectorType& operator [] (const int i) const
251  {
252  assert( valid() );
253  assert( i>=0 && i<corners_ );
254  return coord_[i];
255  }
256 
257  // update geometry coordinates
258  template <class CoordPtrType>
259  inline void update(const CoordPtrType& p0,
260  const CoordPtrType& p1,
261  const CoordPtrType& p2)
262  {
263  copy(p0, coord_[0] );
264  copy(p1, coord_[1] );
265  copy(p2, coord_[2] );
266  status_ = updated;
267  }
268 
269  // return mapping (always up2date)
271  {
272  assert( valid() );
273  if( status_ == buildmapping ) return liMap_;
274 
275  liMap_.buildMapping( coord_[0], coord_[1], coord_[2] );
276  status_ = buildmapping ;
277  return liMap_;
278  }
279 
280  // set status to invalid
281  void invalidate () { status_ = invalid ; }
282 
283  // return true if geometry is valid
284  bool valid () const { return status_ != invalid ; }
285  };
286 
288  //
289  // hexa specializations
290  //
292 
293  // geom impl for cube faces (quadrilaterals)
294  template <int dummy>
295  class GeometryImpl<dummy,2, hexa> : public CoordVecCopy
296  {
297  using CoordVecCopy :: copy ;
298 
299  enum { corners_ = 4 };
300 
302  typedef FieldMatrix<alu3d_ctype, corners_ , cdim> CoordinateMatrixType;
304 
306  CoordinateMatrixType coord_;
307 
308  MappingType biMap_;
309  signed char status_;
310 
311  public:
312  using CoordVecCopy :: update ;
313 
314  // constructor
315  GeometryImpl() : coord_(), biMap_(), status_( invalid ) {}
316  // copy constructor
317  GeometryImpl(const GeometryImpl& other) :
318  coord_(other.coord_),
319  biMap_(other.biMap_),
320  status_(other.status_)
321  {}
322 
323  // return coordinate vector
324  inline const CoordinateVectorType& operator [] (const int i) const
325  {
326  assert( valid() );
327  assert( i>=0 && i<corners_ );
328  return coord_[i];
329  }
330 
331  // update geometry coordinates
332  template <class CoordPtrType>
333  inline void update(const CoordPtrType& p0,
334  const CoordPtrType& p1,
335  const CoordPtrType& p2,
336  const CoordPtrType& p3)
337  {
338  copy(p0, coord_[0] );
339  copy(p1, coord_[1] );
340  copy(p2, coord_[2] );
341  copy(p3, coord_[3] );
342  status_ = updated;
343  }
344 
345  // return mapping (always up2date)
347  {
348  assert( valid() );
349  if( status_ == buildmapping ) return biMap_;
350 
351  biMap_.buildMapping( coord_[0], coord_[1], coord_[2], coord_[3] );
352  status_ = buildmapping ;
353  return biMap_;
354  }
355 
356  // set status to invalid
357  void invalidate () { status_ = invalid ; }
358 
359  // return true if geometry is valid
360  bool valid () const { return status_ != invalid ; }
361  };
362 
363  // geometry impl for hexahedrons
364  template <int dummy>
365  class GeometryImpl<dummy,3, hexa> : public CoordVecCopy
366  {
367  enum { corners_ = 8 };
368 
370  typedef alu3d_ctype CoordPtrType[cdim];
371 
373  typedef FieldMatrix<alu3d_ctype, corners_ , cdim> CoordinateMatrixType;
374 
376 
377  // coordinate pointer vector
378  const alu3d_ctype* coordPtr_[ corners_ ];
379  MappingType triMap_;
380  CoordinateMatrixType* fatherCoord_;
381  signed char status_;
382 
383  public:
384  using CoordVecCopy :: update ;
385  using CoordVecCopy :: copy ;
386 
388  GeometryImpl() : triMap_(),
389  fatherCoord_(0),
390  status_( invalid )
391  {
392  // set initialize coord pointers
393  for( int i=0; i<corners_; ++i )
394  coordPtr_[ i ] = 0;
395  }
396 
397 
399  GeometryImpl(const GeometryImpl& other) :
400  triMap_(other.triMap_),
401  fatherCoord_(0),
402  status_(other.status_)
403  {
404  // copy coord pointers
405  for( int i=0; i<corners_; ++i )
406  coordPtr_[ i ] = other.coordPtr_[ i ];
407 
408  // if father coords are set, then reset coordPtr
409  if( other.fatherCoord_ )
410  {
411  fatherCoord_ = new CoordinateMatrixType(*other.fatherCoord_);
412  CoordinateMatrixType& coord = *fatherCoord_;
413  for(int i=0; i<corners_; ++i)
414  {
415  coordPtr_[i] = (&(coord[i][0]));
416  }
417  }
418  }
419 
420  // desctructor
422  {
423  if( fatherCoord_ ) delete fatherCoord_;
424  }
425 
426  const alu3d_ctype* point( const int i ) const
427  {
428  assert( valid() );
429  assert( i>=0 && i<corners_ );
430  assert( coordPtr_[i] );
431  return coordPtr_[ i ];
432  }
433 
434  // return coordinates
435  inline CoordinateVectorType operator [] (const int i) const
436  {
437  CoordinateVectorType coord ;
438  copy( point( i ), coord );
439  return coord ;
440  }
441 
442  // update geometry coordinates
443  inline void update(const CoordPtrType& p0,
444  const CoordPtrType& p1,
445  const CoordPtrType& p2,
446  const CoordPtrType& p3,
447  const CoordPtrType& p4,
448  const CoordPtrType& p5,
449  const CoordPtrType& p6,
450  const CoordPtrType& p7)
451  {
452  coordPtr_[0] = &p0[ 0 ];
453  coordPtr_[1] = &p1[ 0 ];
454  coordPtr_[2] = &p2[ 0 ];
455  coordPtr_[3] = &p3[ 0 ];
456  coordPtr_[4] = &p4[ 0 ];
457  coordPtr_[5] = &p5[ 0 ];
458  coordPtr_[6] = &p6[ 0 ];
459  coordPtr_[7] = &p7[ 0 ];
460  status_ = updated;
461  }
462 
463  // update geometry in father coordinates
464  template <class GeometryImp>
465  inline void updateInFather(const GeometryImp &fatherGeom ,
466  const GeometryImp &myGeom)
467  {
468  if( fatherCoord_ == 0 )
469  {
470  fatherCoord_ = new CoordinateMatrixType();
471  }
472 
473  CoordinateMatrixType& coord = *fatherCoord_;
474  // compute the local coordinates in father refelem
475  for(int i=0; i < myGeom.corners() ; ++i)
476  {
477  // calculate coordinate
478  coord[i] = fatherGeom.local( myGeom.corner( i ) );
479 
480  // set pointer
481  coordPtr_[i] = (&(coord[i][0]));
482 
483  // to avoid rounding errors
484  for(int j=0; j<cdim; ++j)
485  {
486  if ( coord[i][j] < 1e-16) coord[i][j] = 0.0;
487  }
488  }
489 
490  status_ = updated ;
491  }
492 
493  // return mapping (always up2date)
495  {
496  assert( valid() );
497  if( status_ == buildmapping ) return triMap_;
498 
499  triMap_.buildMapping( point( 0 ), point( 1 ), point( 2 ), point( 3 ),
500  point( 4 ), point( 5 ), point( 6 ), point( 7 ) );
501 
502  status_ = buildmapping;
503  return triMap_;
504  }
505 
506  // set status to invalid
507  void invalidate () { status_ = invalid ; }
508 
509  // return true if geometry is valid
510  bool valid () const { return status_ != invalid ; }
511  };
512 
513 
514  // geometry impl for hexahedrons
515  template <int dummy>
516  class GeometryImpl<dummy,3, tetra> : public CoordVecCopy
517  {
518  enum { corners_ = 4 };
519 
521  typedef alu3d_ctype CoordPtrType[cdim];
522 
524  typedef FieldMatrix<alu3d_ctype, corners_ , cdim> CoordinateMatrixType;
525 
527 
528  // coordinate pointer vector
529  const alu3d_ctype* coordPtr_[ corners_ ];
530  MappingType liMap_;
531  CoordinateMatrixType* fatherCoord_;
532  signed char status_;
533 
534  public:
535  using CoordVecCopy :: update ;
536  using CoordVecCopy :: copy ;
537 
538  // default constructor
539  GeometryImpl() : liMap_(),
540  fatherCoord_(0),
541  status_( invalid )
542  {
543  // set initialize coord pointers
544  for( int i=0; i<corners_; ++i )
545  coordPtr_[ i ] = 0;
546  }
547 
548  // copy constructor
549  GeometryImpl(const GeometryImpl& other) :
550  liMap_(other.liMap_),
551  fatherCoord_(0),
552  status_(other.status_)
553  {
554  // copy coord pointers
555  for( int i=0; i<corners_; ++i )
556  coordPtr_[ i ] = other.coordPtr_[ i ];
557 
558  // if father coords are set, then reset coordPtr
559  if( other.fatherCoord_ )
560  {
561  fatherCoord_ = new CoordinateMatrixType(*other.fatherCoord_);
562  CoordinateMatrixType& coord = *fatherCoord_;
563  for(int i=0; i<corners_; ++i)
564  {
565  coordPtr_[i] = (&(coord[i][0]));
566  }
567  }
568  }
569 
570  // destructor
572  {
573  if( fatherCoord_ ) delete fatherCoord_;
574  }
575 
576  const alu3d_ctype* point( const int i ) const
577  {
578  assert( valid() );
579  assert( i>=0 && i<corners_ );
580  assert( coordPtr_[ i ] );
581  return coordPtr_[ i ];
582  }
583 
584  // return coordinate vector
585  inline CoordinateVectorType operator [] (const int i) const
586  {
587  CoordinateVectorType coord ;
588  copy( point( i ), coord );
589  return coord ;
590  }
591 
592  // update geometry coordinates
593  inline void update(const CoordPtrType& p0,
594  const CoordPtrType& p1,
595  const CoordPtrType& p2,
596  const CoordPtrType& p3)
597  {
598  coordPtr_[0] = &p0[ 0 ];
599  coordPtr_[1] = &p1[ 0 ];
600  coordPtr_[2] = &p2[ 0 ];
601  coordPtr_[3] = &p3[ 0 ];
602  status_ = updated;
603  }
604 
605  // update geometry in father coordinates
606  template <class GeometryImp>
607  inline void updateInFather(const GeometryImp &fatherGeom ,
608  const GeometryImp & myGeom)
609  {
610  if( fatherCoord_ == 0 )
611  {
612  fatherCoord_ = new CoordinateMatrixType();
613  }
614 
615  CoordinateMatrixType& coord = *fatherCoord_;
616  // compute the local coordinates in father refelem
617  for(int i=0; i < myGeom.corners() ; ++i)
618  {
619  // calculate coordinate
620  coord[i] = fatherGeom.local( myGeom.corner( i ) );
621 
622  // set pointer
623  coordPtr_[i] = (&(coord[i][0]));
624 
625  // to avoid rounding errors
626  for(int j=0; j<cdim; ++j)
627  {
628  if ( coord[i][j] < 1e-16) coord[i][j] = 0.0;
629  }
630  }
631 
632  status_ = updated;
633  }
634 
635  // return mapping (always up2date)
637  {
638  assert( valid() );
639  if( status_ == buildmapping ) return liMap_;
640 
641  liMap_.buildMapping( point( 0 ), point( 1 ), point( 2 ), point( 3 ) );
642 
643  status_ = buildmapping;
644  return liMap_;
645  }
646 
647  // set status to invalid
648  void invalidate () { status_ = invalid ; }
649 
650  // return true if geometry is valid
651  bool valid () const { return status_ != invalid ; }
652  };
653  }; // end of class ALUGridGeometryImplementation
654 
655  template <int mydim, int cdim, class GridImp>
656  class ALU3dGridGeometry :
657  public GeometryDefaultImplementation<mydim, cdim, GridImp, ALU3dGridGeometry>
658  {
659  static const ALU3dGridElementType elementType = GridImp::elementType;
660 
661  typedef typename GridImp::MPICommunicatorType Comm;
662 
663  friend class ALU3dGridIntersectionIterator<GridImp>;
664 
665  typedef typename ALU3dImplTraits< elementType, Comm >::IMPLElementType IMPLElementType;
666  typedef typename ALU3dImplTraits< elementType, Comm >::GEOFaceType GEOFaceType;
667  typedef typename ALU3dImplTraits< elementType, Comm >::GEOEdgeType GEOEdgeType;
668  typedef typename ALU3dImplTraits< elementType, Comm >::GEOVertexType GEOVertexType;
669 
670  // interface types
671  typedef typename ALU3dImplTraits< elementType, Comm >::HFaceType HFaceType;
672  typedef typename ALU3dImplTraits< elementType, Comm >::HEdgeType HEdgeType;
673  typedef typename ALU3dImplTraits< elementType, Comm >::VertexType VertexType;
674 
677 
678  enum { corners_ = (elementType == hexa) ? Power_m_p<2,mydim>::power : mydim+1 };
679 
680  // type of specialized geometry implementation
681  typedef typename MyALUGridGeometryImplementation<cdim> ::
682  template GeometryImpl<0, mydim, elementType > GeometryImplType;
683 
684  public:
685  typedef typename GridImp :: ctype ctype;
686 
688  typedef FieldVector<ctype, mydim> LocalCoordinate;
689 
691  typedef FieldVector<ctype, cdim > GlobalCoordinate;
692 
694  typedef FieldMatrix<ctype,cdim,mydim> Jacobian;
695 
697  typedef FieldMatrix< ctype, mydim, cdim > JacobianTransposed;
698 
699  // type of coordinate matrix for faces
700  typedef FieldMatrix<ctype,
702 
706 
709  GeometryType type () const;
710 
712  int corners () const;
713 
715  const GlobalCoordinate& operator[] (int i) const;
716 
718  GlobalCoordinate corner (int i) const;
719 
722  GlobalCoordinate global (const LocalCoordinate& local) const;
723 
726  LocalCoordinate local (const GlobalCoordinate& global) const;
727 
729  ctype integrationElement (const LocalCoordinate& local) const;
730 
733  const Jacobian& jacobianInverseTransposed (const LocalCoordinate& local) const;
734 
736  const JacobianTransposed& jacobianTransposed (const LocalCoordinate& local) const;
737 
739  inline bool affine () const;
740 
742  ctype volume () const;
743 
744  //***********************************************************************
746  //***********************************************************************
748  bool buildGeom(const IMPLElementType & item);
749  bool buildGeom(const HFaceType & item, int twist, int faceNum);
750  bool buildGeom(const HEdgeType & item, int twist, int);
751  bool buildGeom(const VertexType & item, int twist, int);
752 
753  // this method is used by the intersection iterator
754  bool buildGeom(const FaceCoordinatesType& coords);
755 
756  // this method is used by the intersection iterator
757  template <class coord_t>
758  bool buildGeom(const coord_t& p0,
759  const coord_t& p1,
760  const coord_t& p2,
761  const coord_t& p3);
762 
763  // this method is used by the intersection iterator
764  template <class coord_t>
765  bool buildGeom(const coord_t& p0,
766  const coord_t& p1,
767  const coord_t& p2);
768 
770  template <class GeometryType>
771  bool buildGeomInFather(const GeometryType &fatherGeom , const GeometryType & myGeom);
772 
775  void print (std::ostream& ss) const;
776 
778  void invalidate () ;
779 
781  bool valid () const ;
782 
783  private:
784  // implementation of coord and mapping
785  mutable GeometryImplType geoImpl_;
786  // volume
787  mutable ctype volume_;
788  };
789 
790 } // end namespace Dune
791 
792 #include "geometry_imp.cc"
793 
794 #endif