Branch data Line data Source code
1 : : #include "moab/LocalDiscretization/LinearTri.hpp"
2 : : #include "moab/Forward.hpp"
3 : : #include <algorithm>
4 : : #include <math.h>
5 : : #include <limits>
6 : :
7 : : namespace moab
8 : : {
9 : :
10 : : const double LinearTri::corner[3][2] = { { 0, 0 }, { 1, 0 }, { 0, 1 } };
11 : :
12 : 6 : ErrorCode LinearTri::initFcn( const double* verts, const int nverts, double*& work )
13 : : {
14 : : // allocate work array as:
15 : : // work[0..8] = T
16 : : // work[9..17] = Tinv
17 : : // work[18] = detT
18 : : // work[19] = detTinv
19 [ + - ][ - + ]: 6 : assert( nverts == 3 && verts );
20 [ + + ][ + - ]: 6 : if( !work ) work = new double[20];
21 : :
22 : 12 : Matrix3 J( verts[1 * 3 + 0] - verts[0 * 3 + 0], verts[2 * 3 + 0] - verts[0 * 3 + 0], 0.0,
23 : 24 : verts[1 * 3 + 1] - verts[0 * 3 + 1], verts[2 * 3 + 1] - verts[0 * 3 + 1], 0.0,
24 [ + - ]: 6 : verts[1 * 3 + 2] - verts[0 * 3 + 2], verts[2 * 3 + 2] - verts[0 * 3 + 2], 1.0 );
25 [ + - ]: 6 : J *= 0.5;
26 : :
27 [ + - ]: 6 : J.copyto( work );
28 [ + - ][ + - ]: 6 : J.inverse().copyto( work + Matrix3::size );
29 [ + - ]: 6 : work[18] = J.determinant();
30 [ - + ]: 6 : work[19] = ( work[18] < 1e-12 ? std::numeric_limits< double >::max() : 1.0 / work[18] );
31 : :
32 : 6 : return MB_SUCCESS;
33 : : }
34 : :
35 : 5797 : ErrorCode LinearTri::evalFcn( const double* params, const double* field, const int /*ndim*/, const int num_tuples,
36 : : double* /*work*/, double* result )
37 : : {
38 [ + - ][ + - ]: 5797 : assert( params && field && num_tuples > 0 );
[ - + ]
39 : : // convert to [0,1]
40 : 5797 : double p1 = 0.5 * ( 1.0 + params[0] ), p2 = 0.5 * ( 1.0 + params[1] ), p0 = 1.0 - p1 - p2;
41 : :
42 [ + + ]: 23188 : for( int j = 0; j < num_tuples; j++ )
43 : 17391 : result[j] = p0 * field[0 * num_tuples + j] + p1 * field[1 * num_tuples + j] + p2 * field[2 * num_tuples + j];
44 : :
45 : 5797 : return MB_SUCCESS;
46 : : }
47 : :
48 : 2 : ErrorCode LinearTri::integrateFcn( const double* field, const double* /*verts*/, const int nverts, const int /*ndim*/,
49 : : const int num_tuples, double* work, double* result )
50 : : {
51 [ + - ][ - + ]: 2 : assert( field && num_tuples > 0 );
52 [ + - ]: 2 : std::fill( result, result + num_tuples, 0.0 );
53 [ + + ]: 8 : for( int i = 0; i < nverts; ++i )
54 : : {
55 [ + + ]: 18 : for( int j = 0; j < num_tuples; j++ )
56 : 12 : result[j] += field[i * num_tuples + j];
57 : : }
58 : 2 : double tmp = work[18] / 6.0;
59 [ + + ]: 6 : for( int i = 0; i < num_tuples; i++ )
60 : 4 : result[i] *= tmp;
61 : :
62 : 2 : return MB_SUCCESS;
63 : : }
64 : :
65 : 2332 : ErrorCode LinearTri::jacobianFcn( const double*, const double*, const int, const int, double* work, double* result )
66 : : {
67 : : // jacobian is cached in work array
68 [ - + ]: 2332 : assert( work );
69 : 2332 : std::copy( work, work + 9, result );
70 : 2332 : return MB_SUCCESS;
71 : : }
72 : :
73 : 1166 : ErrorCode LinearTri::reverseEvalFcn( EvalFcn eval, JacobianFcn jacob, InsideFcn ins, const double* posn,
74 : : const double* verts, const int nverts, const int ndim, const double iter_tol,
75 : : const double inside_tol, double* work, double* params, int* is_inside )
76 : : {
77 [ + - ][ - + ]: 1166 : assert( posn && verts );
78 : : return evaluate_reverse( eval, jacob, ins, posn, verts, nverts, ndim, iter_tol, inside_tol, work, params,
79 : 1166 : is_inside );
80 : : }
81 : :
82 : 2497 : int LinearTri::insideFcn( const double* params, const int, const double tol )
83 : : {
84 [ + - ][ + - ]: 2497 : return ( params[0] >= -1.0 - tol && params[1] >= -1.0 - tol && params[0] + params[1] <= 1.0 + tol );
[ + + ]
85 : : }
86 : :
87 : 1166 : ErrorCode LinearTri::evaluate_reverse( EvalFcn eval, JacobianFcn jacob, InsideFcn inside_f, const double* posn,
88 : : const double* verts, const int nverts, const int ndim, const double iter_tol,
89 : : const double inside_tol, double* work, double* params, int* inside )
90 : : {
91 : : // TODO: should differentiate between epsilons used for
92 : : // Newton Raphson iteration, and epsilons used for curved boundary geometry errors
93 : : // right now, fix the tolerance used for NR
94 : 1166 : const double error_tol_sqr = iter_tol * iter_tol;
95 : 1166 : CartVect* cvparams = reinterpret_cast< CartVect* >( params );
96 : 1166 : const CartVect* cvposn = reinterpret_cast< const CartVect* >( posn );
97 : :
98 : : // find best initial guess to improve convergence
99 [ + - ][ + - ]: 1166 : CartVect tmp_params[] = { CartVect( -1, -1, -1 ), CartVect( 1, -1, -1 ), CartVect( -1, 1, -1 ) };
[ + - ]
100 : 1166 : double resl = std::numeric_limits< double >::max();
101 [ + - ][ + - ]: 1166 : CartVect new_pos, tmp_pos;
102 : : ErrorCode rval;
103 [ + + ]: 4664 : for( unsigned int i = 0; i < 3; i++ )
104 : : {
105 [ + - ][ + - ]: 3498 : rval = ( *eval )( tmp_params[i].array(), verts, ndim, 3, work, tmp_pos.array() );
[ + - ]
106 [ - + ]: 3498 : if( MB_SUCCESS != rval ) return rval;
107 [ + - ][ + - ]: 3498 : double tmp_resl = ( tmp_pos - *cvposn ).length_squared();
108 [ + + ]: 3498 : if( tmp_resl < resl )
109 : : {
110 : 1980 : *cvparams = tmp_params[i];
111 : 1980 : new_pos = tmp_pos;
112 : 1980 : resl = tmp_resl;
113 : : }
114 : : }
115 : :
116 : : // residual is diff between old and new pos; need to minimize that
117 [ + - ]: 1166 : CartVect res = new_pos - *cvposn;
118 [ + - ]: 1166 : Matrix3 J;
119 [ + - ][ + - ]: 1166 : rval = ( *jacob )( cvparams->array(), verts, nverts, ndim, work, J[0] );
[ + - ]
120 : : #ifndef NDEBUG
121 [ + - ]: 1166 : double det = J.determinant();
122 [ - + ]: 1166 : assert( det > std::numeric_limits< double >::epsilon() );
123 : : #endif
124 [ + - ]: 1166 : Matrix3 Ji = J.inverse();
125 : :
126 : 1166 : int iters = 0;
127 : : // while |res| larger than tol
128 [ + - ][ + + ]: 2299 : while( res % res > error_tol_sqr )
129 : : {
130 [ - + ]: 1133 : if( ++iters > 25 ) return MB_FAILURE;
131 : :
132 : : // new params tries to eliminate residual
133 [ + - ][ + - ]: 1133 : *cvparams -= Ji * res;
134 : :
135 : : // get the new forward-evaluated position, and its difference from the target pt
136 [ + - ][ + - ]: 1133 : rval = ( *eval )( params, verts, ndim, 3, work, new_pos.array() );
137 [ - + ]: 1133 : if( MB_SUCCESS != rval ) return rval;
138 [ + - ]: 1133 : res = new_pos - *cvposn;
139 : : }
140 : :
141 [ + - ][ + - ]: 1166 : if( inside ) *inside = ( *inside_f )( params, ndim, inside_tol );
142 : :
143 : 1166 : return MB_SUCCESS;
144 : : } // Map::evaluate_reverse()
145 : :
146 : : /* ErrorCode LinearTri::get_normal( int facet, double *work, double *normal)
147 : : {
148 : : ErrorCode error;
149 : : //Get the local vertex ids of local edge
150 : : int id1 = ledges[facet][0];
151 : : int id2 = ledges[facet][1];
152 : :
153 : : //Find the normal to the face
154 : : double face_normal[3];
155 : :
156 : :
157 : : }*/
158 : :
159 : 8 : ErrorCode LinearTri::normalFcn( const int ientDim, const int facet, const int nverts, const double* verts,
160 : : double normal[3] )
161 : : {
162 : : // assert(facet < 3 && ientDim == 1 && nverts==3);
163 [ - + ][ # # ]: 8 : if( nverts != 3 ) MB_SET_ERR( MB_FAILURE, "Incorrect vertex count for passed triangle :: expected value = 3 " );
[ # # ][ # # ]
[ # # ][ # # ]
164 [ - + ][ # # ]: 8 : if( ientDim != 1 ) MB_SET_ERR( MB_FAILURE, "Requesting normal for unsupported dimension :: expected value = 1 " );
[ # # ][ # # ]
[ # # ][ # # ]
165 [ + - ][ - + ]: 8 : if( facet > 3 || facet < 0 ) MB_SET_ERR( MB_FAILURE, "Incorrect local edge id :: expected value = one of 0-2" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
166 : :
167 : : // Get the local vertex ids of local edge
168 : 8 : int id0 = CN::mConnectivityMap[MBTRI][ientDim - 1].conn[facet][0];
169 : 8 : int id1 = CN::mConnectivityMap[MBTRI][ientDim - 1].conn[facet][1];
170 : :
171 : : // Find a vector along the edge
172 : : double edge[3];
173 [ + + ]: 32 : for( int i = 0; i < 3; i++ )
174 : : {
175 : 24 : edge[i] = verts[3 * id1 + i] - verts[3 * id0 + i];
176 : : }
177 : : // Find the normal of the face
178 : : double x0[3], x1[3], fnrm[3];
179 [ + + ]: 32 : for( int i = 0; i < 3; i++ )
180 : : {
181 : 24 : x0[i] = verts[3 * 1 + i] - verts[3 * 0 + i];
182 : 24 : x1[i] = verts[3 * 2 + i] - verts[3 * 0 + i];
183 : : }
184 : 8 : fnrm[0] = x0[1] * x1[2] - x1[1] * x0[2];
185 : 8 : fnrm[1] = x1[0] * x0[2] - x0[0] * x1[2];
186 : 8 : fnrm[2] = x0[0] * x1[1] - x1[0] * x0[1];
187 : :
188 : : // Find the normal of the edge as the cross product of edge and face normal
189 : :
190 : 8 : double a = edge[1] * fnrm[2] - fnrm[1] * edge[2];
191 : 8 : double b = edge[2] * fnrm[0] - fnrm[2] * edge[0];
192 : 8 : double c = edge[0] * fnrm[1] - fnrm[0] * edge[1];
193 : 8 : double nrm = sqrt( a * a + b * b + c * c );
194 : :
195 [ + - ]: 8 : if( nrm > std::numeric_limits< double >::epsilon() )
196 : : {
197 : 8 : normal[0] = a / nrm;
198 : 8 : normal[1] = b / nrm;
199 : 8 : normal[2] = c / nrm;
200 : : }
201 : 8 : return MB_SUCCESS;
202 : : }
203 : :
204 [ + - ][ + - ]: 228 : } // namespace moab
|
