Branch data Line data Source code
1 : : #include "moab/LocalDiscretization/LinearHex.hpp"
2 : : #include "moab/Matrix3.hpp"
3 : : #include "moab/Forward.hpp"
4 : : #include <math.h>
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 : 9580 : ErrorCode LinearHex::jacobianFcn( const double* params, const double* verts, const int /*nverts*/, const int ndim,
20 : : double*, double* result )
21 : : {
22 [ + - ][ - + ]: 9580 : assert( params && verts );
23 : 9580 : Matrix3* J = reinterpret_cast< Matrix3* >( result );
24 [ + - ]: 9580 : *J = Matrix3( 0.0 );
25 [ + + ]: 86220 : for( unsigned i = 0; i < 8; ++i )
26 : : {
27 : 76640 : const double params_p = 1 + params[0] * corner[i][0];
28 : 76640 : const double eta_p = 1 + params[1] * corner[i][1];
29 : 76640 : const double zeta_p = 1 + params[2] * corner[i][2];
30 : 76640 : const double dNi_dparams = corner[i][0] * eta_p * zeta_p;
31 : 76640 : const double dNi_deta = corner[i][1] * params_p * zeta_p;
32 : 76640 : const double dNi_dzeta = corner[i][2] * params_p * eta_p;
33 : 76640 : ( *J )( 0, 0 ) += dNi_dparams * verts[i * ndim + 0];
34 : 76640 : ( *J )( 1, 0 ) += dNi_dparams * verts[i * ndim + 1];
35 : 76640 : ( *J )( 2, 0 ) += dNi_dparams * verts[i * ndim + 2];
36 : 76640 : ( *J )( 0, 1 ) += dNi_deta * verts[i * ndim + 0];
37 : 76640 : ( *J )( 1, 1 ) += dNi_deta * verts[i * ndim + 1];
38 : 76640 : ( *J )( 2, 1 ) += dNi_deta * verts[i * ndim + 2];
39 : 76640 : ( *J )( 0, 2 ) += dNi_dzeta * verts[i * ndim + 0];
40 : 76640 : ( *J )( 1, 2 ) += dNi_dzeta * verts[i * ndim + 1];
41 : 76640 : ( *J )( 2, 2 ) += dNi_dzeta * verts[i * ndim + 2];
42 : : }
43 : 9580 : ( *J ) *= 0.125;
44 : 9580 : return MB_SUCCESS;
45 : : } // LinearHex::jacobian()
46 : :
47 : 17828 : ErrorCode LinearHex::evalFcn( const double* params, const double* field, const int /*ndim*/, const int num_tuples,
48 : : double*, double* result )
49 : : {
50 [ + - ][ + - ]: 17828 : assert( params && field && num_tuples != -1 );
[ - + ]
51 [ + + ]: 71310 : for( int i = 0; i < num_tuples; i++ )
52 : 53482 : result[i] = 0.0;
53 [ + + ]: 160452 : for( unsigned i = 0; i < 8; ++i )
54 : : {
55 : : const double N_i =
56 : 142624 : ( 1 + params[0] * corner[i][0] ) * ( 1 + params[1] * corner[i][1] ) * ( 1 + params[2] * corner[i][2] );
57 [ + + ]: 570480 : for( int j = 0; j < num_tuples; j++ )
58 : 427856 : result[j] += N_i * field[i * num_tuples + j];
59 : : }
60 [ + + ]: 71310 : for( int i = 0; i < num_tuples; i++ )
61 : 53482 : result[i] *= 0.125;
62 : :
63 : 17828 : return MB_SUCCESS;
64 : : }
65 : :
66 : 2 : ErrorCode LinearHex::integrateFcn( const double* field, const double* verts, const int nverts, const int ndim,
67 : : const int num_tuples, double* work, double* result )
68 : : {
69 [ + - ][ + - ]: 2 : assert( field && verts && num_tuples != -1 );
[ - + ]
70 : : double tmp_result[8];
71 : 2 : ErrorCode rval = MB_SUCCESS;
72 [ + + ]: 6 : for( int i = 0; i < num_tuples; i++ )
73 : 4 : result[i] = 0.0;
74 [ + - ]: 2 : CartVect x;
75 [ + - ]: 2 : Matrix3 J;
76 [ + + ]: 4 : for( unsigned int j1 = 0; j1 < LinearHex::gauss_count; ++j1 )
77 : : {
78 [ + - ]: 2 : x[0] = LinearHex::gauss[j1][1];
79 : 2 : double w1 = LinearHex::gauss[j1][0];
80 [ + + ]: 4 : for( unsigned int j2 = 0; j2 < LinearHex::gauss_count; ++j2 )
81 : : {
82 [ + - ]: 2 : x[1] = LinearHex::gauss[j2][1];
83 : 2 : double w2 = LinearHex::gauss[j2][0];
84 [ + + ]: 4 : for( unsigned int j3 = 0; j3 < LinearHex::gauss_count; ++j3 )
85 : : {
86 [ + - ]: 2 : x[2] = LinearHex::gauss[j3][1];
87 : 2 : double w3 = LinearHex::gauss[j3][0];
88 [ + - ][ + - ]: 2 : rval = evalFcn( x.array(), field, ndim, num_tuples, NULL, tmp_result );
89 [ - + ]: 2 : if( MB_SUCCESS != rval ) return rval;
90 [ + - ][ + - ]: 2 : rval = jacobianFcn( x.array(), verts, nverts, ndim, work, J[0] );
[ + - ]
91 [ - + ]: 2 : if( MB_SUCCESS != rval ) return rval;
92 [ + - ]: 2 : double tmp_det = w1 * w2 * w3 * J.determinant();
93 [ + + ]: 6 : for( int i = 0; i < num_tuples; i++ )
94 : 4 : result[i] += tmp_result[i] * tmp_det;
95 : : }
96 : : }
97 : : }
98 : :
99 : 2 : return MB_SUCCESS;
100 : : } // LinearHex::integrate_vector()
101 : :
102 : 8248 : ErrorCode LinearHex::reverseEvalFcn( EvalFcn eval, JacobianFcn jacob, InsideFcn ins, const double* posn,
103 : : const double* verts, const int nverts, const int ndim, const double iter_tol,
104 : : const double inside_tol, double* work, double* params, int* is_inside )
105 : : {
106 [ + - ][ - + ]: 8248 : assert( posn && verts );
107 : : return EvalSet::evaluate_reverse( eval, jacob, ins, posn, verts, nverts, ndim, iter_tol, inside_tol, work, params,
108 : 8248 : is_inside );
109 : : }
110 : :
111 : 9579 : int LinearHex::insideFcn( const double* params, const int ndim, const double tol )
112 : : {
113 : 9579 : return EvalSet::inside_function( params, ndim, tol );
114 : : }
115 : :
116 : 8 : ErrorCode LinearHex::normalFcn( const int ientDim, const int facet, const int nverts, const double* verts,
117 : : double normal[3] )
118 : : {
119 : : // assert(facet < 6 && ientDim == 2 && nverts == 8);
120 [ - + ][ # # ]: 8 : if( nverts != 8 ) MB_SET_ERR( MB_FAILURE, "Incorrect vertex count for passed hex :: expected value = 8 " );
[ # # ][ # # ]
[ # # ][ # # ]
121 [ - + ][ # # ]: 8 : if( ientDim != 2 ) MB_SET_ERR( MB_FAILURE, "Requesting normal for unsupported dimension :: expected value = 2 " );
[ # # ][ # # ]
[ # # ][ # # ]
122 [ + - ][ - + ]: 8 : if( facet > 6 || facet < 0 ) MB_SET_ERR( MB_FAILURE, "Incorrect local face id :: expected value = one of 0-5" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
123 : :
124 : 8 : int id0 = CN::mConnectivityMap[MBHEX][ientDim - 1].conn[facet][0];
125 : 8 : int id1 = CN::mConnectivityMap[MBHEX][ientDim - 1].conn[facet][1];
126 : 8 : int id2 = CN::mConnectivityMap[MBHEX][ientDim - 1].conn[facet][3];
127 : :
128 : : double x0[3], x1[3];
129 : :
130 [ + + ]: 32 : for( int i = 0; i < 3; i++ )
131 : : {
132 : 24 : x0[i] = verts[3 * id1 + i] - verts[3 * id0 + i];
133 : 24 : x1[i] = verts[3 * id2 + i] - verts[3 * id0 + i];
134 : : }
135 : :
136 : 8 : double a = x0[1] * x1[2] - x1[1] * x0[2];
137 : 8 : double b = x1[0] * x0[2] - x0[0] * x1[2];
138 : 8 : double c = x0[0] * x1[1] - x1[0] * x0[1];
139 : 8 : double nrm = sqrt( a * a + b * b + c * c );
140 : :
141 [ + - ]: 8 : if( nrm > std::numeric_limits< double >::epsilon() )
142 : : {
143 : 8 : normal[0] = a / nrm;
144 : 8 : normal[1] = b / nrm;
145 : 8 : normal[2] = c / nrm;
146 : : }
147 : 8 : return MB_SUCCESS;
148 : : }
149 : :
150 [ + - ][ + - ]: 228 : } // namespace moab
|
