linear_hex_map.hpp

Go to the documentation of this file.

    1 #ifndef MOAB_LINEAR_HEX_HPP
    2 #define MOAB_LINEAR_HEX_HPP
    3  
    4 #include "moab/Matrix3.hpp"
    5 #include "moab/CartVect.hpp"
    6 #include <sstream>
    7 #include <iomanip>
    8 #include <iostream>
    9  
   10 namespace moab
   11 {
   12  
   13 namespace element_utility
   14 {
   15  
   16     class Linear_hex_map
   17     {
   18       public:
   19         // Constructor
   20         Linear_hex_map() {}
   21         // Copy constructor
   22         Linear_hex_map( const Linear_hex_map& ) {}
   23  
   24       public:
   25         // Natural coordinates
   26         template < typename Moab, typename Entity_handle, typename Points, typename Point >
   27         std::pair< bool, CartVect > evaluate_reverse( const double* verts,
   28                                                       const double* eval_point,
   29                                                       const double tol = 1.e-6 ) const
   30         {
   31             CartVect params( 3, 0.0 );
   32             solve_inverse( eval_point, params, verts );
   33             bool point_found = solve_inverse( eval_point, params, verts, tol ) && is_contained( params, tol );
   34             return std::make_pair( point_found, params );
   35         }
   36  
   37       private:
   38         // This is a hack to avoid a .cpp file and C++11
   39         // reference_points(i,j) will be a 1 or -1;
   40         // This should unroll..
   41         inline const CartVect& reference_points( const std::size_t& i ) const
   42         {
   43             const CartVect rpts[8] = { CartVect( -1, -1, -1 ), CartVect( 1, -1, -1 ), CartVect( 1, 1, -1 ),
   44                                        CartVect( -1, 1, -1 ),  CartVect( -1, -1, 1 ), CartVect( 1, -1, 1 ),
   45                                        CartVect( 1, 1, 1 ),    CartVect( -1, 1, 1 ) };
   46             return rpts[i];
   47         }
   48  
   49         bool is_contained( const CartVect& params, const double tol ) const
   50         {
   51             // just look at the box+tol here
   52             return ( params[0] >= -1. - tol ) && ( params[0] <= 1. + tol ) && ( params[1] >= -1. - tol ) &&
   53                    ( params[1] <= 1. + tol ) && ( params[2] >= -1. - tol ) && ( params[2] <= 1. + tol );
   54         }
   55  
   56         bool solve_inverse( const CartVect& point,
   57                             CartVect& params,
   58                             const CartVect* verts,
   59                             const double tol = 1.e-6 ) const
   60         {
   61             const double error_tol_sqr = tol * tol;
   62             CartVect delta( 0.0, 0.0, 0.0 );
   63             params = delta;
   64             evaluate_forward( params, verts, delta );
   65             delta -= point;
   66             std::size_t num_iterations = 0;
   67 #ifdef LINEAR_HEX_DEBUG
   68             std::stringstream ss;
   69             ss << "CartVect: ";
   70             ss << point[0] << ", " << point[1] << ", " << point[2] << std::endl;
   71             ss << "Hex: ";
   72             for( int i = 0; i < 8; ++i )
   73                 ss << points[i][0] << ", " << points[i][1] << ", " << points[i][2] << std::endl;
   74             ss << std::endl;
   75 #endif
   76             while( delta.length_squared() > error_tol_sqr )
   77             {
   78 #ifdef LINEAR_HEX_DEBUG
   79                 ss << "Iter #: " << num_iterations << " Err: " << delta.length() << " Iterate: ";
   80                 ss << params[0] << ", " << params[1] << ", " << params[2] << std::endl;
   81 #endif
   82                 if( ++num_iterations >= 5 )
   83                 {
   84                     return false;
   85                 }
   86                 Matrix3 J;
   87                 jacobian( params, verts, J );
   88                 double det = moab::Matrix3::determinant3( J );
   89                 if( fabs( det ) < 1.e-10 )
   90                 {
   91 #ifdef LINEAR_HEX_DEBUG
   92                     std::cerr << ss.str();
   93 #endif
   94 #ifndef LINEAR_HEX_DEBUG
   95                     std::cerr << x[0] << ", " << x[1] << ", " << x[2] << std::endl;
   96 #endif
   97                     std::cerr << "inverse solve failure: det: " << det << std::endl;
   98                     exit( -1 );
   99                 }
  100                 params -= moab::Matrix3::inverse( J, 1.0 / det ) * delta;
  101                 evaluate_forward( params, points, delta );
  102                 delta -= x;
  103             }
  104             return true;
  105         }
  106  
  107         void evaluate_forward( const CartVect& p, const CartVect* verts, CartVect& f ) const
  108         {
  109             typedef typename Points::value_type Vector;
  110             f.set( 0.0, 0.0, 0.0 );
  111             for( unsigned i = 0; i < 8; ++i )
  112             {
  113                 const double N_i = ( 1 + p[0] * reference_points( i )[0] ) * ( 1 + p[1] * reference_points( i )[1] ) *
  114                                    ( 1 + p[2] * reference_points( i )[2] );
  115                 f += N_i * verts[i];
  116             }
  117             f *= 0.125;
  118             return f;
  119         }
  120  
  121         double integrate_scalar_field( const CartVect* points, const double* field_values ) const {}
  122  
  123         template < typename Point, typename Points >
  124         Matrix3& jacobian( const Point& p, const Points& points, Matrix3& J ) const
  125         {
  126         }
  127  
  128       private:
  129     };  // Class Linear_hex_map
  130  
  131 }  // namespace element_utility
  132  
  133 }  // namespace moab
  134 #endif  // MOAB_LINEAR_HEX_nPP