Mesh Oriented datABase  (version 5.4.1)
Array-based unstructured mesh datastructure
LinearTet.cpp
Go to the documentation of this file.
00001 #include "moab/LocalDiscretization/LinearTet.hpp"
00002 #include "moab/Forward.hpp"
00003 #include <algorithm>
00004 #include <cmath>
00005 #include <limits>
00006 
00007 namespace moab
00008 {
00009 
00010 const double LinearTet::corner[4][3] = { { 0, 0, 0 }, { 1, 0, 0 }, { 0, 1, 0 }, { 0, 0, 1 } };
00011 
00012 ErrorCode LinearTet::initFcn( const double* verts, const int nverts, double*& work )
00013 {
00014     // allocate work array as:
00015     // work[0..8] = T
00016     // work[9..17] = Tinv
00017     // work[18] = detT
00018     // work[19] = detTinv
00019     if( nverts != 4 )
00020     {
00021         std::cout << "Invalid Tetrahedron. Expected 4 vertices.\n";
00022         return MB_FAILURE;
00023     }
00024 
00025     assert( verts );
00026 
00027     Matrix3 J( verts[1 * 3 + 0] - verts[0 * 3 + 0], verts[2 * 3 + 0] - verts[0 * 3 + 0],
00028                verts[3 * 3 + 0] - verts[0 * 3 + 0], verts[1 * 3 + 1] - verts[0 * 3 + 1],
00029                verts[2 * 3 + 1] - verts[0 * 3 + 1], verts[3 * 3 + 1] - verts[0 * 3 + 1],
00030                verts[1 * 3 + 2] - verts[0 * 3 + 2], verts[2 * 3 + 2] - verts[0 * 3 + 2],
00031                verts[3 * 3 + 2] - verts[0 * 3 + 2] );
00032 
00033     // Update the work array
00034     if( !work ) work = new double[20];
00035 
00036     J.copyto( work );
00037     J.inverse().copyto( work + Matrix3::size );
00038     work[18] = J.determinant();
00039     work[19] = ( work[18] < 1e-12 ? std::numeric_limits< double >::max() : 1.0 / work[18] );
00040 
00041     return MB_SUCCESS;
00042 }
00043 
00044 ErrorCode LinearTet::evalFcn( const double* params,
00045                               const double* field,
00046                               const int /*ndim*/,
00047                               const int num_tuples,
00048                               double* /*work*/,
00049                               double* result )
00050 {
00051     assert( params && field && num_tuples > 0 );
00052     std::vector< double > f0( num_tuples );
00053     std::copy( field, field + num_tuples, f0.begin() );
00054     std::copy( field, field + num_tuples, result );
00055 
00056     for( unsigned i = 1; i < 4; ++i )
00057     {
00058         double p = 0.5 * ( params[i - 1] + 1 );  // transform from -1 <= p <= 1 to 0 <= p <= 1
00059         for( int j = 0; j < num_tuples; j++ )
00060             result[j] += ( field[i * num_tuples + j] - f0[j] ) * p;
00061     }
00062 
00063     return MB_SUCCESS;
00064 }
00065 
00066 ErrorCode LinearTet::integrateFcn( const double* field,
00067                                    const double* /*verts*/,
00068                                    const int nverts,
00069                                    const int /*ndim*/,
00070                                    const int num_tuples,
00071                                    double* work,
00072                                    double* result )
00073 {
00074     assert( field && num_tuples > 0 );
00075     std::fill( result, result + num_tuples, 0.0 );
00076     for( int i = 0; i < nverts; ++i )
00077     {
00078         for( int j = 0; j < num_tuples; j++ )
00079             result[j] += field[i * num_tuples + j];
00080     }
00081     double tmp = work[18] / 24.0;
00082     for( int i = 0; i < num_tuples; i++ )
00083         result[i] *= tmp;
00084 
00085     return MB_SUCCESS;
00086 }
00087 
00088 ErrorCode LinearTet::jacobianFcn( const double*, const double*, const int, const int, double* work, double* result )
00089 {
00090     // jacobian is cached in work array
00091     assert( work );
00092     std::copy( work, work + 9, result );
00093     return MB_SUCCESS;
00094 }
00095 
00096 ErrorCode LinearTet::reverseEvalFcn( EvalFcn eval,
00097                                      JacobianFcn jacob,
00098                                      InsideFcn ins,
00099                                      const double* posn,
00100                                      const double* verts,
00101                                      const int nverts,
00102                                      const int ndim,
00103                                      const double iter_tol,
00104                                      const double inside_tol,
00105                                      double* work,
00106                                      double* params,
00107                                      int* is_inside )
00108 {
00109     assert( posn && verts );
00110     return evaluate_reverse( eval, jacob, ins, posn, verts, nverts, ndim, iter_tol, inside_tol, work, params,
00111                              is_inside );
00112 }
00113 
00114 int LinearTet::insideFcn( const double* params, const int, const double tol )
00115 {
00116     return ( params[0] >= -1.0 - tol && params[1] >= -1.0 - tol && params[2] >= -1.0 - tol &&
00117              params[0] + params[1] + params[2] <= 1.0 + tol );
00118 }
00119 
00120 ErrorCode LinearTet::evaluate_reverse( EvalFcn eval,
00121                                        JacobianFcn jacob,
00122                                        InsideFcn inside_f,
00123                                        const double* posn,
00124                                        const double* verts,
00125                                        const int nverts,
00126                                        const int ndim,
00127                                        const double iter_tol,
00128                                        const double inside_tol,
00129                                        double* work,
00130                                        double* params,
00131                                        int* inside )
00132 {
00133     // TODO: should differentiate between epsilons used for
00134     // Newton Raphson iteration, and epsilons used for curved boundary geometry errors
00135     // right now, fix the tolerance used for NR
00136     const double error_tol_sqr = iter_tol * iter_tol;
00137     CartVect* cvparams         = reinterpret_cast< CartVect* >( params );
00138     const CartVect* cvposn     = reinterpret_cast< const CartVect* >( posn );
00139 
00140     // find best initial guess to improve convergence
00141     CartVect tmp_params[] = { CartVect( -1, -1, -1 ), CartVect( 1, -1, -1 ), CartVect( -1, 1, -1 ),
00142                               CartVect( -1, -1, 1 ) };
00143     double resl           = std::numeric_limits< double >::max();
00144     CartVect new_pos, tmp_pos;
00145     ErrorCode rval;
00146     for( unsigned int i = 0; i < 4; i++ )
00147     {
00148         rval = ( *eval )( tmp_params[i].array(), verts, ndim, ndim, work, tmp_pos.array() );
00149         if( MB_SUCCESS != rval ) return rval;
00150         double tmp_resl = ( tmp_pos - *cvposn ).length_squared();
00151         if( tmp_resl < resl )
00152         {
00153             *cvparams = tmp_params[i];
00154             new_pos   = tmp_pos;
00155             resl      = tmp_resl;
00156         }
00157     }
00158 
00159     // residual is diff between old and new pos; need to minimize that
00160     CartVect res = new_pos - *cvposn;
00161     Matrix3 J;
00162     rval = ( *jacob )( cvparams->array(), verts, nverts, ndim, work, J.array() );
00163 #ifndef NDEBUG
00164     double det = J.determinant();
00165     assert( det > std::numeric_limits< double >::epsilon() );
00166 #endif
00167     Matrix3 Ji = J.inverse();
00168 
00169     int iters = 0;
00170     // while |res| larger than tol
00171     int dum, *tmp_inside = ( inside ? inside : &dum );
00172     while( res % res > error_tol_sqr )
00173     {
00174         if( ++iters > 25 )
00175         {
00176             // if we haven't converged but we're outside, that's defined as success
00177             *tmp_inside = ( *inside_f )( params, ndim, inside_tol );
00178             if( !( *tmp_inside ) )
00179                 return MB_SUCCESS;
00180             else
00181                 return MB_INDEX_OUT_OF_RANGE;
00182         }
00183 
00184         // new params tries to eliminate residual
00185         *cvparams -= Ji * res;
00186 
00187         // get the new forward-evaluated position, and its difference from the target pt
00188         rval = ( *eval )( params, verts, ndim, ndim, work, new_pos.array() );
00189         if( MB_SUCCESS != rval ) return rval;
00190         res = new_pos - *cvposn;
00191     }
00192 
00193     if( inside ) *inside = ( *inside_f )( params, ndim, inside_tol );
00194 
00195     return MB_SUCCESS;
00196 }  // Map::evaluate_reverse()
00197 
00198 ErrorCode LinearTet::normalFcn( const int ientDim,
00199                                 const int facet,
00200                                 const int nverts,
00201                                 const double* verts,
00202                                 double normal[3] )
00203 {
00204     // assert(facet < 4 && ientDim == 2 && nverts == 4);
00205     if( nverts != 4 ) MB_SET_ERR( MB_FAILURE, "Incorrect vertex count for passed tet :: expected value = 4 " );
00206     if( ientDim != 2 ) MB_SET_ERR( MB_FAILURE, "Requesting normal for unsupported dimension :: expected value = 2 " );
00207     if( facet > 4 || facet < 0 ) MB_SET_ERR( MB_FAILURE, "Incorrect local face id :: expected value = one of 0-3" );
00208 
00209     int id0 = CN::mConnectivityMap[MBTET][ientDim - 1].conn[facet][0];
00210     int id1 = CN::mConnectivityMap[MBTET][ientDim - 1].conn[facet][1];
00211     int id2 = CN::mConnectivityMap[MBTET][ientDim - 1].conn[facet][2];
00212 
00213     double x0[3], x1[3];
00214 
00215     for( int i = 0; i < 3; i++ )
00216     {
00217         x0[i] = verts[3 * id1 + i] - verts[3 * id0 + i];
00218         x1[i] = verts[3 * id2 + i] - verts[3 * id0 + i];
00219     }
00220 
00221     double a   = x0[1] * x1[2] - x1[1] * x0[2];
00222     double b   = x1[0] * x0[2] - x0[0] * x1[2];
00223     double c   = x0[0] * x1[1] - x1[0] * x0[1];
00224     double nrm = sqrt( a * a + b * b + c * c );
00225 
00226     if( nrm > std::numeric_limits< double >::epsilon() )
00227     {
00228         normal[0] = a / nrm;
00229         normal[1] = b / nrm;
00230         normal[2] = c / nrm;
00231     }
00232     return MB_SUCCESS;
00233 }
00234 
00235 }  // namespace moab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines