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