1 #include "moab/LocalDiscretization/LinearHex.hpp" 2 #include "moab/Matrix3.hpp" 3 #include "moab/Forward.hpp" 4 #include <math.h> 5 #include <limits> 6 7 namespace moab 8 { 9 10 const double LinearHex::corner[8][3] = { { -1, -1, -1 }, 11 { 1, -1, -1 }, 12 { 1, 1, -1 }, 13 { -1, 1, -1 }, 14 { -1, -1, 1 }, 15 { 1, -1, 1 }, 16 { 1, 1, 1 }, 17 { -1, 1, 1 } }; 18 19 /* For each point, its weight and location are stored as an array. 20 Hence, the inner dimension is 2, the outer dimension is gauss_count. 21 We use a one-point Gaussian quadrature, since it integrates linear functions exactly. 22 */ 23 const double LinearHex::gauss[1][2] = { { 2.0, 0.0 } }; 24 jacobianFcn(const double * params,const double * verts,const int,const int ndim,double *,double * result)25 ErrorCode LinearHex::jacobianFcn(const double *params, const double *verts, const int /*nverts*/, const int ndim, 26 double *, double *result) 27 { 28 assert(params && verts); 29 Matrix3 *J = reinterpret_cast<Matrix3*>(result); 30 *J = Matrix3(0.0); 31 for (unsigned i = 0; i < 8; ++i) { 32 const double params_p = 1 + params[0]*corner[i][0]; 33 const double eta_p = 1 + params[1]*corner[i][1]; 34 const double zeta_p = 1 + params[2]*corner[i][2]; 35 const double dNi_dparams = corner[i][0] * eta_p * zeta_p; 36 const double dNi_deta = corner[i][1] * params_p * zeta_p; 37 const double dNi_dzeta = corner[i][2] * params_p * eta_p; 38 (*J)(0,0) += dNi_dparams * verts[i*ndim+0]; 39 (*J)(1,0) += dNi_dparams * verts[i*ndim+1]; 40 (*J)(2,0) += dNi_dparams * verts[i*ndim+2]; 41 (*J)(0,1) += dNi_deta * verts[i*ndim+0]; 42 (*J)(1,1) += dNi_deta * verts[i*ndim+1]; 43 (*J)(2,1) += dNi_deta * verts[i*ndim+2]; 44 (*J)(0,2) += dNi_dzeta * verts[i*ndim+0]; 45 (*J)(1,2) += dNi_dzeta * verts[i*ndim+1]; 46 (*J)(2,2) += dNi_dzeta * verts[i*ndim+2]; 47 } 48 (*J) *= 0.125; 49 return MB_SUCCESS; 50 }// LinearHex::jacobian() 51 evalFcn(const double * params,const double * field,const int,const int num_tuples,double *,double * result)52 ErrorCode LinearHex::evalFcn(const double *params, const double *field, const int /*ndim*/, const int num_tuples, 53 double *, double *result) 54 { 55 assert(params && field && num_tuples != -1); 56 for (int i = 0; i < num_tuples; i++) result[i] = 0.0; 57 for (unsigned i = 0; i < 8; ++i) { 58 const double N_i = (1 + params[0]*corner[i][0]) 59 * (1 + params[1]*corner[i][1]) 60 * (1 + params[2]*corner[i][2]); 61 for (int j = 0; j < num_tuples; j++) result[j] += N_i * field[i*num_tuples+j]; 62 } 63 for (int i = 0; i < num_tuples; i++) result[i] *= 0.125; 64 65 return MB_SUCCESS; 66 } 67 integrateFcn(const double * field,const double * verts,const int nverts,const int ndim,const int num_tuples,double * work,double * result)68 ErrorCode LinearHex::integrateFcn(const double *field, const double *verts, const int nverts, const int ndim, 69 const int num_tuples, double *work, double *result) 70 { 71 assert(field && verts && num_tuples != -1); 72 double tmp_result[8]; 73 ErrorCode rval = MB_SUCCESS; 74 for (int i = 0; i < num_tuples; i++) result[i] = 0.0; 75 CartVect x; 76 Matrix3 J; 77 for(unsigned int j1 = 0; j1 < LinearHex::gauss_count; ++j1) { 78 x[0] = LinearHex::gauss[j1][1]; 79 double w1 = LinearHex::gauss[j1][0]; 80 for(unsigned int j2 = 0; j2 < LinearHex::gauss_count; ++j2) { 81 x[1] = LinearHex::gauss[j2][1]; 82 double w2 = LinearHex::gauss[j2][0]; 83 for(unsigned int j3 = 0; j3 < LinearHex::gauss_count; ++j3) { 84 x[2] = LinearHex::gauss[j3][1]; 85 double w3 = LinearHex::gauss[j3][0]; 86 rval = evalFcn(x.array(),field, ndim, num_tuples, NULL, tmp_result); 87 if (MB_SUCCESS != rval) return rval; 88 rval = jacobianFcn(x.array(), verts, nverts, ndim, work, J[0]); 89 if (MB_SUCCESS != rval) return rval; 90 double tmp_det = w1*w2*w3*J.determinant(); 91 for (int i = 0; i < num_tuples; i++) result[i] += tmp_result[i]*tmp_det; 92 } 93 } 94 } 95 96 return MB_SUCCESS; 97 }// LinearHex::integrate_vector() 98 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)99 ErrorCode LinearHex::reverseEvalFcn(EvalFcn eval, JacobianFcn jacob, InsideFcn ins, 100 const double *posn, const double *verts, const int nverts, const int ndim, 101 const double iter_tol, const double inside_tol, double *work, 102 double *params, int *is_inside) 103 { 104 assert(posn && verts); 105 return EvalSet::evaluate_reverse(eval, jacob, ins, posn, verts, nverts, ndim, iter_tol, inside_tol, work, 106 params, is_inside); 107 } 108 insideFcn(const double * params,const int ndim,const double tol)109 int LinearHex::insideFcn(const double *params, const int ndim, const double tol) 110 { 111 return EvalSet::inside_function(params, ndim, tol); 112 } 113 normalFcn(const int ientDim,const int facet,const int nverts,const double * verts,double normal[3])114 ErrorCode LinearHex::normalFcn(const int ientDim, const int facet, const int nverts, const double *verts, double normal[3]) 115 { 116 //assert(facet < 6 && ientDim == 2 && nverts == 8); 117 if (nverts != 8) 118 MB_SET_ERR(MB_FAILURE, "Incorrect vertex count for passed hex :: expected value = 8 "); 119 if (ientDim != 2) 120 MB_SET_ERR(MB_FAILURE, "Requesting normal for unsupported dimension :: expected value = 2 "); 121 if (facet >6 || facet < 0) 122 MB_SET_ERR(MB_FAILURE, "Incorrect local face id :: expected value = one of 0-5"); 123 124 int id0 = CN::mConnectivityMap[MBHEX][ientDim-1].conn[facet][0]; 125 int id1 = CN::mConnectivityMap[MBHEX][ientDim-1].conn[facet][1]; 126 int id2 = CN::mConnectivityMap[MBHEX][ientDim-1].conn[facet][3]; 127 128 double x0[3], x1[3]; 129 130 for (int i=0; i<3; i++) 131 { 132 x0[i] = verts[3*id1+i] - verts[3*id0+i]; 133 x1[i] = verts[3*id2+i] - verts[3*id0+i]; 134 } 135 136 double a = x0[1]*x1[2] - x1[1]*x0[2]; 137 double b = x1[0]*x0[2] - x0[0]*x1[2]; 138 double c = x0[0]*x1[1] - x1[0]*x0[1]; 139 double nrm = sqrt(a*a+b*b+c*c); 140 141 if (nrm > std::numeric_limits<double>::epsilon()) { 142 normal[0] = a/nrm; 143 normal[1] = b/nrm; 144 normal[2] = c/nrm; 145 } 146 return MB_SUCCESS; 147 } 148 149 150 } // namespace moab 151