linear_hex_map.hpp

Go to the documentation of this file.

#ifndef MOAB_LINEAR_HEX_HPP
#define MOAB_LINEAR_HEX_HPP
 
#include "moab/Matrix3.hpp"
#include "moab/CartVect.hpp"
#include <sstream>
#include <iomanip>
#include <iostream>
 
namespace moab
{
 
namespace element_utility
{
 
 class Linear_hex_map
 {
 public:
 // Constructor
 Linear_hex_map() {}
 // Copy constructor
 Linear_hex_map( const Linear_hex_map& ) {}
 
 public:
 // Natural coordinates
 template < typename Moab, typename Entity_handle, typename Points, typename Point >
 std::pair< bool, CartVect > evaluate_reverse( const double* verts,
 const double* eval_point,
 const double tol = 1.e-6 ) const
 {
 CartVect params( 3, 0.0 );
 solve_inverse( eval_point, params, verts );
 bool point_found = solve_inverse( eval_point, params, verts, tol ) && is_contained( params, tol );
 return std::make_pair( point_found, params );
 }
 
 private:
 // This is a hack to avoid a .cpp file and C++11
 // reference_points(i,j) will be a 1 or -1;
 // This should unroll..
 inline const CartVect& reference_points( const std::size_t& i ) const
 {
 const CartVect rpts[8] = { CartVect( -1, -1, -1 ), CartVect( 1, -1, -1 ), CartVect( 1, 1, -1 ),
 CartVect( -1, 1, -1 ), CartVect( -1, -1, 1 ), CartVect( 1, -1, 1 ),
 CartVect( 1, 1, 1 ), CartVect( -1, 1, 1 ) };
 return rpts[i];
 }
 
 bool is_contained( const CartVect& params, const double tol ) const
 {
 // just look at the box+tol here
 return ( params[0] >= -1. - tol ) && ( params[0] <= 1. + tol ) && ( params[1] >= -1. - tol ) &&
 ( params[1] <= 1. + tol ) && ( params[2] >= -1. - tol ) && ( params[2] <= 1. + tol );
 }
 
 bool solve_inverse( const CartVect& point,
 CartVect& params,
 const CartVect* verts,
 const double tol = 1.e-6 ) const
 {
 const double error_tol_sqr = tol * tol;
 CartVect delta( 0.0, 0.0, 0.0 );
 params = delta;
 evaluate_forward( params, verts, delta );
 delta -= point;
 std::size_t num_iterations = 0;
#ifdef LINEAR_HEX_DEBUG
 std::stringstream ss;
 ss << "CartVect: ";
 ss << point[0] << ", " << point[1] << ", " << point[2] << std::endl;
 ss << "Hex: ";
 for( int i = 0; i < 8; ++i )
 ss << points[i][0] << ", " << points[i][1] << ", " << points[i][2] << std::endl;
 ss << std::endl;
#endif
 while( delta.length_squared() > error_tol_sqr )
 {
#ifdef LINEAR_HEX_DEBUG
 ss << "Iter #: " << num_iterations << " Err: " << delta.length() << " Iterate: ";
 ss << params[0] << ", " << params[1] << ", " << params[2] << std::endl;
#endif
 if( ++num_iterations >= 5 )
 {
 return false;
 }
 Matrix3 J;
 jacobian( params, verts, J );
 double det = moab::Matrix3::determinant3( J );
 if( fabs( det ) < 1.e-10 )
 {
#ifdef LINEAR_HEX_DEBUG
 std::cerr << ss.str();
#endif
#ifndef LINEAR_HEX_DEBUG
 std::cerr << x[0] << ", " << x[1] << ", " << x[2] << std::endl;
#endif
 std::cerr << "inverse solve failure: det: " << det << std::endl;
 exit( -1 );
 }
 params -= moab::Matrix3::inverse( J, 1.0 / det ) * delta;
 evaluate_forward( params, points, delta );
 delta -= x;
 }
 return true;
 }
 
 void evaluate_forward( const CartVect& p, const CartVect* verts, CartVect& f ) const
 {
 typedef typename Points::value_type Vector;
 f.set( 0.0, 0.0, 0.0 );
 for( unsigned i = 0; i < 8; ++i )
 {
 const double N_i = ( 1 + p[0] * reference_points( i )[0] ) * ( 1 + p[1] * reference_points( i )[1] ) *
 ( 1 + p[2] * reference_points( i )[2] );
 f += N_i * verts[i];
 }
 f *= 0.125;
 return f;
 }
 
 double integrate_scalar_field( const CartVect* points, const double* field_values ) const {}
 
 template < typename Point, typename Points >
 Matrix3& jacobian( const Point& p, const Points& points, Matrix3& J ) const
 {
 }
 
 private:
 }; // Class Linear_hex_map
 
} // namespace element_utility
 
} // namespace moab
#endif // MOAB_LINEAR_HEX_nPP