 |
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
|
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
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
1.9.1.