moab
|
00001 #include "moab/LinearHex.hpp" 00002 #include "moab/Matrix3.hpp" 00003 #include "moab/Forward.hpp" 00004 00005 namespace moab 00006 { 00007 00008 const double LinearHex::corner[8][3] = { { -1, -1, -1 }, 00009 { 1, -1, -1 }, 00010 { 1, 1, -1 }, 00011 { -1, 1, -1 }, 00012 { -1, -1, 1 }, 00013 { 1, -1, 1 }, 00014 { 1, 1, 1 }, 00015 { -1, 1, 1 } }; 00016 00017 /* For each point, its weight and location are stored as an array. 00018 Hence, the inner dimension is 2, the outer dimension is gauss_count. 00019 We use a one-point Gaussian quadrature, since it integrates linear functions exactly. 00020 */ 00021 const double LinearHex::gauss[1][2] = { { 2.0, 0.0 } }; 00022 00023 ErrorCode LinearHex::jacobianFcn(const double *params, const double *verts, const int /*nverts*/, const int ndim, 00024 double *, double *result) 00025 { 00026 assert(params && verts); 00027 Matrix3 *J = reinterpret_cast<Matrix3*>(result); 00028 *J = Matrix3(0.0); 00029 for (unsigned i = 0; i < 8; ++i) { 00030 const double params_p = 1 + params[0]*corner[i][0]; 00031 const double eta_p = 1 + params[1]*corner[i][1]; 00032 const double zeta_p = 1 + params[2]*corner[i][2]; 00033 const double dNi_dparams = corner[i][0] * eta_p * zeta_p; 00034 const double dNi_deta = corner[i][1] * params_p * zeta_p; 00035 const double dNi_dzeta = corner[i][2] * params_p * eta_p; 00036 (*J)(0,0) += dNi_dparams * verts[i*ndim+0]; 00037 (*J)(1,0) += dNi_dparams * verts[i*ndim+1]; 00038 (*J)(2,0) += dNi_dparams * verts[i*ndim+2]; 00039 (*J)(0,1) += dNi_deta * verts[i*ndim+0]; 00040 (*J)(1,1) += dNi_deta * verts[i*ndim+1]; 00041 (*J)(2,1) += dNi_deta * verts[i*ndim+2]; 00042 (*J)(0,2) += dNi_dzeta * verts[i*ndim+0]; 00043 (*J)(1,2) += dNi_dzeta * verts[i*ndim+1]; 00044 (*J)(2,2) += dNi_dzeta * verts[i*ndim+2]; 00045 } 00046 (*J) *= 0.125; 00047 return MB_SUCCESS; 00048 }// LinearHex::jacobian() 00049 00050 ErrorCode LinearHex::evalFcn(const double *params, const double *field, const int /*ndim*/, const int num_tuples, 00051 double *, double *result) 00052 { 00053 assert(params && field && num_tuples != -1); 00054 for (int i = 0; i < num_tuples; i++) result[i] = 0.0; 00055 for (unsigned i = 0; i < 8; ++i) { 00056 const double N_i = (1 + params[0]*corner[i][0]) 00057 * (1 + params[1]*corner[i][1]) 00058 * (1 + params[2]*corner[i][2]); 00059 for (int j = 0; j < num_tuples; j++) result[j] += N_i * field[i*num_tuples+j]; 00060 } 00061 for (int i = 0; i < num_tuples; i++) result[i] *= 0.125; 00062 00063 return MB_SUCCESS; 00064 } 00065 00066 ErrorCode LinearHex::integrateFcn(const double *field, const double *verts, const int nverts, const int ndim, 00067 const int num_tuples, double *work, double *result) 00068 { 00069 assert(field && verts && num_tuples != -1); 00070 double tmp_result[8]; 00071 ErrorCode rval = MB_SUCCESS; 00072 for (int i = 0; i < num_tuples; i++) result[i] = 0.0; 00073 CartVect x; 00074 Matrix3 J; 00075 for(unsigned int j1 = 0; j1 < LinearHex::gauss_count; ++j1) { 00076 x[0] = LinearHex::gauss[j1][1]; 00077 double w1 = LinearHex::gauss[j1][0]; 00078 for(unsigned int j2 = 0; j2 < LinearHex::gauss_count; ++j2) { 00079 x[1] = LinearHex::gauss[j2][1]; 00080 double w2 = LinearHex::gauss[j2][0]; 00081 for(unsigned int j3 = 0; j3 < LinearHex::gauss_count; ++j3) { 00082 x[2] = LinearHex::gauss[j3][1]; 00083 double w3 = LinearHex::gauss[j3][0]; 00084 rval = evalFcn(x.array(),field, ndim, num_tuples, NULL, tmp_result); 00085 if (MB_SUCCESS != rval) return rval; 00086 rval = jacobianFcn(x.array(), verts, nverts, ndim, work, J[0]); 00087 if (MB_SUCCESS != rval) return rval; 00088 double tmp_det = w1*w2*w3*J.determinant(); 00089 for (int i = 0; i < num_tuples; i++) result[i] += tmp_result[i]*tmp_det; 00090 } 00091 } 00092 } 00093 00094 return MB_SUCCESS; 00095 }// LinearHex::integrate_vector() 00096 00097 ErrorCode LinearHex::reverseEvalFcn(EvalFcn eval, JacobianFcn jacob, InsideFcn ins, 00098 const double *posn, const double *verts, const int nverts, const int ndim, 00099 const double iter_tol, const double inside_tol, double *work, 00100 double *params, int *is_inside) 00101 { 00102 assert(posn && verts); 00103 return EvalSet::evaluate_reverse(eval, jacob, ins, posn, verts, nverts, ndim, iter_tol, inside_tol, work, 00104 params, is_inside); 00105 } 00106 00107 int LinearHex::insideFcn(const double *params, const int ndim, const double tol) 00108 { 00109 return EvalSet::inside_function(params, ndim, tol); 00110 } 00111 00112 } // namespace moab