QuadraticHex.cpp

Go to the documentation of this file.

#include "moab/LocalDiscretization/QuadraticHex.hpp"
#include "moab/Forward.hpp"

namespace moab
{

// those are not just the corners, but for simplicity, keep this name
//
const int QuadraticHex::corner[27][3] = { { -1, -1, -1 }, { 1, -1, -1 }, { 1, 1, -1 },  // corner nodes: 0-7
                                          { -1, 1, -1 },                                // mid-edge nodes: 8-19
                                          { -1, -1, 1 },                 // center-face nodes 20-25  center node  26
                                          { 1, -1, 1 },                  //
                                          { 1, 1, 1 },    { -1, 1, 1 },  //                    4   ----- 19   -----  7
                                          { 0, -1, -1 },                 //                .   |                 .   |
                                          { 1, 0, -1 },                  //            16         25         18      |
                                          { 0, 1, -1 },                  //         .          |          .          |
                                          { -1, 0, -1 },                 //      5   ----- 17   -----  6             |
                                          { -1, -1, 0 },                 //      |            12       | 23         15
                                          { 1, -1, 0 },                  //      |                     |             |
                                          { 1, 1, 0 },                   //      |     20      |  26   |     22      |
                                          { -1, 1, 0 },                  //      |                     |             |
                                          { 0, -1, 1 },                  //     13         21  |      14             |
                                          { 1, 0, 1 },                   //      |             0   ----- 11   -----  3
                                          { 0, 1, 1 },                   //      |         .           |         .
                                          { -1, 0, 1 },                  //      |      8         24   |     10
                                          { 0, -1, 0 },                  //      |  .                  |  .
                                          { 1, 0, 0 },                   //      1   -----  9   -----  2
                                          { 0, 1, 0 },                   //
                                          { -1, 0, 0 },   { 0, 0, -1 },  { 0, 0, 1 },  { 0, 0, 0 } };

double QuadraticHex::SH( const int i, const double params )
{
    switch( i )
    {
        case -1:
            return ( params * params - params ) / 2;
        case 0:
            return 1 - params * params;
        case 1:
            return ( params * params + params ) / 2;
        default:
            return 0.;
    }
}
double QuadraticHex::DSH( const int i, const double params )
{
    switch( i )
    {
        case -1:
            return params - 0.5;
        case 0:
            return -2 * params;
        case 1:
            return params + 0.5;
        default:
            return 0.;
    }
}

ErrorCode QuadraticHex::evalFcn( const double* params,
                                 const double* field,
                                 const int /*ndim*/,
                                 const int num_tuples,
                                 double* /*work*/,
                                 double* result )
{
    assert( params && field && num_tuples > 0 );
    std::fill( result, result + num_tuples, 0.0 );
    for( int i = 0; i < 27; i++ )
    {
        const double sh = SH( corner[i][0], params[0] ) * SH( corner[i][1], params[1] ) * SH( corner[i][2], params[2] );
        for( int j = 0; j < num_tuples; j++ )
            result[j] += sh * field[num_tuples * i + j];
    }

    return MB_SUCCESS;
}

ErrorCode QuadraticHex::jacobianFcn( const double* params,
                                     const double* verts,
                                     const int nverts,
                                     const int ndim,
                                     double* /*work*/,
                                     double* result )
{
    assert( 27 == nverts && params && verts );
    if( 27 != nverts ) return MB_FAILURE;
    Matrix3* J = reinterpret_cast< Matrix3* >( result );
    for( int i = 0; i < 27; i++ )
    {
        const double sh[3]  = { SH( corner[i][0], params[0] ), SH( corner[i][1], params[1] ),
                               SH( corner[i][2], params[2] ) };
        const double dsh[3] = { DSH( corner[i][0], params[0] ), DSH( corner[i][1], params[1] ),
                                DSH( corner[i][2], params[2] ) };

        for( int j = 0; j < 3; j++ )
        {
            ( *J )( j, 0 ) += dsh[0] * sh[1] * sh[2] * verts[ndim * i + j];  // dxj/dr first column
            ( *J )( j, 1 ) += sh[0] * dsh[1] * sh[2] * verts[ndim * i + j];  // dxj/ds
            ( *J )( j, 2 ) += sh[0] * sh[1] * dsh[2] * verts[ndim * i + j];  // dxj/dt
        }
    }

    return MB_SUCCESS;
}

ErrorCode QuadraticHex::integrateFcn( const double* /*field*/,
                                      const double* /*verts*/,
                                      const int /*nverts*/,
                                      const int /*ndim*/,
                                      const int /*num_tuples*/,
                                      double* /*work*/,
                                      double* /*result*/ )
{
    return MB_NOT_IMPLEMENTED;
}

ErrorCode QuadraticHex::reverseEvalFcn( EvalFcn eval,
                                        JacobianFcn jacob,
                                        InsideFcn ins,
                                        const double* posn,
                                        const double* verts,
                                        const int nverts,
                                        const int ndim,
                                        const double iter_tol,
                                        const double inside_tol,
                                        double* work,
                                        double* params,
                                        int* is_inside )
{
    assert( posn && verts );
    return EvalSet::evaluate_reverse( eval, jacob, ins, posn, verts, nverts, ndim, iter_tol, inside_tol, work, params,
                                      is_inside );
}

int QuadraticHex::insideFcn( const double* params, const int ndim, const double tol )
{
    return EvalSet::inside_function( params, ndim, tol );
}

ErrorCode QuadraticHex::normalFcn( const int /*ientDim*/,
                                   const int /*facet*/,
                                   const int /*nverts*/,
                                   const double* /*verts*/,
                                   double* /*normal[3]*/ )
{
    return MB_NOT_IMPLEMENTED;
}
}  // namespace moab