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