MOAB: Mesh Oriented datABase  (version 5.4.1)
LinearQuad.cpp
Go to the documentation of this file.
00001 #include "moab/LocalDiscretization/LinearQuad.hpp"
00002 #include "moab/Matrix3.hpp"
00003 #include "moab/Forward.hpp"
00004 #include <cmath>
00005 #include <limits>
00006 
00007 namespace moab
00008 {
00009 
00010 const double LinearQuad::corner[4][2] = { { -1, -1 }, { 1, -1 }, { 1, 1 }, { -1, 1 } };
00011 
00012 /* For each point, its weight and location are stored as an array.
00013    Hence, the inner dimension is 2, the outer dimension is gauss_count.
00014    We use a one-point Gaussian quadrature, since it integrates linear functions exactly.
00015 */
00016 const double LinearQuad::gauss[1][2] = { { 2.0, 0.0 } };
00017 
00018 ErrorCode LinearQuad::jacobianFcn( const double* params,
00019                                    const double* verts,
00020                                    const int /*nverts*/,
00021                                    const int /*ndim*/,
00022                                    double*,
00023                                    double* result )
00024 {
00025     Matrix3* J = reinterpret_cast< Matrix3* >( result );
00026     *J         = Matrix3( 0.0 );
00027     for( unsigned i = 0; i < 4; ++i )
00028     {
00029         const double xi_p     = 1 + params[0] * corner[i][0];
00030         const double eta_p    = 1 + params[1] * corner[i][1];
00031         const double dNi_dxi  = corner[i][0] * eta_p;
00032         const double dNi_deta = corner[i][1] * xi_p;
00033         ( *J )( 0, 0 ) += dNi_dxi * verts[i * 3 + 0];
00034         ( *J )( 1, 0 ) += dNi_dxi * verts[i * 3 + 1];
00035         ( *J )( 0, 1 ) += dNi_deta * verts[i * 3 + 0];
00036         ( *J )( 1, 1 ) += dNi_deta * verts[i * 3 + 1];
00037     }
00038     ( *J ) *= 0.25;
00039     ( *J )( 2, 2 ) = 1.0; /* to make sure the Jacobian determinant is non-zero */
00040     return MB_SUCCESS;
00041 }  // LinearQuad::jacobian()
00042 
00043 ErrorCode LinearQuad::evalFcn( const double* params,
00044                                const double* field,
00045                                const int /*ndim*/,
00046                                const int num_tuples,
00047                                double*,
00048                                double* result )
00049 {
00050     for( int i = 0; i < num_tuples; i++ )
00051         result[i] = 0.0;
00052     for( unsigned i = 0; i < 4; ++i )
00053     {
00054         const double N_i = ( 1 + params[0] * corner[i][0] ) * ( 1 + params[1] * corner[i][1] );
00055         for( int j = 0; j < num_tuples; j++ )
00056             result[j] += N_i * field[i * num_tuples + j];
00057     }
00058     for( int i = 0; i < num_tuples; i++ )
00059         result[i] *= 0.25;
00060 
00061     return MB_SUCCESS;
00062 }
00063 
00064 ErrorCode LinearQuad::integrateFcn( const double* field,
00065                                     const double* verts,
00066                                     const int nverts,
00067                                     const int ndim,
00068                                     const int num_tuples,
00069                                     double* work,
00070                                     double* result )
00071 {
00072     double tmp_result[4];
00073     ErrorCode rval = MB_SUCCESS;
00074     for( int i = 0; i < num_tuples; i++ )
00075         result[i] = 0.0;
00076     CartVect x;
00077     Matrix3 J;
00078     for( unsigned int j1 = 0; j1 < LinearQuad::gauss_count; ++j1 )
00079     {
00080         x[0]      = LinearQuad::gauss[j1][1];
00081         double w1 = LinearQuad::gauss[j1][0];
00082         for( unsigned int j2 = 0; j2 < LinearQuad::gauss_count; ++j2 )
00083         {
00084             x[1]      = LinearQuad::gauss[j2][1];
00085             double w2 = LinearQuad::gauss[j2][0];
00086             rval      = evalFcn( x.array(), field, ndim, num_tuples, NULL, tmp_result );
00087             if( MB_SUCCESS != rval ) return rval;
00088             rval = jacobianFcn( x.array(), verts, nverts, ndim, work, J[0] );
00089             if( MB_SUCCESS != rval ) return rval;
00090             double tmp_det = w1 * w2 * J.determinant();
00091             for( int i = 0; i < num_tuples; i++ )
00092                 result[i] += tmp_result[i] * tmp_det;
00093         }
00094     }
00095     return MB_SUCCESS;
00096 }  // LinearHex::integrate_vector()
00097 
00098 ErrorCode LinearQuad::reverseEvalFcn( EvalFcn eval,
00099                                       JacobianFcn jacob,
00100                                       InsideFcn ins,
00101                                       const double* posn,
00102                                       const double* verts,
00103                                       const int nverts,
00104                                       const int ndim,
00105                                       const double iter_tol,
00106                                       const double inside_tol,
00107                                       double* work,
00108                                       double* params,
00109                                       int* is_inside )
00110 {
00111     return EvalSet::evaluate_reverse( eval, jacob, ins, posn, verts, nverts, ndim, iter_tol, inside_tol, work, params,
00112                                       is_inside );
00113 }
00114 
00115 int LinearQuad::insideFcn( const double* params, const int ndim, const double tol )
00116 {
00117     return EvalSet::inside_function( params, ndim, tol );
00118 }
00119 
00120 ErrorCode LinearQuad::normalFcn( const int ientDim,
00121                                  const int facet,
00122                                  const int nverts,
00123                                  const double* verts,
00124                                  double normal[3] )
00125 {
00126     // assert(facet <4 && ientDim == 1 && nverts==4);
00127     if( nverts != 4 ) MB_SET_ERR( MB_FAILURE, "Incorrect vertex count for passed quad :: expected value = 4" );
00128     if( ientDim != 1 ) MB_SET_ERR( MB_FAILURE, "Requesting normal for unsupported dimension :: expected value = 1 " );
00129     if( facet > 4 || facet < 0 ) MB_SET_ERR( MB_FAILURE, "Incorrect local edge id :: expected value = one of 0-3" );
00130 
00131     // Get the local vertex ids of  local edge
00132     int id0 = CN::mConnectivityMap[MBQUAD][ientDim - 1].conn[facet][0];
00133     int id1 = CN::mConnectivityMap[MBQUAD][ientDim - 1].conn[facet][1];
00134 
00135     // Find a vector along the edge
00136     double edge[3];
00137     for( int i = 0; i < 3; i++ )
00138     {
00139         edge[i] = verts[3 * id1 + i] - verts[3 * id0 + i];
00140     }
00141     // Find the normal of the face
00142     double x0[3], x1[3], fnrm[3];
00143     for( int i = 0; i < 3; i++ )
00144     {
00145         x0[i] = verts[3 * 1 + i] - verts[3 * 0 + i];
00146         x1[i] = verts[3 * 3 + i] - verts[3 * 0 + i];
00147     }
00148     fnrm[0] = x0[1] * x1[2] - x1[1] * x0[2];
00149     fnrm[1] = x1[0] * x0[2] - x0[0] * x1[2];
00150     fnrm[2] = x0[0] * x1[1] - x1[0] * x0[1];
00151 
00152     // Find the normal of the edge as the cross product of edge and face normal
00153 
00154     double a   = edge[1] * fnrm[2] - fnrm[1] * edge[2];
00155     double b   = edge[2] * fnrm[0] - fnrm[2] * edge[0];
00156     double c   = edge[0] * fnrm[1] - fnrm[0] * edge[1];
00157     double nrm = sqrt( a * a + b * b + c * c );
00158 
00159     if( nrm > std::numeric_limits< double >::epsilon() )
00160     {
00161         normal[0] = a / nrm;
00162         normal[1] = b / nrm;
00163         normal[2] = c / nrm;
00164     }
00165     return MB_SUCCESS;
00166 }
00167 
00168 }  // namespace moab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines