|
moab
|
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
1.7.6.1 
