dune-grid  2.2.0
alugrid/2d/geometry.hh
Go to the documentation of this file.
1 #ifndef DUNE_ALU2DGRIDGEOMETRY_HH
2 #define DUNE_ALU2DGRIDGEOMETRY_HH
3 
4 // Dune includes
5 #include <dune/common/misc.hh>
7 #include <dune/geometry/genericgeometry/topologytypes.hh>
8 
11 
12 namespace Dune
13 {
14 
15  // Forward declarations
16  template<int cd, int dim, class GridImp>
17  class ALU2dGridEntity;
18  template<int cd, class GridImp >
19  class ALU2dGridEntityPointer;
20  template<int mydim, int cdim, class GridImp>
21  class ALU2dGridGeometry;
22  template< int dim, int dimworld, ALU2DSPACE ElementType eltype >
23  class ALU2dGrid;
24 
25 
26  template< int mydim, int cdim, ALU2DSPACE ElementType eltype >
27  class MyALU2dGridGeometryImpl;
28 
29  // geometry implementation for vertices
30  template< int cdim, ALU2DSPACE ElementType eltype >
31  class MyALU2dGridGeometryImpl< 0, cdim, eltype >
32  {
34 
35  typedef typename MappingType::ctype ctype;
36 
37  typedef typename MappingType::map_t map_t;
38  typedef typename MappingType::world_t world_t;
39 
40  typedef typename MappingType::matrix_t matrix_t;
41  typedef typename MappingType::inv_t inv_t;
42 
43  MappingType mapping_;
44  bool valid_ ;
45 
46  const MappingType& mapping() const
47  {
48  assert( valid() );
49  return mapping_;
50  }
51 
52  public:
53  MyALU2dGridGeometryImpl() : mapping_(), valid_( false ) {}
54 
55  // returns true if goemetry info is valid
56  bool valid () const { return valid_; }
57 
58  // reset geometry status
59  void invalidate() { valid_ = false ; }
60 
61  bool affine() const
62  {
63  assert( valid() );
64  return mapping().affine();
65  }
66 
67  int corners () const
68  {
69  return 1;
70  }
71 
72  GeometryType type () const
73  {
74  return GeometryType(
75  (eltype == ALU2DSPACE triangle ?
76  GenericGeometry :: SimplexTopology< 0 > :: type :: id :
77  GenericGeometry :: CubeTopology < 0 > :: type :: id),
78  0 );
79  }
80 
81  void map2world ( const map_t &m, world_t &w ) const
82  {
83  return mapping().map2world( m, w );
84  }
85 
86  void world2map ( const world_t &w, map_t &m ) const
87  {
88  return mapping().world2map( w, m );
89  }
90 
91  const matrix_t &jacobianTransposed ( const map_t &m ) const
92  {
93  return mapping().jacobianTransposed( m );
94  }
95 
96  const inv_t &jacobianInverseTransposed ( const map_t &m ) const
97  {
98  return mapping().jacobianInverseTransposed( m );
99  }
100 
101  ctype det ( const map_t &m ) const
102  {
103  return mapping().det( m );
104  }
105 
106  // update geometry coordinates
107  template< class Vector >
108  void update ( const Vector &p0 )
109  {
110  mapping_.buildMapping( p0 );
111  valid_ = true ;
112  }
113  };
114 
115  // geometry implementation for lines
116  template< int cdim, ALU2DSPACE ElementType eltype >
117  class MyALU2dGridGeometryImpl< 1, cdim, eltype >
118  {
119  static const int ncorners = 2;
120 
122 
123  typedef typename MappingType::ctype ctype;
124 
125  typedef typename MappingType::map_t map_t;
126  typedef typename MappingType::world_t world_t;
127 
128  typedef typename MappingType::matrix_t matrix_t;
129  typedef typename MappingType::inv_t inv_t;
130 
131  MappingType mapping_;
132  bool valid_;
133 
134  const MappingType& mapping() const
135  {
136  assert( valid() );
137  return mapping_;
138  }
139 
140  public:
141  MyALU2dGridGeometryImpl() : mapping_(), valid_( false ) {}
142 
143  // returns true if goemetry info is valid
144  bool valid () const { return valid_; }
145 
146  // reset geometry status
147  void invalidate() { valid_ = false ; }
148 
149  bool affine() const
150  {
151  return mapping().affine();
152  }
153 
154  int corners () const
155  {
156  return ncorners;
157  }
158 
159  GeometryType type () const
160  {
161  return GeometryType(
162  (eltype == ALU2DSPACE triangle ?
163  GenericGeometry :: SimplexTopology< 1 > :: type :: id :
164  GenericGeometry :: CubeTopology < 1 > :: type :: id),
165  1 );
166  }
167 
168  void map2world ( const map_t &m, world_t &w ) const
169  {
170  return mapping().map2world( m, w );
171  }
172 
173  void world2map ( const world_t &w, map_t &m ) const
174  {
175  return mapping().world2map( w, m );
176  }
177 
178  const matrix_t &jacobianTransposed ( const map_t &m ) const
179  {
180  return mapping().jacobianTransposed( m );
181  }
182 
183  const inv_t &jacobianInverseTransposed ( const map_t &m ) const
184  {
185  return mapping().jacobianInverseTransposed( m );
186  }
187 
188  ctype det ( const map_t &m ) const
189  {
190  return mapping().det( m );
191  }
192 
193  // update geometry in father coordinates
194  template< class Geo, class LocalGeo >
195  void updateLocal ( const Geo &geo, const LocalGeo &localGeo )
196  {
197  assert( localGeo.corners() == ncorners );
198  // compute the local coordinates in father refelem
199  FieldMatrix< alu2d_ctype, ncorners, cdim > coord;
200  for( int i = 0; i < ncorners; ++i )
201  {
202  // calculate coordinate
203  coord[ i ] = geo.local( localGeo.corner( i ) );
204  // to avoid rounding errors
205  for( int j = 0; j < cdim; ++j )
206  coord[ i ][ j ] = (coord[ i ][ j ] < 1e-14 ? 0 : coord[ i ][ j ]);
207  }
208  mapping_.buildMapping( coord[ 0 ], coord[ 1 ] );
209  valid_ = true ;
210  }
211 
212  // update geometry coordinates
213  template< class Vector >
214  void update ( const Vector &p0, const Vector &p1 )
215  {
216  mapping_.buildMapping( p0, p1 );
217  valid_ = true ;
218  }
219  };
220 
221  // geometry implementation for triangles
222  template< int cdim >
223  class MyALU2dGridGeometryImpl< 2, cdim, ALU2DSPACE triangle >
224  {
225  static const int ncorners = 3;
226 
228 
229  typedef typename MappingType::ctype ctype;
230 
231  typedef typename MappingType::map_t map_t;
232  typedef typename MappingType::world_t world_t;
233 
234  typedef typename MappingType::matrix_t matrix_t;
235  typedef typename MappingType::inv_t inv_t;
236 
237  MappingType mapping_;
238  bool valid_;
239 
240  const MappingType& mapping() const
241  {
242  assert( valid() );
243  return mapping_;
244  }
245 
246  public:
247  MyALU2dGridGeometryImpl() : mapping_(), valid_( false ) {}
248 
249  // returns true if goemetry info is valid
250  bool valid () const { return valid_; }
251 
252  // reset geometry status
253  void invalidate() { valid_ = false ; }
254 
255  bool affine () const
256  {
257  return mapping().affine();
258  }
259 
260  int corners () const
261  {
262  return ncorners;
263  }
264 
265  GeometryType type () const
266  {
267  return GeometryType( GenericGeometry :: SimplexTopology< 2 > :: type :: id , 2 );
268  }
269 
270  void map2world ( const map_t &m, world_t &w ) const
271  {
272  return mapping().map2world( m, w );
273  }
274 
275  void world2map ( const world_t &w, map_t &m ) const
276  {
277  return mapping().world2map( w, m );
278  }
279 
280  const matrix_t &jacobianTransposed ( const map_t &m ) const
281  {
282  return mapping().jacobianTransposed( m );
283  }
284 
285  const inv_t &jacobianInverseTransposed ( const map_t &m ) const
286  {
287  return mapping().jacobianInverseTransposed( m );
288  }
289 
290  ctype det ( const map_t &m ) const
291  {
292  return mapping().det( m );
293  }
294 
295  // update geometry in father coordinates
296  template< class Geo, class LocalGeo >
297  void updateLocal ( const Geo &geo, const LocalGeo &localGeo )
298  {
299  assert( localGeo.corners() == ncorners );
300  // compute the local coordinates in father refelem
301  FieldMatrix< alu2d_ctype, ncorners, cdim > coord;
302  for( int i = 0; i < ncorners; ++i )
303  {
304  // calculate coordinate
305  coord[ i ] = geo.local( localGeo.corner( i ) );
306  // to avoid rounding errors
307  for( int j = 0; j < cdim; ++j )
308  coord[ i ][ j ] = (coord[ i ][ j ] < 1e-14 ? 0 : coord[ i ][ j ]);
309  }
310  mapping_.buildMapping( coord[ 0 ], coord[ 1 ], coord[ 2 ] );
311  valid_ = true ;
312  }
313 
314  template< class HElement >
315  void update ( const HElement &item )
316  {
317  mapping_.buildMapping( item.getVertex( 0 )->coord(), item.getVertex( 1 )->coord(),
318  item.getVertex( 2 )->coord() );
319  valid_ = true ;
320  }
321  };
322 
323  // geometry implementation for quadrilaterals
324  template< int cdim >
325  class MyALU2dGridGeometryImpl< 2, cdim, ALU2DSPACE quadrilateral >
326  {
327  static const int ncorners = 4;
328 
330 
331  typedef typename MappingType::ctype ctype;
332 
333  typedef typename MappingType::map_t map_t;
334  typedef typename MappingType::world_t world_t;
335 
336  typedef typename MappingType::matrix_t matrix_t;
337  typedef typename MappingType::inv_t inv_t;
338 
339  MappingType mapping_;
340  bool valid_ ;
341 
342  const MappingType& mapping() const
343  {
344  assert( valid() );
345  return mapping_;
346  }
347 
348  public:
349  MyALU2dGridGeometryImpl() : mapping_(), valid_( false ) {}
350 
351  // returns true if goemetry info is valid
352  bool valid () const { return valid_; }
353 
354  // reset geometry status
355  void invalidate() { valid_ = false ; }
356 
357  bool affine () const
358  {
359  return mapping().affine();
360  }
361 
362  int corners () const
363  {
364  return ncorners;
365  }
366 
367  GeometryType type () const
368  {
369  return GeometryType( GeometryType::cube, 2 );
370  }
371 
372  void map2world ( const map_t &m, world_t &w ) const
373  {
374  return mapping().map2world( m, w );
375  }
376 
377  void world2map ( const world_t &w, map_t &m ) const
378  {
379  return mapping().world2map( w, m );
380  }
381 
382  const matrix_t &jacobianTransposed ( const map_t &m ) const
383  {
384  return mapping().jacobianTransposed( m );
385  }
386 
387  const inv_t &jacobianInverseTransposed ( const map_t &m ) const
388  {
389  return mapping().jacobianInverseTransposed( m );
390  }
391 
392  ctype det ( const map_t &m ) const
393  {
394  return mapping().det( m );
395  }
396 
397  // update geometry in father coordinates
398  template< class Geo, class LocalGeo >
399  void updateLocal ( const Geo &geo, const LocalGeo &localGeo )
400  {
401  assert( localGeo.corners() == ncorners );
402  // compute the local coordinates in father refelem
403  FieldMatrix< alu2d_ctype, ncorners, cdim > coord;
404  for( int i = 0; i < ncorners; ++i )
405  {
406  // calculate coordinate
407  coord[ i ] = geo.local( localGeo.corner( i ) );
408  // to avoid rounding errors
409  for( int j = 0; j < cdim; ++j )
410  coord[ i ][ j ] = (coord[ i ][ j ] < 1e-14 ? 0 : coord[ i ][ j ]);
411  }
412  mapping_.buildMapping( coord[ 0 ], coord[ 1 ], coord[ 2 ], coord[ 3 ] );
413  valid_ = true ;
414  }
415 
416  template< class HElement >
417  void update ( const HElement &item )
418  {
419  mapping_.buildMapping( item.getVertex( 0 )->coord(), item.getVertex( 1 )->coord(),
420  item.getVertex( 3 )->coord(), item.getVertex( 2 )->coord() );
421  valid_ = true ;
422  }
423  };
424 
425  // geometry implementation for triangles
426  template< int cdim >
427  class MyALU2dGridGeometryImpl< 2, cdim, ALU2DSPACE mixed >
428  {
431 
432  typedef typename LinearMapping::ctype ctype;
433 
434  typedef typename LinearMapping::map_t map_t;
435  typedef typename LinearMapping::world_t world_t;
436 
437  typedef typename LinearMapping::matrix_t matrix_t;
438  typedef typename LinearMapping::inv_t inv_t;
439 
440  static const int lms = sizeof( LinearMapping );
441  static const int bms = sizeof( BilinearMapping );
442 
443  char mapping_[ lms > bms ? lms : bms ];
444  int corners_;
445  bool valid_ ;
446 
447  public:
448  MyALU2dGridGeometryImpl () : corners_( 0 ), valid_( false ) {}
449 
450  MyALU2dGridGeometryImpl ( const MyALU2dGridGeometryImpl &other )
451  : corners_( other.corners() ), valid_( other.valid_ )
452  {
453  if( corners_ == 3 )
454  new( &mapping_ ) LinearMapping( other.linearMapping() );
455  if( corners_ == 4 )
456  new( &mapping_ ) BilinearMapping( other.bilinearMapping() );
457  }
458 
459  // returns true if goemetry info is valid
460  bool valid () const { return valid_; }
461 
462  // reset geometry status
463  void invalidate() { valid_ = false ; }
464 
465  bool affine () const
466  {
467  return (corners() == 3 ? linearMapping().affine() : bilinearMapping().affine());
468  }
469 
470  int corners () const
471  {
472  return corners_;
473  }
474 
475  GeometryType type () const
476  {
477  return GeometryType( (corners_ == 3 ?
478  GenericGeometry :: SimplexTopology< 2 > :: type :: id :
479  GenericGeometry :: CubeTopology < 2 > :: type :: id), 2);
480  }
481 
482  void map2world ( const map_t &m, world_t &w ) const
483  {
484  if( corners() == 3 )
485  linearMapping().map2world( m, w );
486  else
487  bilinearMapping().map2world( m, w );
488  }
489 
490  void world2map ( const world_t &w, map_t &m ) const
491  {
492  if( corners() == 3 )
493  linearMapping().world2map( w, m );
494  else
495  bilinearMapping().world2map( w, m );
496  }
497 
498  const matrix_t &jacobianTransposed ( const map_t &m ) const
499  {
500  return (corners() == 3 ? linearMapping().jacobianTransposed( m ) : bilinearMapping().jacobianTransposed( m ));
501  }
502 
503  const inv_t &jacobianInverseTransposed ( const map_t &m ) const
504  {
505  return (corners() == 3 ? linearMapping().jacobianInverseTransposed( m ) : bilinearMapping().jacobianInverseTransposed( m ));
506  }
507 
508  ctype det ( const map_t &m ) const
509  {
510  return (corners() == 3 ? linearMapping().det( m ) : bilinearMapping().det( m ));
511  }
512 
513  // update geometry in father coordinates
514  template< class Geo, class LocalGeo >
515  void updateLocal ( const Geo &geo, const LocalGeo &localGeo )
516  {
517  const int corners = localGeo.corners();
518 
519  // compute the local coordinates in father refelem
520  FieldMatrix< alu2d_ctype, 4, cdim > coord;
521  for( int i = 0; i < corners; ++i )
522  {
523  // calculate coordinate
524  coord[ i ] = geo.local( localGeo.corner( i ) );
525  // to avoid rounding errors
526  for( int j = 0; j < cdim; ++j )
527  coord[ i ][ j ] = (coord[ i ][ j ] < 1e-14 ? 0 : coord[ i ][ j ]);
528  }
529 
530  updateMapping( corners );
531  if( corners == 3 )
532  linearMapping().buildMapping( coord[ 0 ], coord[ 1 ], coord[ 2 ] );
533  else
534  bilinearMapping().buildMapping( coord[ 0 ], coord[ 1 ], coord[ 2 ], coord[ 3 ] );
535 
536  valid_ = true ;
537  }
538 
539  template< class HElement >
540  void update ( const HElement &item )
541  {
542  const int corners = item.numvertices();
543  updateMapping( corners );
544  if( corners == 3 )
545  linearMapping().buildMapping( item.getVertex( 0 )->coord(), item.getVertex( 1 )->coord(),
546  item.getVertex( 2 )->coord() );
547  else
548  bilinearMapping().buildMapping( item.getVertex( 0 )->coord(), item.getVertex( 1 )->coord(),
549  item.getVertex( 3 )->coord(), item.getVertex( 2 )->coord() );
550 
551  valid_ = true ;
552  }
553 
554  private:
555  MyALU2dGridGeometryImpl &operator= ( const MyALU2dGridGeometryImpl &other );
556 
557  const LinearMapping &linearMapping () const
558  {
559  assert( valid() );
560  return static_cast< const LinearMapping * >( &mapping_ );
561  }
562 
563  LinearMapping &linearMapping ()
564  {
565  assert( valid() );
566  return static_cast< LinearMapping * >( &mapping_ );
567  }
568 
569  const BilinearMapping &bilinearMapping () const
570  {
571  assert( valid() );
572  return static_cast< const BilinearMapping * >( &mapping_ );
573  }
574 
575  BilinearMapping &bilinearMapping ()
576  {
577  assert( valid() );
578  return static_cast< BilinearMapping * >( &mapping_ );
579  }
580 
581  void updateMapping ( const int corners )
582  {
583  assert( (corners == 3) || (corners == 4) );
584  if( corners != corners_ )
585  {
586  destroyMapping();
587  corners = corners_;
588  if( corners == 3 )
589  new( &mapping_ ) LinearMapping;
590  else
591  new( &mapping_ ) BilinearMapping;
592  }
593  }
594 
595  void destroyMapping ()
596  {
597  if( corners() == 3 )
598  linearMapping().~LinearMapping();
599  else if( corners() == 4 )
600  bilinearMapping().~BilinearMapping();
601  }
602  };
603 
604 
605  //**********************************************************************
606  //
607  // --ALU2dGridGeometry
608  // --Geometry
609  //**********************************************************************
622 
623 
624  template< int mydim, int cdim, class GridImp >
626  : public GeometryDefaultImplementation< mydim, cdim, GridImp, ALU2dGridGeometry >
627  {
628  static const ALU2DSPACE ElementType eltype = GridImp::elementType;
629 
631  typedef typename GridImp::template Codim<0>::Geometry Geometry;
635  enum { dimbary=mydim+1};
636 
638  typedef typename ALU2dImplInterface< 0, GridImp::dimensionworld, eltype >::Type VertexType;
639 
640  // type of specialized geometry implementation
641  typedef MyALU2dGridGeometryImpl< mydim, cdim, eltype > GeometryImplType;
642 
643  public:
644  typedef FieldVector< alu2d_ctype, cdim > GlobalCoordinate;
645  typedef FieldVector< alu2d_ctype, mydim > LocalCoordinate;
646 
647  public:
651 
654  const GeometryType type () const { return geoImpl_.type(); }
655 
657  int corners () const { return geoImpl_.corners(); }
658 
660  const GlobalCoordinate &operator[] ( int i ) const;
661 
663  GlobalCoordinate corner ( int i ) const;
664 
667  GlobalCoordinate global ( const LocalCoordinate& local ) const;
668 
672 
675 
677  alu2d_ctype volume () const;
678 
680  bool affine() const { return geoImpl_.affine(); }
681 
683  const FieldMatrix<alu2d_ctype,cdim,mydim>& jacobianInverseTransposed (const LocalCoordinate& local) const;
684 
686  const FieldMatrix<alu2d_ctype,mydim,cdim>& jacobianTransposed (const LocalCoordinate& local) const;
687 
688  //***********************************************************************
690  //***********************************************************************
692  // method for elements
693  bool buildGeom(const HElementType & item);
694  // method for edges
695  bool buildGeom(const HElementType & item, const int aluFace);
696  // method for vertices
697  bool buildGeom(const VertexType & item, const int );
698 
701  template <class GeometryType, class LocalGeomType >
702  bool buildLocalGeom(const GeometryType & geo , const LocalGeomType & lg);
703 
705  bool buildLocalGeometry(const int faceNumber, const int twist,const int coorns);
706 
708  GlobalCoordinate& getCoordVec (int i);
709 
711  void print (std::ostream& ss) const;
712 
714  inline bool buildGeomInFather(const Geometry &fatherGeom,
715  const Geometry & myGeom );
716 
717  // returns true if geometry information is valid
718  inline bool valid() const { return geoImpl_.valid(); }
719 
720  // invalidate geometry information
721  inline void invalidate() const { geoImpl_.invalidate(); }
722 
723  protected:
724  // return reference coordinates of the alu triangle
725  static std::pair< FieldMatrix< alu2d_ctype, 4, 2 >, FieldVector< alu2d_ctype, 4 > >
726  calculateReferenceCoords ( const int corners );
727 
728  // implementation of coord and mapping
729  mutable GeometryImplType geoImpl_;
730 
731  // determinant
732  mutable alu2d_ctype det_;
733  };
734 
735 } // end namespace Dune
736 
737 #include "geometry_imp.cc"
738 
739 #endif