moab
linear_tet_map.hpp
Go to the documentation of this file.
00001 #ifndef MOAB_LINEAR_TET_HPP
00002 #define MOAB_LINEAR_TET_HPP
00003 
00004 #include "moab/Matrix3.hpp"
00005 
00006 namespace moab { 
00007 namespace element_utility {
00008 
00009 template<  typename Entity_handle, typename Matrix>
00010 class Linear_tet_map {
00011   private:
00012     typedef Linear_tet_map< Entity_handle, Matrix> Self;
00013   public: 
00014     //Constructor
00015     Linear_tet_map() : Tinv(), eh() {}
00016     //Copy constructor
00017     Linear_tet_map( const Self & f ) : Tinv( f.Tinv), eh( f.eh){}
00018     //Natural coordinates
00019     template< typename Moab, typename Points, typename Point>
00020     std::pair< bool, Point> operator()( const Moab & moab,
00021                     const Entity_handle _eh,
00022                         const Points & v, 
00023                     const Point & p, 
00024                     const double tol=1e-6) {
00025       // Remove the warning about unused parameter
00026       if (NULL != &moab) {}
00027 
00028       set_tet( _eh, v);
00029       //TODO: Make sure this is correct
00030       Point result = Tinv*p;
00031       return std::make_pair( is_contained( result, tol), result);
00032     }
00033   
00034     private:
00035     template< typename Point>
00036     bool is_contained( const Point & result, const double tol=1e-6){
00037     double sum=0.0;
00038     for( std::size_t i = 0; i < 3; ++i){ 
00039         sum += result[ i]; 
00040         if( result[ i] < -tol){ return false; } 
00041     }
00042     return sum < 1.0+tol;
00043     }
00044     template< typename Point, typename Field>
00045     double evaluate_scalar_field( const Point & p , 
00046                   const Field & field_values) const{
00047     double f0 = field_values[ 0];
00048     double f = f0;
00049     for(std::size_t i = 1; i < 5; ++i){f+=(field_values[ i] - f0)*p[ i -1];}
00050     return f;
00051     }
00052     template< typename Points, typename Field>
00053     double integrate_scalar_field(const Points & v, const Field & field_values) const {
00054       double I(0.0);
00055       for(unsigned int i = 0; i < 4; ++i) { I += field_values[i]; }
00056       double det = Matrix( v[1][0]-v[0][0], v[2][0]-v[0][0], 
00057                    v[3][0]-v[0][0],
00058                    v[1][1]-v[0][1], v[2][1]-v[0][1], 
00059                    v[3][1]-v[0][1],
00060                    v[1][2]-v[0][2], v[2][2]-v[0][2], 
00061                    v[3][2]-v[0][2]).determinant();
00062       I *= det/24.0; 
00063       return I;
00064     }
00065     
00066     template< typename Points>
00067     void set_tet( const Entity_handle _eh, const Points & v){
00068         if (eh != _eh){
00069            eh = _eh;
00070            Tinv = moab::Matrix::inverse( 
00071                Matrix( v[1][0]-v[0][0], v[2][0]-v[0][0], 
00072                    v[3][0]-v[0][0],
00073                                v[1][1]-v[0][1], v[2][1]-v[0][1], 
00074                    v[3][1]-v[0][1],
00075                                v[1][2]-v[0][2], v[2][2]-v[0][2], 
00076                    v[3][2]-v[0][2]) );
00077         }
00078     }
00079   private:
00080     Matrix Tinv;
00081     Entity_handle eh;
00082     /* We don't need this!
00083     static const double reference_points[ 4][ 3] = { {0,0,0},
00084                                  {1,0,0},
00085                                  {0,1,0},
00086                                  {0,0,1} };*/
00087     
00088 }; //Class Linear_tet_map
00089 
00090 }// namespace element_utility
00091 
00092 } // namespace moab
00093 #endif //MOAB_LINEAR_TET_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines