moab
LinearQuad.cpp
Go to the documentation of this file.
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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines