moab
QuadraticHex.cpp
Go to the documentation of this file.
00001 #include "moab/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] = {
00010         { -1, -1, -1 },
00011         {  1, -1, -1 },
00012         {  1,  1, -1 },  // corner nodes: 0-7
00013         { -1,  1, -1 },  // mid-edge nodes: 8-19
00014         { -1, -1,  1 },  // center-face nodes 20-25  center node  26
00015         {  1, -1,  1 },  //
00016         {  1,  1,  1 },
00017         { -1,  1,  1 }, //                    4   ----- 19   -----  7
00018         {  0, -1, -1 }, //                .   |                 .   |
00019         {  1,  0, -1 }, //            16         25         18      |
00020         {  0,  1, -1 }, //         .          |          .          |
00021         { -1,  0, -1 }, //      5   ----- 17   -----  6             |
00022         { -1, -1,  0 }, //      |            12       | 23         15
00023         {  1, -1,  0 }, //      |                     |             |
00024         {  1,  1,  0 }, //      |     20      |  26   |     22      |
00025         { -1,  1,  0 }, //      |                     |             |
00026         {  0, -1,  1 }, //     13         21  |      14             |
00027         {  1,  0,  1 }, //      |             0   ----- 11   -----  3
00028         {  0,  1,  1 }, //      |         .           |         .
00029         { -1,  0,  1 }, //      |      8         24   |     10
00030         {  0, -1,  0 }, //      |  .                  |  .
00031         {  1,  0,  0 }, //      1   -----  9   -----  2
00032         {  0,  1,  0 }, //
00033         { -1,  0,  0 },
00034         {  0,  0, -1 },
00035         {  0,  0,  1 },
00036         {  0,  0,  0 }
00037     };
00038 
00039     double QuadraticHex::SH(const int i, const double params)
00040     {
00041       switch (i)
00042       {
00043         case -1: return (params*params-params)/2;
00044         case 0: return 1-params*params;
00045         case 1: return (params*params+params)/2;
00046         default: return 0.;
00047       }
00048     }
00049     double QuadraticHex::DSH(const int i, const double params)
00050     {
00051       switch (i)
00052       {
00053         case -1: return params-0.5;
00054         case 0: return -2*params;
00055         case 1: return params+0.5;
00056         default: 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]),
00084                              SH(corner[i][1], params[1]),
00085                              SH(corner[i][2], params[2]) };
00086         const double dsh[3]={ DSH(corner[i][0], params[0]),
00087                               DSH(corner[i][1], params[1]),
00088                               DSH(corner[i][2], params[2]) };
00089 
00090 
00091         for (int j=0; j<3; j++)
00092         {
00093           (*J)(j,0)+=dsh[0]*sh[1]*sh[2]*verts[ndim*i+j]; // dxj/dr first column
00094           (*J)(j,1)+=sh[0]*dsh[1]*sh[2]*verts[ndim*i+j]; // dxj/ds
00095           (*J)(j,2)+=sh[0]*sh[1]*dsh[2]*verts[ndim*i+j]; // dxj/dt
00096         }
00097       }
00098       
00099       return MB_SUCCESS;
00100     }
00101 
00102     ErrorCode QuadraticHex::integrateFcn(const double */*field*/, const double */*verts*/, const int /*nverts*/, const int /*ndim*/, const int /*num_tuples*/,
00103                                          double */*work*/, double */*result*/)
00104     {
00105       return MB_NOT_IMPLEMENTED;
00106     }
00107 
00108     ErrorCode QuadraticHex::reverseEvalFcn(EvalFcn eval, JacobianFcn jacob, InsideFcn ins, 
00109                                            const double *posn, const double *verts, const int nverts, const int ndim,
00110                                            const double iter_tol, const double inside_tol, double *work, 
00111                                            double *params, int *is_inside) 
00112     {
00113       assert(posn && verts);
00114       return EvalSet::evaluate_reverse(eval, jacob, ins, posn, verts, nverts, ndim, iter_tol, inside_tol, 
00115                                        work, params, is_inside);
00116     } 
00117 
00118     int QuadraticHex::insideFcn(const double *params, const int ndim, const double tol) 
00119     {
00120       return EvalSet::inside_function(params, ndim, tol);
00121     }
00122 } // namespace moab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines