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