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