dune-grid  2.2.0
mappings.hh
Go to the documentation of this file.
00001 #ifndef DUNE_ALU3DGRIDMAPPINGS_HH
00002 #define DUNE_ALU3DGRIDMAPPINGS_HH
00003 
00004 // System includes
00005 #include <limits>
00006 #include <cmath>
00007 
00008 // Dune includes
00009 #include <dune/common/fvector.hh>
00010 #include <dune/common/fmatrix.hh>
00011 #include <dune/common/exceptions.hh>
00012 
00013 // Local includes
00014 #include "alu3dinclude.hh"
00015 
00016 namespace Dune {
00017 
00018   static const alu3d_ctype ALUnumericEpsilon = 10.0 * std::numeric_limits< alu3d_ctype >::epsilon();
00019 
00020   template<int mydim, int coorddim, class GridImp>
00021   class ALU3dGridGeometry;
00022 
00023   template< ALU3dGridElementType, class >
00024   class ALU3dGrid;
00025 
00028   class TrilinearMapping 
00029   {
00030   public:  
00031     typedef alu3d_ctype double_t[3];
00032     typedef FieldVector<alu3d_ctype, 3> coord_t;
00033     typedef FieldMatrix<alu3d_ctype, 3, 3> mat_t;
00034   private:  
00035     static const double _epsilon ;
00036 
00037     // the internal mapping 
00038     alu3d_ctype a [8][3] ;
00039     mat_t Df; 
00040     mat_t Dfi; 
00041     mat_t invTransposed_; 
00042     alu3d_ctype DetDf ;
00043 
00044     bool calcedDet_; 
00045     bool calcedLinear_;
00046     bool calcedInv_; 
00047     bool affine_;
00048 
00049     void linear (const alu3d_ctype, const alu3d_ctype, const alu3d_ctype) ;
00050     void linear (const coord_t&) ;
00051     void inverse (const coord_t&) ;
00052   public :
00053     TrilinearMapping (const coord_t&, const coord_t&,
00054                       const coord_t&, const coord_t&,
00055                       const coord_t&, const coord_t&,
00056                       const coord_t&, const coord_t&);
00057     
00058     // only to call from geometry class 
00059     TrilinearMapping () {}
00060     
00061     TrilinearMapping (const TrilinearMapping &) ;
00062     
00063     ~TrilinearMapping () {}
00064     alu3d_ctype det (const coord_t&) ;
00065     const mat_t& jacobianInverseTransposed(const coord_t&);
00066     const mat_t& jacobianTransposed(const coord_t&);
00067     void map2world (const coord_t&, coord_t&) const ;
00068     void map2world (const alu3d_ctype , const alu3d_ctype , const alu3d_ctype , 
00069                     coord_t&) const ;
00070     void world2map (const coord_t&, coord_t&) ;
00071 
00072     template <class vector_t>
00073     void buildMapping(const vector_t&, const vector_t&,
00074                       const vector_t&, const vector_t&,  
00075                       const vector_t&, const vector_t&,  
00076                       const vector_t&, const vector_t&); 
00077 
00078     // returns true if mapping is affine 
00079     inline bool affine () const { return affine_; }
00080   };
00081 
00083   // NOTE: this class is different to the BilinearSurfaceMapping in
00084   // ALUGrid, for example the reference elements differ 
00085   // here we have [0,1]^2 and in ALUGrid its [-1,1]^2
00086   // also the point numbering is different 
00087   class SurfaceNormalCalculator 
00088   {
00089   public:  
00090     // our coordinate types  
00091     typedef FieldVector<alu3d_ctype, 3> coord3_t;
00092     typedef FieldVector<alu3d_ctype, 2> coord2_t;
00093 
00094     // type of coordinate vectors from elements 
00095     typedef alu3d_ctype double3_t[3];
00096   protected:  
00097 
00098     alu3d_ctype _n [3][3] ;
00099 
00100     static const double _epsilon ;
00101 
00102     bool _affine;
00103 
00104   public :
00106     SurfaceNormalCalculator();
00107 
00108     SurfaceNormalCalculator (const SurfaceNormalCalculator &) ;
00109     ~SurfaceNormalCalculator () {}
00110 
00111     // returns true if mapping is affine 
00112     inline bool affine () const { return _affine ; }
00113 
00114     // calcuates normal 
00115     void normal(const coord2_t&, coord3_t&) const ;
00116     void normal(const alu3d_ctype, const alu3d_ctype, coord3_t&)const;
00117     
00118     void negativeNormal(const coord2_t&, coord3_t&) const ;
00119     void negativeNormal(const alu3d_ctype, const alu3d_ctype, coord3_t&)const;
00120     
00121   public:
00122     // builds _b and _n, called from the constructors 
00123     // public because also used in faceutility 
00124     template <class vector_t>
00125     void buildMapping (const vector_t & , const vector_t & ,
00126                        const vector_t & , const vector_t & );
00127   protected:  
00128     // builds _b and _n, called from the constructors 
00129     // public because also used in faceutility 
00130     template <class vector_t>
00131     void buildMapping (const vector_t & , const vector_t & ,
00132                        const vector_t & , const vector_t & , 
00133                        alu3d_ctype (&_b)[4][3] );
00134   } ;  
00135 
00136 
00138   // NOTE: this class is different to the BilinearSurfaceMapping in
00139   // ALUGrid, for example the reference elements differ 
00140   // here we have [0,1]^2 and in ALUGrid its [-1,1]^2
00141   // also the point numbering is different 
00142   class BilinearSurfaceMapping : public SurfaceNormalCalculator 
00143   {
00144   protected:  
00145     typedef SurfaceNormalCalculator BaseType; 
00146     
00147     using BaseType :: _n;
00148     static const double _epsilon;
00149 
00150     // our coordinate types  
00151     typedef FieldVector<alu3d_ctype, 3> coord3_t;
00152     typedef FieldVector<alu3d_ctype, 2> coord2_t;
00153 
00154     // type of coordinate vectors from elements 
00155     typedef alu3d_ctype double3_t[3];
00156 
00157     // type for helper matrices 
00158     typedef FieldMatrix<alu3d_ctype,3,3> mat3_t;
00159 
00160     // type for inverse matrices 
00161     typedef FieldMatrix<alu3d_ctype,2,3> matrix_t;
00162 
00163     // type for inverse matrices 
00164     typedef FieldMatrix<alu3d_ctype,3,2> inv_t;
00165 
00166     alu3d_ctype _b [4][3] ;
00167     
00168     mutable mat3_t Df,Dfi;
00169     mutable inv_t invTransposed_;
00170     mutable matrix_t matrix_;
00171     mutable alu3d_ctype DetDf;
00172 
00173     mutable coord3_t normal_;
00174     mutable coord3_t tmp_;
00175 
00176     mutable bool _calcedInv;
00177     mutable bool _calcedTransposed;
00178     mutable bool _calcedMatrix;
00179 
00180   public :
00182     BilinearSurfaceMapping ();
00183 
00185     BilinearSurfaceMapping (const coord3_t&, const coord3_t&,
00186                             const coord3_t&, const coord3_t&) ;
00188     BilinearSurfaceMapping (const double3_t &, const double3_t &,
00189                             const double3_t &, const double3_t &) ;
00190     BilinearSurfaceMapping (const BilinearSurfaceMapping &) ;
00191     ~BilinearSurfaceMapping () {}
00192 
00193     void inverse (const coord3_t&) const;
00194     const inv_t& jacobianInverseTransposed(const coord2_t&) const ;
00195 
00196     const matrix_t& jacobianTransposed(const coord2_t&) const ;
00197 
00198     // calculates determinant of face mapping using the normal 
00199     alu3d_ctype det(const coord2_t&) const;
00200     
00201     // maps from local coordinates to global coordinates 
00202     void world2map(const coord3_t &, coord2_t & ) const;
00203 
00204     // maps form global coordinates to local (within reference element)
00205     // coordinates 
00206     void map2world(const coord2_t&, coord3_t&) const ;
00207     void map2world(const alu3d_ctype ,const alu3d_ctype , coord3_t&) const ;
00208     
00209   private:  
00210     void map2worldnormal(const alu3d_ctype, const alu3d_ctype, const alu3d_ctype , coord3_t&)const;
00211     void map2worldlinear(const alu3d_ctype, const alu3d_ctype, const alu3d_ctype ) const;
00212 
00213   public:
00214     // builds _b and _n, called from the constructors 
00215     // public because also used in faceutility 
00216     template <class vector_t>
00217     void buildMapping (const vector_t & , const vector_t & ,
00218                        const vector_t & , const vector_t & );
00219   } ;  
00220 
00221 
00222 
00224   template< int cdim >
00225   class BilinearMapping
00226   {
00227   public:  
00228     typedef alu3d_ctype ctype;
00229 
00230     typedef FieldVector< ctype, cdim > world_t;
00231     typedef FieldVector< ctype, 2 > map_t;
00232 
00233     typedef FieldMatrix< ctype, 2, cdim > matrix_t;
00234     typedef FieldMatrix< ctype, cdim, 2 > inv_t;
00235 
00236   protected:
00237     ctype _b [4][cdim];
00238 
00239     mutable ctype det_;
00240     mutable matrix_t matrix_;
00241     mutable inv_t invTransposed_;
00242 
00243     mutable bool affine_;
00244     mutable bool calcedMatrix_;
00245     mutable bool calcedDet_;
00246     mutable bool calcedInv_;
00247 
00248   public:  
00249     BilinearMapping ();
00250     BilinearMapping ( const world_t &p0, const world_t &p1,
00251                       const world_t &p2, const world_t &p3 );
00252     BilinearMapping ( const ctype (&p0)[ cdim ], const ctype (&p1)[ cdim ],
00253                       const ctype (&p2)[ cdim ], const ctype (&p3)[ cdim ] );
00254 
00255     bool affine () const;
00256 
00257     void world2map ( const world_t &, map_t & ) const;
00258     void map2world ( const ctype x, const ctype y, world_t &w ) const;
00259     void map2world ( const map_t &, world_t & ) const;
00260 
00261     ctype det ( const map_t & ) const;
00262 
00263     const matrix_t &jacobianTransposed ( const map_t & ) const;
00264     const inv_t &jacobianInverseTransposed ( const map_t & ) const;
00265 
00266     template< class vector_t >
00267     void buildMapping ( const vector_t &, const vector_t &,
00268                         const vector_t &, const vector_t & );
00269 
00270   protected:
00271     static void multTransposedMatrix ( const matrix_t &, FieldMatrix< ctype, 2, 2 > & );
00272     static void multMatrix ( const matrix_t &, const FieldMatrix< ctype, 2, 2 > &, inv_t & );
00273 
00274     void map2worldlinear ( const ctype, const ctype ) const;
00275     void inverse ( const map_t & ) const;
00276   };
00277 
00278 
00279 
00281   template< int cdim, int mydim >
00282   class LinearMapping 
00283   {
00284   public:  
00285     typedef alu3d_ctype ctype;
00286 
00287     typedef ctype double_t[ cdim ]; 
00288 
00289     typedef FieldVector< ctype, cdim > world_t;
00290     typedef FieldVector< ctype, mydim > map_t;
00291 
00292     typedef FieldMatrix< ctype, mydim, cdim > matrix_t;
00293     typedef FieldMatrix< ctype, cdim, mydim > inv_t;
00294 
00295   protected:  
00296     matrix_t      _matrix;        
00297     mutable inv_t _invTransposed; 
00298     world_t _p0;                  
00299 
00301     mutable ctype _det;                            
00302 
00304     mutable bool _calcedInv; 
00305 
00307     mutable bool _calcedDet; 
00308 
00309   public:  
00311     LinearMapping ();
00312 
00314     LinearMapping (const LinearMapping &) ;
00315     
00316     // returns true if mapping is affine (which is always true)
00317     inline bool affine () const { return true ; }
00318 
00319     // return reference to transposed jacobian  
00320     const matrix_t& jacobianTransposed(const map_t &) const ;
00321 
00322     // return reference to transposed jacobian inverse 
00323     const inv_t& jacobianInverseTransposed(const map_t &) const ;
00324 
00325     // calculates determinant of mapping 
00326     ctype det(const map_t&) const;
00327     
00328     // maps from local coordinates to global coordinates 
00329     void world2map(const world_t &, map_t &) const;
00330 
00331     // maps form global coordinates to local (within reference element)
00332     // coordinates 
00333     void map2world(const map_t &, world_t &) const ;
00334 
00335   protected:  
00336     // calculate inverse 
00337     void inverse (const map_t&) const;
00338     
00339     // calculate inverse one codim one entity 
00340     void inverseCodimOne (const map_t&) const;
00341     
00342     // calculate determinant  
00343     void calculateDeterminant (const map_t&) const;
00344 
00345     void multTransposedMatrix(const matrix_t& matrix,
00346                               FieldMatrix<ctype, mydim, mydim>& result) const;
00347 
00348     void multMatrix ( const matrix_t& A,
00349                       const FieldMatrix< ctype, mydim, mydim> &B,
00350                       inv_t& ret ) const ;
00351     
00352   public:
00353     // builds _b and _n, called from the constructors 
00354     // public because also used in faceutility 
00355     template <class vector_t>
00356     void buildMapping (const vector_t & , const vector_t & ,
00357                        const vector_t & , const vector_t & );
00358 
00359     // builds _b and _n, called from the constructors 
00360     // public because also used in faceutility 
00361     template <class vector_t>
00362     void buildMapping (const vector_t & , const vector_t & ,
00363                        const vector_t & );
00364 
00365     // builds _b and _n, called from the constructors 
00366     // public because also used in faceutility 
00367     template <class vector_t>
00368     void buildMapping (const vector_t & , const vector_t & );
00369 
00370     template <class vector_t>
00371     void buildMapping (const vector_t & );
00372   } ;  
00373 
00374 
00376   //
00377   // NonConforming Mappings 
00378   //
00380 
00381 
00384   template< ALU3dGridElementType type, class Comm >
00385   class NonConformingFaceMapping;
00386 
00388   template< class Comm >
00389   struct NonConformingFaceMapping< tetra, Comm > 
00390   {
00391     typedef FieldVector< alu3d_ctype, 3 > CoordinateType;
00392     typedef typename ALU3dImplTraits< tetra, Comm >::HfaceRuleType RefinementRuleType;
00393     
00394     NonConformingFaceMapping ( RefinementRuleType rule, int nChild )
00395     : rule_( rule ), nChild_( nChild )
00396     {}
00397     
00398     void child2parent ( const CoordinateType &childCoordinates,
00399                         CoordinateType &parentCoordinates) const;
00400 
00401     CoordinateType child2parent ( const FieldVector< alu3d_ctype, 2 > &childCoordinates ) const;
00402 
00403   private:
00404     void child2parentNosplit(const CoordinateType& childCoordinates,
00405                              CoordinateType& parentCoordinates) const;
00406     void child2parentE01(const CoordinateType& childCoordinates,
00407                          CoordinateType& parentCoordinates) const;
00408     void child2parentE12(const CoordinateType& childCoordinates,
00409                          CoordinateType& parentCoordinates) const;
00410     void child2parentE20(const CoordinateType& childCoordinates,
00411                          CoordinateType& parentCoordinates) const;
00412     void child2parentIso4(const CoordinateType& childCoordinates,
00413                           CoordinateType& parentCoordinates) const;
00414 
00415     RefinementRuleType rule_;
00416     int nChild_;
00417   };
00418 
00420   template< class Comm >
00421   struct NonConformingFaceMapping< hexa, Comm >
00422   {
00423     typedef FieldVector< alu3d_ctype, 2 > CoordinateType;
00424     typedef typename ALU3dImplTraits< hexa, Comm >::HfaceRuleType RefinementRuleType;
00425 
00426     NonConformingFaceMapping ( RefinementRuleType rule, int nChild )
00427     : rule_( rule ), nChild_( nChild )
00428     {}
00429 
00430     void child2parent ( const CoordinateType &childCoordinates,
00431                         CoordinateType &parentCoordinates) const;
00432 
00433     CoordinateType child2parent ( const FieldVector< alu3d_ctype, 2 > &childCoordinates ) const;
00434 
00435   private:
00436     void child2parentNosplit(const CoordinateType& childCoordinates,
00437                              CoordinateType& parentCoordinates) const;
00438     void child2parentIso4(const CoordinateType& childCoordinates,
00439                           CoordinateType& parentCoordinates) const;
00440 
00441     RefinementRuleType rule_;
00442     int nChild_; 
00443   };
00444 
00445 } // end namespace Dune
00446 
00447 #if COMPILE_ALUGRID_INLINE
00448   #include "mappings_imp.cc"
00449 #endif
00450 #endif