Mesh Oriented datABase  (version 5.4.1)
Array-based unstructured mesh datastructure
QuadraticHex.cpp
Go to the documentation of this file.
00001 #include "moab/LocalDiscretization/QuadraticHex.hpp"
00002 #include "moab/Forward.hpp"
00003 
00004 namespace moab
00005 {
00006 
00007 // those are not just the corners, but for simplicity, keep this name
00008 //
00009 const int QuadraticHex::corner[27][3] = { { -1, -1, -1 }, { 1, -1, -1 }, { 1, 1, -1 },  // corner nodes: 0-7
00010                                           { -1, 1, -1 },                                // mid-edge nodes: 8-19
00011                                           { -1, -1, 1 },                 // center-face nodes 20-25  center node  26
00012                                           { 1, -1, 1 },                  //
00013                                           { 1, 1, 1 },    { -1, 1, 1 },  //                    4   ----- 19   -----  7
00014                                           { 0, -1, -1 },                 //                .   |                 .   |
00015                                           { 1, 0, -1 },                  //            16         25         18      |
00016                                           { 0, 1, -1 },                  //         .          |          .          |
00017                                           { -1, 0, -1 },                 //      5   ----- 17   -----  6             |
00018                                           { -1, -1, 0 },                 //      |            12       | 23         15
00019                                           { 1, -1, 0 },                  //      |                     |             |
00020                                           { 1, 1, 0 },                   //      |     20      |  26   |     22      |
00021                                           { -1, 1, 0 },                  //      |                     |             |
00022                                           { 0, -1, 1 },                  //     13         21  |      14             |
00023                                           { 1, 0, 1 },                   //      |             0   ----- 11   -----  3
00024                                           { 0, 1, 1 },                   //      |         .           |         .
00025                                           { -1, 0, 1 },                  //      |      8         24   |     10
00026                                           { 0, -1, 0 },                  //      |  .                  |  .
00027                                           { 1, 0, 0 },                   //      1   -----  9   -----  2
00028                                           { 0, 1, 0 },                   //
00029                                           { -1, 0, 0 },   { 0, 0, -1 },  { 0, 0, 1 },  { 0, 0, 0 } };
00030 
00031 double QuadraticHex::SH( const int i, const double params )
00032 {
00033     switch( i )
00034     {
00035         case -1:
00036             return ( params * params - params ) / 2;
00037         case 0:
00038             return 1 - params * params;
00039         case 1:
00040             return ( params * params + params ) / 2;
00041         default:
00042             return 0.;
00043     }
00044 }
00045 double QuadraticHex::DSH( const int i, const double params )
00046 {
00047     switch( i )
00048     {
00049         case -1:
00050             return params - 0.5;
00051         case 0:
00052             return -2 * params;
00053         case 1:
00054             return params + 0.5;
00055         default:
00056             return 0.;
00057     }
00058 }
00059 
00060 ErrorCode QuadraticHex::evalFcn( const double* params,
00061                                  const double* field,
00062                                  const int /*ndim*/,
00063                                  const int num_tuples,
00064                                  double* /*work*/,
00065                                  double* result )
00066 {
00067     assert( params && field && num_tuples > 0 );
00068     std::fill( result, result + num_tuples, 0.0 );
00069     for( int i = 0; i < 27; i++ )
00070     {
00071         const double sh = SH( corner[i][0], params[0] ) * SH( corner[i][1], params[1] ) * SH( corner[i][2], params[2] );
00072         for( int j = 0; j < num_tuples; j++ )
00073             result[j] += sh * field[num_tuples * i + j];
00074     }
00075 
00076     return MB_SUCCESS;
00077 }
00078 
00079 ErrorCode QuadraticHex::jacobianFcn( const double* params,
00080                                      const double* verts,
00081                                      const int nverts,
00082                                      const int ndim,
00083                                      double* /*work*/,
00084                                      double* result )
00085 {
00086     assert( 27 == nverts && params && verts );
00087     if( 27 != nverts ) return MB_FAILURE;
00088     Matrix3* J = reinterpret_cast< Matrix3* >( result );
00089     for( int i = 0; i < 27; i++ )
00090     {
00091         const double sh[3]  = { SH( corner[i][0], params[0] ), SH( corner[i][1], params[1] ),
00092                                SH( corner[i][2], params[2] ) };
00093         const double dsh[3] = { DSH( corner[i][0], params[0] ), DSH( corner[i][1], params[1] ),
00094                                 DSH( corner[i][2], params[2] ) };
00095 
00096         for( int j = 0; j < 3; j++ )
00097         {
00098             ( *J )( j, 0 ) += dsh[0] * sh[1] * sh[2] * verts[ndim * i + j];  // dxj/dr first column
00099             ( *J )( j, 1 ) += sh[0] * dsh[1] * sh[2] * verts[ndim * i + j];  // dxj/ds
00100             ( *J )( j, 2 ) += sh[0] * sh[1] * dsh[2] * verts[ndim * i + j];  // dxj/dt
00101         }
00102     }
00103 
00104     return MB_SUCCESS;
00105 }
00106 
00107 ErrorCode QuadraticHex::integrateFcn( const double* /*field*/,
00108                                       const double* /*verts*/,
00109                                       const int /*nverts*/,
00110                                       const int /*ndim*/,
00111                                       const int /*num_tuples*/,
00112                                       double* /*work*/,
00113                                       double* /*result*/ )
00114 {
00115     return MB_NOT_IMPLEMENTED;
00116 }
00117 
00118 ErrorCode QuadraticHex::reverseEvalFcn( EvalFcn eval,
00119                                         JacobianFcn jacob,
00120                                         InsideFcn ins,
00121                                         const double* posn,
00122                                         const double* verts,
00123                                         const int nverts,
00124                                         const int ndim,
00125                                         const double iter_tol,
00126                                         const double inside_tol,
00127                                         double* work,
00128                                         double* params,
00129                                         int* is_inside )
00130 {
00131     assert( posn && verts );
00132     return EvalSet::evaluate_reverse( eval, jacob, ins, posn, verts, nverts, ndim, iter_tol, inside_tol, work, params,
00133                                       is_inside );
00134 }
00135 
00136 int QuadraticHex::insideFcn( const double* params, const int ndim, const double tol )
00137 {
00138     return EvalSet::inside_function( params, ndim, tol );
00139 }
00140 
00141 ErrorCode QuadraticHex::normalFcn( const int /*ientDim*/,
00142                                    const int /*facet*/,
00143                                    const int /*nverts*/,
00144                                    const double* /*verts*/,
00145                                    double* /*normal[3]*/ )
00146 {
00147     return MB_NOT_IMPLEMENTED;
00148 }
00149 }  // namespace moab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines