1 #ifndef DUNE_ALU2DGRIDGEOMETRY_HH
2 #define DUNE_ALU2DGRIDGEOMETRY_HH
5 #include <dune/common/misc.hh>
7 #include <dune/geometry/genericgeometry/topologytypes.hh>
16 template<
int cd,
int dim,
class Gr
idImp>
17 class ALU2dGridEntity;
18 template<
int cd,
class Gr
idImp >
19 class ALU2dGridEntityPointer;
20 template<
int mydim,
int cdim,
class Gr
idImp>
21 class ALU2dGridGeometry;
22 template<
int dim,
int dimworld, ALU2DSPACE ElementType eltype >
26 template<
int mydim,
int cdim, ALU2DSPACE ElementType eltype >
27 class MyALU2dGridGeometryImpl;
30 template<
int cdim, ALU2DSPACE ElementType eltype >
31 class MyALU2dGridGeometryImpl< 0, cdim, eltype >
56 bool valid ()
const {
return valid_; }
64 return mapping().affine();
76 GenericGeometry :: SimplexTopology< 0 > :: type :: id :
77 GenericGeometry :: CubeTopology < 0 > :: type :: id),
81 void map2world (
const map_t &m, world_t &w )
const
83 return mapping().map2world( m, w );
86 void world2map (
const world_t &w, map_t &m )
const
88 return mapping().world2map( w, m );
91 const matrix_t &jacobianTransposed (
const map_t &m )
const
93 return mapping().jacobianTransposed( m );
96 const inv_t &jacobianInverseTransposed (
const map_t &m )
const
98 return mapping().jacobianInverseTransposed( m );
101 ctype det (
const map_t &m )
const
103 return mapping().det( m );
107 template<
class Vector >
108 void update (
const Vector &p0 )
110 mapping_.buildMapping( p0 );
116 template<
int cdim, ALU2DSPACE ElementType eltype >
117 class MyALU2dGridGeometryImpl< 1, cdim, eltype >
119 static const int ncorners = 2;
144 bool valid ()
const {
return valid_; }
151 return mapping().affine();
163 GenericGeometry :: SimplexTopology< 1 > :: type :: id :
164 GenericGeometry :: CubeTopology < 1 > :: type :: id),
168 void map2world (
const map_t &m, world_t &w )
const
170 return mapping().map2world( m, w );
173 void world2map (
const world_t &w, map_t &m )
const
175 return mapping().world2map( w, m );
178 const matrix_t &jacobianTransposed (
const map_t &m )
const
180 return mapping().jacobianTransposed( m );
183 const inv_t &jacobianInverseTransposed (
const map_t &m )
const
185 return mapping().jacobianInverseTransposed( m );
188 ctype det (
const map_t &m )
const
190 return mapping().det( m );
194 template<
class Geo,
class LocalGeo >
195 void updateLocal (
const Geo &geo,
const LocalGeo &localGeo )
197 assert( localGeo.corners() == ncorners );
199 FieldMatrix< alu2d_ctype, ncorners, cdim > coord;
200 for(
int i = 0; i < ncorners; ++i )
203 coord[ i ] = geo.local( localGeo.corner( i ) );
205 for(
int j = 0; j < cdim; ++j )
206 coord[ i ][ j ] = (coord[ i ][ j ] < 1e-14 ? 0 : coord[ i ][ j ]);
208 mapping_.buildMapping( coord[ 0 ], coord[ 1 ] );
213 template<
class Vector >
214 void update (
const Vector &p0,
const Vector &p1 )
216 mapping_.buildMapping( p0, p1 );
225 static const int ncorners = 3;
250 bool valid ()
const {
return valid_; }
257 return mapping().affine();
267 return GeometryType( GenericGeometry :: SimplexTopology< 2 > :: type :: id , 2 );
270 void map2world (
const map_t &m, world_t &w )
const
272 return mapping().map2world( m, w );
275 void world2map (
const world_t &w, map_t &m )
const
277 return mapping().world2map( w, m );
280 const matrix_t &jacobianTransposed (
const map_t &m )
const
282 return mapping().jacobianTransposed( m );
285 const inv_t &jacobianInverseTransposed (
const map_t &m )
const
287 return mapping().jacobianInverseTransposed( m );
290 ctype det (
const map_t &m )
const
292 return mapping().det( m );
296 template<
class Geo,
class LocalGeo >
297 void updateLocal (
const Geo &geo,
const LocalGeo &localGeo )
299 assert( localGeo.corners() == ncorners );
301 FieldMatrix< alu2d_ctype, ncorners, cdim > coord;
302 for(
int i = 0; i < ncorners; ++i )
305 coord[ i ] = geo.local( localGeo.corner( i ) );
307 for(
int j = 0; j < cdim; ++j )
308 coord[ i ][ j ] = (coord[ i ][ j ] < 1e-14 ? 0 : coord[ i ][ j ]);
310 mapping_.buildMapping( coord[ 0 ], coord[ 1 ], coord[ 2 ] );
314 template<
class HElement >
315 void update (
const HElement &item )
317 mapping_.buildMapping( item.getVertex( 0 )->coord(), item.getVertex( 1 )->coord(),
318 item.getVertex( 2 )->coord() );
327 static const int ncorners = 4;
352 bool valid ()
const {
return valid_; }
359 return mapping().affine();
372 void map2world (
const map_t &m, world_t &w )
const
374 return mapping().map2world( m, w );
377 void world2map (
const world_t &w, map_t &m )
const
379 return mapping().world2map( w, m );
382 const matrix_t &jacobianTransposed (
const map_t &m )
const
384 return mapping().jacobianTransposed( m );
387 const inv_t &jacobianInverseTransposed (
const map_t &m )
const
389 return mapping().jacobianInverseTransposed( m );
392 ctype det (
const map_t &m )
const
394 return mapping().det( m );
398 template<
class Geo,
class LocalGeo >
399 void updateLocal (
const Geo &geo,
const LocalGeo &localGeo )
401 assert( localGeo.corners() == ncorners );
403 FieldMatrix< alu2d_ctype, ncorners, cdim > coord;
404 for(
int i = 0; i < ncorners; ++i )
407 coord[ i ] = geo.local( localGeo.corner( i ) );
409 for(
int j = 0; j < cdim; ++j )
410 coord[ i ][ j ] = (coord[ i ][ j ] < 1e-14 ? 0 : coord[ i ][ j ]);
412 mapping_.buildMapping( coord[ 0 ], coord[ 1 ], coord[ 2 ], coord[ 3 ] );
416 template<
class HElement >
417 void update (
const HElement &item )
419 mapping_.buildMapping( item.getVertex( 0 )->coord(), item.getVertex( 1 )->coord(),
420 item.getVertex( 3 )->coord(), item.getVertex( 2 )->coord() );
443 char mapping_[ lms > bms ? lms : bms ];
450 MyALU2dGridGeometryImpl (
const MyALU2dGridGeometryImpl &other )
451 : corners_( other.corners() ), valid_( other.valid_ )
460 bool valid ()
const {
return valid_; }
467 return (corners() == 3 ? linearMapping().affine() : bilinearMapping().affine());
478 GenericGeometry :: SimplexTopology< 2 > :: type :: id :
479 GenericGeometry :: CubeTopology < 2 > :: type :: id), 2);
482 void map2world (
const map_t &m, world_t &w )
const
485 linearMapping().map2world( m, w );
487 bilinearMapping().map2world( m, w );
490 void world2map (
const world_t &w, map_t &m )
const
493 linearMapping().world2map( w, m );
495 bilinearMapping().world2map( w, m );
498 const matrix_t &jacobianTransposed (
const map_t &m )
const
500 return (corners() == 3 ? linearMapping().jacobianTransposed( m ) : bilinearMapping().jacobianTransposed( m ));
503 const inv_t &jacobianInverseTransposed (
const map_t &m )
const
505 return (corners() == 3 ? linearMapping().jacobianInverseTransposed( m ) : bilinearMapping().jacobianInverseTransposed( m ));
508 ctype det (
const map_t &m )
const
510 return (corners() == 3 ? linearMapping().det( m ) : bilinearMapping().det( m ));
514 template<
class Geo,
class LocalGeo >
515 void updateLocal (
const Geo &geo,
const LocalGeo &localGeo )
517 const int corners = localGeo.corners();
520 FieldMatrix< alu2d_ctype, 4, cdim > coord;
521 for(
int i = 0; i < corners; ++i )
524 coord[ i ] = geo.local( localGeo.corner( i ) );
526 for(
int j = 0; j < cdim; ++j )
527 coord[ i ][ j ] = (coord[ i ][ j ] < 1e-14 ? 0 : coord[ i ][ j ]);
530 updateMapping( corners );
532 linearMapping().buildMapping( coord[ 0 ], coord[ 1 ], coord[ 2 ] );
534 bilinearMapping().buildMapping( coord[ 0 ], coord[ 1 ], coord[ 2 ], coord[ 3 ] );
539 template<
class HElement >
540 void update (
const HElement &item )
542 const int corners = item.numvertices();
543 updateMapping( corners );
545 linearMapping().buildMapping( item.getVertex( 0 )->coord(), item.getVertex( 1 )->coord(),
546 item.getVertex( 2 )->coord() );
548 bilinearMapping().buildMapping( item.getVertex( 0 )->coord(), item.getVertex( 1 )->coord(),
549 item.getVertex( 3 )->coord(), item.getVertex( 2 )->coord() );
555 MyALU2dGridGeometryImpl &operator= (
const MyALU2dGridGeometryImpl &other );
563 LinearMapping &linearMapping ()
566 return static_cast< LinearMapping *
>( &mapping_ );
569 const BilinearMapping &bilinearMapping ()
const
572 return static_cast< const BilinearMapping *
>( &mapping_ );
575 BilinearMapping &bilinearMapping ()
578 return static_cast< BilinearMapping *
>( &mapping_ );
581 void updateMapping (
const int corners )
583 assert( (corners == 3) || (corners == 4) );
584 if( corners != corners_ )
589 new( &mapping_ ) LinearMapping;
591 new( &mapping_ ) BilinearMapping;
595 void destroyMapping ()
598 linearMapping().~LinearMapping();
599 else if( corners() == 4 )
600 bilinearMapping().~BilinearMapping();
624 template<
int mydim,
int cdim,
class Gr
idImp >
631 typedef typename GridImp::template Codim<0>::Geometry Geometry;
635 enum { dimbary=mydim+1};
638 typedef typename ALU2dImplInterface< 0, GridImp::dimensionworld, eltype >::Type VertexType;
641 typedef MyALU2dGridGeometryImpl< mydim, cdim, eltype > GeometryImplType;
693 bool buildGeom(
const HElementType & item);
695 bool buildGeom(
const HElementType & item,
const int aluFace);
697 bool buildGeom(
const VertexType & item,
const int );
701 template <
class GeometryType,
class LocalGeomType >
711 void print (std::ostream& ss)
const;
725 static std::pair< FieldMatrix< alu2d_ctype, 4, 2 >, FieldVector< alu2d_ctype, 4 > >
737 #include "geometry_imp.cc"