Branch data Line data Source code
1 : : #include "moab/LocalDiscretization/LinearQuad.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 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 : 2653 : ErrorCode LinearQuad::jacobianFcn( const double* params, const double* verts, const int /*nverts*/, const int /*ndim*/,
19 : : double*, double* result )
20 : : {
21 : 2653 : Matrix3* J = reinterpret_cast< Matrix3* >( result );
22 [ + - ]: 2653 : *J = Matrix3( 0.0 );
23 [ + + ]: 13265 : for( unsigned i = 0; i < 4; ++i )
24 : : {
25 : 10612 : const double xi_p = 1 + params[0] * corner[i][0];
26 : 10612 : const double eta_p = 1 + params[1] * corner[i][1];
27 : 10612 : const double dNi_dxi = corner[i][0] * eta_p;
28 : 10612 : const double dNi_deta = corner[i][1] * xi_p;
29 : 10612 : ( *J )( 0, 0 ) += dNi_dxi * verts[i * 3 + 0];
30 : 10612 : ( *J )( 1, 0 ) += dNi_dxi * verts[i * 3 + 1];
31 : 10612 : ( *J )( 0, 1 ) += dNi_deta * verts[i * 3 + 0];
32 : 10612 : ( *J )( 1, 1 ) += dNi_deta * verts[i * 3 + 1];
33 : : }
34 : 2653 : ( *J ) *= 0.25;
35 : 2653 : ( *J )( 2, 2 ) = 1.0; /* to make sure the Jacobian determinant is non-zero */
36 : 2653 : return MB_SUCCESS;
37 : : } // LinearQuad::jacobian()
38 : :
39 : 3984 : ErrorCode LinearQuad::evalFcn( const double* params, const double* field, const int /*ndim*/, const int num_tuples,
40 : : double*, double* result )
41 : : {
42 [ + + ]: 15934 : for( int i = 0; i < num_tuples; i++ )
43 : 11950 : result[i] = 0.0;
44 [ + + ]: 19920 : for( unsigned i = 0; i < 4; ++i )
45 : : {
46 : 15936 : const double N_i = ( 1 + params[0] * corner[i][0] ) * ( 1 + params[1] * corner[i][1] );
47 [ + + ]: 63736 : for( int j = 0; j < num_tuples; j++ )
48 : 47800 : result[j] += N_i * field[i * num_tuples + j];
49 : : }
50 [ + + ]: 15934 : for( int i = 0; i < num_tuples; i++ )
51 : 11950 : result[i] *= 0.25;
52 : :
53 : 3984 : return MB_SUCCESS;
54 : : }
55 : :
56 : 2 : ErrorCode LinearQuad::integrateFcn( const double* field, const double* verts, const int nverts, const int ndim,
57 : : const int num_tuples, double* work, double* result )
58 : : {
59 : : double tmp_result[4];
60 : 2 : ErrorCode rval = MB_SUCCESS;
61 [ + + ]: 6 : for( int i = 0; i < num_tuples; i++ )
62 : 4 : result[i] = 0.0;
63 [ + - ]: 2 : CartVect x;
64 [ + - ]: 2 : Matrix3 J;
65 [ + + ]: 4 : for( unsigned int j1 = 0; j1 < LinearQuad::gauss_count; ++j1 )
66 : : {
67 [ + - ]: 2 : x[0] = LinearQuad::gauss[j1][1];
68 : 2 : double w1 = LinearQuad::gauss[j1][0];
69 [ + + ]: 4 : for( unsigned int j2 = 0; j2 < LinearQuad::gauss_count; ++j2 )
70 : : {
71 [ + - ]: 2 : x[1] = LinearQuad::gauss[j2][1];
72 : 2 : double w2 = LinearQuad::gauss[j2][0];
73 [ + - ][ + - ]: 2 : rval = evalFcn( x.array(), field, ndim, num_tuples, NULL, tmp_result );
74 [ - + ]: 2 : if( MB_SUCCESS != rval ) return rval;
75 [ + - ][ + - ]: 2 : rval = jacobianFcn( x.array(), verts, nverts, ndim, work, J[0] );
[ + - ]
76 [ - + ]: 2 : if( MB_SUCCESS != rval ) return rval;
77 [ + - ]: 2 : double tmp_det = w1 * w2 * J.determinant();
78 [ + + ]: 6 : for( int i = 0; i < num_tuples; i++ )
79 : 4 : result[i] += tmp_result[i] * tmp_det;
80 : : }
81 : : }
82 : 2 : return MB_SUCCESS;
83 : : } // LinearHex::integrate_vector()
84 : :
85 : 1331 : ErrorCode LinearQuad::reverseEvalFcn( EvalFcn eval, JacobianFcn jacob, InsideFcn ins, const double* posn,
86 : : const double* verts, const int nverts, const int ndim, const double iter_tol,
87 : : const double inside_tol, double* work, double* params, int* is_inside )
88 : : {
89 : : return EvalSet::evaluate_reverse( eval, jacob, ins, posn, verts, nverts, ndim, iter_tol, inside_tol, work, params,
90 : 1331 : is_inside );
91 : : }
92 : :
93 : 2662 : int LinearQuad::insideFcn( const double* params, const int ndim, const double tol )
94 : : {
95 : 2662 : return EvalSet::inside_function( params, ndim, tol );
96 : : }
97 : :
98 : 8 : ErrorCode LinearQuad::normalFcn( const int ientDim, const int facet, const int nverts, const double* verts,
99 : : double normal[3] )
100 : : {
101 : : // assert(facet <4 && ientDim == 1 && nverts==4);
102 [ - + ][ # # ]: 8 : if( nverts != 4 ) MB_SET_ERR( MB_FAILURE, "Incorrect vertex count for passed quad :: expected value = 4" );
[ # # ][ # # ]
[ # # ][ # # ]
103 [ - + ][ # # ]: 8 : if( ientDim != 1 ) MB_SET_ERR( MB_FAILURE, "Requesting normal for unsupported dimension :: expected value = 1 " );
[ # # ][ # # ]
[ # # ][ # # ]
104 [ + - ][ - + ]: 8 : if( facet > 4 || facet < 0 ) MB_SET_ERR( MB_FAILURE, "Incorrect local edge id :: expected value = one of 0-3" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
105 : :
106 : : // Get the local vertex ids of local edge
107 : 8 : int id0 = CN::mConnectivityMap[MBQUAD][ientDim - 1].conn[facet][0];
108 : 8 : int id1 = CN::mConnectivityMap[MBQUAD][ientDim - 1].conn[facet][1];
109 : :
110 : : // Find a vector along the edge
111 : : double edge[3];
112 [ + + ]: 32 : for( int i = 0; i < 3; i++ )
113 : : {
114 : 24 : edge[i] = verts[3 * id1 + i] - verts[3 * id0 + i];
115 : : }
116 : : // Find the normal of the face
117 : : double x0[3], x1[3], fnrm[3];
118 [ + + ]: 32 : for( int i = 0; i < 3; i++ )
119 : : {
120 : 24 : x0[i] = verts[3 * 1 + i] - verts[3 * 0 + i];
121 : 24 : x1[i] = verts[3 * 3 + i] - verts[3 * 0 + i];
122 : : }
123 : 8 : fnrm[0] = x0[1] * x1[2] - x1[1] * x0[2];
124 : 8 : fnrm[1] = x1[0] * x0[2] - x0[0] * x1[2];
125 : 8 : fnrm[2] = x0[0] * x1[1] - x1[0] * x0[1];
126 : :
127 : : // Find the normal of the edge as the cross product of edge and face normal
128 : :
129 : 8 : double a = edge[1] * fnrm[2] - fnrm[1] * edge[2];
130 : 8 : double b = edge[2] * fnrm[0] - fnrm[2] * edge[0];
131 : 8 : double c = edge[0] * fnrm[1] - fnrm[0] * edge[1];
132 : 8 : double nrm = sqrt( a * a + b * b + c * c );
133 : :
134 [ + - ]: 8 : if( nrm > std::numeric_limits< double >::epsilon() )
135 : : {
136 : 8 : normal[0] = a / nrm;
137 : 8 : normal[1] = b / nrm;
138 : 8 : normal[2] = c / nrm;
139 : : }
140 : 8 : return MB_SUCCESS;
141 : : }
142 : :
143 [ + - ][ + - ]: 228 : } // namespace moab
|
