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