MOAB: Mesh Oriented datABase  (version 5.2.1)
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, const double* field, const int /*ndim*/, const int num_tuples,
00061                                  double* /*work*/, double* result )
00062 {
00063     assert( params && field && num_tuples > 0 );
00064     std::fill( result, result + num_tuples, 0.0 );
00065     for( int i = 0; i < 27; i++ )
00066     {
00067         const double sh = SH( corner[i][0], params[0] ) * SH( corner[i][1], params[1] ) * SH( corner[i][2], params[2] );
00068         for( int j = 0; j < num_tuples; j++ )
00069             result[j] += sh * field[num_tuples * i + j];
00070     }
00071 
00072     return MB_SUCCESS;
00073 }
00074 
00075 ErrorCode QuadraticHex::jacobianFcn( const double* params, const double* verts, const int nverts, const int ndim,
00076                                      double* /*work*/, double* result )
00077 {
00078     assert( 27 == nverts && params && verts );
00079     if( 27 != nverts ) return MB_FAILURE;
00080     Matrix3* J = reinterpret_cast< Matrix3* >( result );
00081     for( int i = 0; i < 27; i++ )
00082     {
00083         const double sh[3]  = { SH( corner[i][0], params[0] ), SH( corner[i][1], params[1] ),
00084                                SH( corner[i][2], params[2] ) };
00085         const double dsh[3] = { DSH( corner[i][0], params[0] ), DSH( corner[i][1], params[1] ),
00086                                 DSH( corner[i][2], params[2] ) };
00087 
00088         for( int j = 0; j < 3; j++ )
00089         {
00090             ( *J )( j, 0 ) += dsh[0] * sh[1] * sh[2] * verts[ndim * i + j];  // dxj/dr first column
00091             ( *J )( j, 1 ) += sh[0] * dsh[1] * sh[2] * verts[ndim * i + j];  // dxj/ds
00092             ( *J )( j, 2 ) += sh[0] * sh[1] * dsh[2] * verts[ndim * i + j];  // dxj/dt
00093         }
00094     }
00095 
00096     return MB_SUCCESS;
00097 }
00098 
00099 ErrorCode QuadraticHex::integrateFcn( const double* /*field*/, const double* /*verts*/, const int /*nverts*/,
00100                                       const int /*ndim*/, const int /*num_tuples*/, double* /*work*/,
00101                                       double* /*result*/ )
00102 {
00103     return MB_NOT_IMPLEMENTED;
00104 }
00105 
00106 ErrorCode QuadraticHex::reverseEvalFcn( EvalFcn eval, JacobianFcn jacob, InsideFcn ins, const double* posn,
00107                                         const double* verts, const int nverts, const int ndim, const double iter_tol,
00108                                         const double inside_tol, double* work, double* params, int* is_inside )
00109 {
00110     assert( posn && verts );
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 QuadraticHex::insideFcn( const double* params, const int ndim, const double tol )
00116 {
00117     return EvalSet::inside_function( params, ndim, tol );
00118 }
00119 
00120 ErrorCode QuadraticHex::normalFcn( const int /*ientDim*/, const int /*facet*/, const int /*nverts*/,
00121                                    const double* /*verts*/, double* /*normal[3]*/ )
00122 {
00123     return MB_NOT_IMPLEMENTED;
00124 }
00125 }  // namespace moab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines