MOAB: Mesh Oriented datABase  (version 5.2.1)
linear_hex_map.hpp
Go to the documentation of this file.
00001 #ifndef MOAB_LINEAR_HEX_HPP
00002 #define MOAB_LINEAR_HEX_HPP
00003 
00004 #include "moab/Matrix3.hpp"
00005 #include "moab/CartVect.hpp"
00006 #include <sstream>
00007 #include <iomanip>
00008 #include <iostream>
00009 
00010 namespace moab
00011 {
00012 
00013 namespace element_utility
00014 {
00015 
00016     class Linear_hex_map
00017     {
00018       public:
00019         // Constructor
00020         Linear_hex_map() {}
00021         // Copy constructor
00022         Linear_hex_map( const Linear_hex_map& ) {}
00023 
00024       public:
00025         // Natural coordinates
00026         template < typename Moab, typename Entity_handle, typename Points, typename Point >
00027         std::pair< bool, CartVect > evaluate_reverse( const double* verts, const double* eval_point,
00028                                                       const double tol = 1.e-6 ) const
00029         {
00030             CartVect params( 3, 0.0 );
00031             solve_inverse( eval_point, params, verts );
00032             bool point_found = solve_inverse( eval_point, params, verts, tol ) && is_contained( params, tol );
00033             return std::make_pair( point_found, params );
00034         }
00035 
00036       private:
00037         // This is a hack to avoid a .cpp file and C++11
00038         // reference_points(i,j) will be a 1 or -1;
00039         // This should unroll..
00040         inline const CartVect& reference_points( const std::size_t& i ) const
00041         {
00042             const CartVect rpts[8] = { CartVect( -1, -1, -1 ), CartVect( 1, -1, -1 ), CartVect( 1, 1, -1 ),
00043                                        CartVect( -1, 1, -1 ),  CartVect( -1, -1, 1 ), CartVect( 1, -1, 1 ),
00044                                        CartVect( 1, 1, 1 ),    CartVect( -1, 1, 1 ) };
00045             return rpts[i];
00046         }
00047 
00048         bool is_contained( const CartVect& params, const double tol ) const
00049         {
00050             // just look at the box+tol here
00051             return ( params[0] >= -1. - tol ) && ( params[0] <= 1. + tol ) && ( params[1] >= -1. - tol ) &&
00052                    ( params[1] <= 1. + tol ) && ( params[2] >= -1. - tol ) && ( params[2] <= 1. + tol );
00053         }
00054 
00055         bool solve_inverse( const CartVect& point, CartVect& params, const CartVect* verts,
00056                             const double tol = 1.e-6 ) const
00057         {
00058             const double error_tol_sqr = tol * tol;
00059             CartVect delta( 0.0, 0.0, 0.0 );
00060             params = delta;
00061             evaluate_forward( params, verts, delta );
00062             delta -= point;
00063             std::size_t num_iterations = 0;
00064 #ifdef LINEAR_HEX_DEBUG
00065             std::stringstream ss;
00066             ss << "CartVect: ";
00067             ss << point[0] << ", " << point[1] << ", " << point[2] << std::endl;
00068             ss << "Hex: ";
00069             for( int i = 0; i < 8; ++i )
00070                 ss << points[i][0] << ", " << points[i][1] << ", " << points[i][2] << std::endl;
00071             ss << std::endl;
00072 #endif
00073             while( delta.length_squared() > error_tol_sqr )
00074             {
00075 #ifdef LINEAR_HEX_DEBUG
00076                 ss << "Iter #: " << num_iterations << " Err: " << delta.length() << " Iterate: ";
00077                 ss << params[0] << ", " << params[1] << ", " << params[2] << std::endl;
00078 #endif
00079                 if( ++num_iterations >= 5 ) { return false; }
00080                 Matrix3 J;
00081                 jacobian( params, verts, J );
00082                 double det = moab::Matrix3::determinant3( J );
00083                 if( fabs( det ) < 1.e-10 )
00084                 {
00085 #ifdef LINEAR_HEX_DEBUG
00086                     std::cerr << ss.str();
00087 #endif
00088 #ifndef LINEAR_HEX_DEBUG
00089                     std::cerr << x[0] << ", " << x[1] << ", " << x[2] << std::endl;
00090 #endif
00091                     std::cerr << "inverse solve failure: det: " << det << std::endl;
00092                     exit( -1 );
00093                 }
00094                 params -= moab::Matrix3::inverse( J, 1.0 / det ) * delta;
00095                 evaluate_forward( params, points, delta );
00096                 delta -= x;
00097             }
00098             return true;
00099         }
00100 
00101         void evaluate_forward( const CartVect& p, const CartVect* verts, CartVect& f ) const
00102         {
00103             typedef typename Points::value_type Vector;
00104             f.set( 0.0, 0.0, 0.0 );
00105             for( unsigned i = 0; i < 8; ++i )
00106             {
00107                 const double N_i = ( 1 + p[0] * reference_points( i )[0] ) * ( 1 + p[1] * reference_points( i )[1] ) *
00108                                    ( 1 + p[2] * reference_points( i )[2] );
00109                 f += N_i * verts[i];
00110             }
00111             f *= 0.125;
00112             return f;
00113         }
00114 
00115         double integrate_scalar_field( const CartVect* points, const double* field_values ) const {}
00116 
00117         template < typename Point, typename Points >
00118         Matrix3& jacobian( const Point& p, const Points& points, Matrix3& J ) const
00119         {
00120         }
00121 
00122       private:
00123     };  // Class Linear_hex_map
00124 
00125 }  // namespace element_utility
00126 
00127 }  // namespace moab
00128 #endif  // MOAB_LINEAR_HEX_nPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines