Branch data Line data Source code
1 : : /* *****************************************************************
2 : : MESQUITE -- The Mesh Quality Improvement Toolkit
3 : :
4 : : Copyright 2006 Sandia National Laboratories. Developed at the
5 : : University of Wisconsin--Madison under SNL contract number
6 : : 624796. The U.S. Government and the University of Wisconsin
7 : : retain certain rights to this software.
8 : :
9 : : This library is free software; you can redistribute it and/or
10 : : modify it under the terms of the GNU Lesser General Public
11 : : License as published by the Free Software Foundation; either
12 : : version 2.1 of the License, or (at your option) any later version.
13 : :
14 : : This library is distributed in the hope that it will be useful,
15 : : but WITHOUT ANY WARRANTY; without even the implied warranty of
16 : : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 : : Lesser General Public License for more details.
18 : :
19 : : You should have received a copy of the GNU Lesser General Public License
20 : : (lgpl.txt) along with this library; if not, write to the Free Software
21 : : Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
22 : :
23 : : (2006) [email protected]
24 : :
25 : : ***************************************************************** */
26 : :
27 : : /** \file TMPQualityMetric.cpp
28 : : * \brief
29 : : * \author Jason Kraftcheck
30 : : */
31 : :
32 : : #undef PRINT_INFO
33 : :
34 : : #include "Mesquite.hpp"
35 : : #include "TMPQualityMetric.hpp"
36 : : #include "MsqMatrix.hpp"
37 : : #include "ElementQM.hpp"
38 : : #include "MsqError.hpp"
39 : : #include "Vector3D.hpp"
40 : : #include "PatchData.hpp"
41 : : #include "MappingFunction.hpp"
42 : : #include "WeightCalculator.hpp"
43 : : #include "TargetCalculator.hpp"
44 : : #include "TargetMetricUtil.hpp"
45 : :
46 : : #ifdef PRINT_INFO
47 : : #include <iostream>
48 : : #endif
49 : :
50 : : #include <functional>
51 : : #include <algorithm>
52 : :
53 : : namespace MBMesquite
54 : : {
55 : :
56 : 436593 : int TMPQualityMetric::get_negate_flag() const
57 : : {
58 : 436593 : return 1;
59 : : }
60 : :
61 : 216552 : void TMPQualityMetric::get_evaluations( PatchData& pd, std::vector< size_t >& handles, bool free, MsqError& err )
62 : : {
63 : 216552 : get_sample_pt_evaluations( pd, handles, free, err );
64 : 216552 : }
65 : :
66 : 23 : void TMPQualityMetric::get_patch_evaluations( PatchData& pd, std::vector< size_t >& handles, bool free, MsqError& err )
67 : : {
68 : 23 : get_sample_pt_evaluations( pd, handles, free, err );
69 : 23 : }
70 : :
71 : 5267379 : void TMPQualityMetric::get_element_evaluations( PatchData& pd, size_t p_elem, std::vector< size_t >& handles,
72 : : MsqError& err )
73 : : {
74 : 5267379 : get_elem_sample_points( pd, p_elem, handles, err );
75 : 5267379 : }
76 : :
77 : 9767364 : bool TMPQualityMetric::evaluate( PatchData& pd, size_t p_handle, double& value, MsqError& err )
78 : : {
79 : : size_t num_idx;
80 [ + - ]: 9767364 : bool valid = evaluate_internal( pd, p_handle, value, mIndices, num_idx, err );
81 [ + - ][ + + ]: 9767364 : if( MSQ_CHKERR( err ) || !valid ) return false;
[ + - ][ - + ]
[ + + ][ + + ]
82 : :
83 : : // apply target weight to value
84 [ + + ]: 9766133 : if( weightCalc )
85 : : {
86 [ + - ]: 437678 : const Sample s = ElemSampleQM::sample( p_handle );
87 [ + - ]: 437678 : const size_t e = ElemSampleQM::elem( p_handle );
88 [ + - ]: 437678 : double ck = weightCalc->get_weight( pd, e, s, err );
89 [ + - ][ - + ]: 437678 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
90 : 437678 : value *= ck;
91 : : }
92 : 9767364 : return true;
93 : : }
94 : :
95 : 0 : bool TMPQualityMetric::evaluate_with_indices( PatchData& pd, size_t p_handle, double& value,
96 : : std::vector< size_t >& indices, MsqError& err )
97 : : {
98 [ # # ]: 0 : indices.resize( MAX_ELEM_NODES );
99 : 0 : size_t num_idx = 0;
100 [ # # ][ # # ]: 0 : bool result = evaluate_internal( pd, p_handle, value, arrptr( indices ), num_idx, err );
101 [ # # ][ # # ]: 0 : if( MSQ_CHKERR( err ) || !result ) return false;
[ # # ][ # # ]
[ # # ][ # # ]
102 : :
103 [ # # ]: 0 : indices.resize( num_idx );
104 : :
105 : : // apply target weight to value
106 [ # # ]: 0 : if( weightCalc )
107 : : {
108 [ # # ]: 0 : const Sample s = ElemSampleQM::sample( p_handle );
109 [ # # ]: 0 : const size_t e = ElemSampleQM::elem( p_handle );
110 [ # # ]: 0 : double ck = weightCalc->get_weight( pd, e, s, err );
111 [ # # ][ # # ]: 0 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ # # ]
112 : 0 : value *= ck;
113 : : }
114 : :
115 : 0 : return true;
116 : : }
117 : :
118 : 96 : static void get_u_perp( const MsqVector< 3 >& u, MsqVector< 3 >& u_perp )
119 : : {
120 : 96 : double a = sqrt( u[0] * u[0] + u[1] * u[1] );
121 [ + + ]: 96 : if( a < 1e-10 )
122 : : {
123 : 4 : u_perp[0] = 1.0;
124 : 4 : u_perp[1] = u_perp[2] = 0.0;
125 : : }
126 : : else
127 : : {
128 : 92 : double b = -u[2] / a;
129 : 92 : u_perp[0] = u[0] * b;
130 : 92 : u_perp[1] = u[1] * b;
131 : 92 : u_perp[2] = a;
132 : : }
133 : 96 : }
134 : :
135 : : /* Do transform M_hat = S_a M_{3x2}, M_{2x2} Theta^-1 M_hat
136 : : * where the plane into which we are projecting is orthogonal
137 : : * to the passed u vector.
138 : : */
139 : 96 : static inline bool project_to_perp_plane( MsqMatrix< 3, 2 > J, const MsqVector< 3 >& u, const MsqVector< 3 >& u_perp,
140 : : MsqMatrix< 2, 2 >& A, MsqMatrix< 3, 2 >& S_a_transpose_Theta )
141 : : {
142 [ + - ][ + - ]: 96 : MsqVector< 3 > n_a = J.column( 0 ) * J.column( 1 );
[ + - ][ + - ]
143 [ + - ]: 96 : double sc, len = length( n_a );
144 [ + - ][ + - ]: 96 : if( !divide( 1.0, len, sc ) ) return false;
145 [ # # ]: 0 : n_a *= sc;
146 [ # # ]: 0 : double ndot = n_a % u;
147 [ # # ]: 0 : double sigma = ( ndot < 0.0 ) ? -1 : 1;
148 : 0 : double cosphi = sigma * ndot;
149 [ # # ][ # # ]: 0 : MsqVector< 3 > cross = n_a * u;
150 [ # # ]: 0 : double sinphi = length( cross );
151 : :
152 [ # # ]: 0 : MsqMatrix< 3, 2 > Theta;
153 [ # # ]: 0 : Theta.set_column( 0, u_perp );
154 [ # # ][ # # ]: 0 : Theta.set_column( 1, u * u_perp );
155 : :
156 : : // If columns of J are not in plane orthogonal to u, then
157 : : // rotate J such that they are.
158 [ # # ]: 0 : if( sinphi > 1e-12 )
159 : : {
160 [ # # ][ # # ]: 0 : MsqVector< 3 > m = sigma * cross;
161 [ # # ][ # # ]: 0 : MsqVector< 3 > n = ( 1 / sinphi ) * m;
162 [ # # ][ # # ]: 0 : MsqVector< 3 > p = ( 1 - cosphi ) * n;
163 [ # # ][ # # ]: 0 : double s_a[] = { p[0] * n[0] + cosphi, p[0] * n[1] - m[2], p[0] * n[2] + m[1],
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
164 [ # # ][ # # ]: 0 : p[1] * n[0] + m[2], p[1] * n[1] + cosphi, p[1] * n[2] - m[0],
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
165 [ # # ][ # # ]: 0 : p[2] * n[0] - m[1], p[2] * n[1] + m[0], p[2] * n[2] + cosphi };
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
166 [ # # ]: 0 : MsqMatrix< 3, 3 > S_a( s_a );
167 [ # # ]: 0 : J = S_a * J;
168 [ # # ][ # # ]: 0 : S_a_transpose_Theta = transpose( S_a ) * Theta;
169 : : }
170 : : else
171 : : {
172 : 0 : S_a_transpose_Theta = Theta;
173 : : // J *= sigma;
174 : : }
175 : :
176 : : // Project to get 2x2 A from A_hat (which might be equal to J)
177 [ # # ][ # # ]: 0 : A = transpose( Theta ) * J;
178 : 96 : return true;
179 : : }
180 : :
181 : : /* Do transform M_hat = S_a M_{3x2}, M_{2x2} Theta^-1 M_hat
182 : : * where the plane into which we are projecting is the cross
183 : : * product of the columns of M, such that S_a is I. Use the
184 : : * first column of M as u_perp.
185 : : *
186 : : * Also pass back the cross product of the columns of M as u,
187 : : * and the first column of M as u_perp, both normalized.
188 : : */
189 : 11987566 : static inline void project_to_matrix_plane( const MsqMatrix< 3, 2 >& M_in, MsqMatrix< 2, 2 >& M_out, MsqVector< 3 >& u,
190 : : MsqVector< 3 >& u_perp )
191 : : {
192 [ + - ][ + - ]: 11987566 : u = M_in.column( 0 ) * M_in.column( 1 );
[ + - ][ + - ]
193 [ + - ][ + - ]: 11987566 : u_perp = M_in.column( 0 );
194 [ + - ]: 11987566 : double len0 = length( u_perp );
195 [ + - ]: 11987566 : double u_len = length( u );
196 : : double d_perp, d_u;
197 [ + - ][ - + ]: 11987566 : if( !divide( 1.0, len0, d_perp ) )
198 : : {
199 : : // try the other column
200 [ # # ][ # # ]: 0 : u_perp = M_in.column( 1 );
201 [ # # ]: 0 : len0 = length( u_perp );
202 [ # # ][ # # ]: 0 : if( !divide( 1.0, len0, d_perp ) )
203 : : {
204 : : // matrix is all zeros
205 [ # # ]: 0 : u[0] = 0;
206 [ # # ]: 0 : u[1] = 0;
207 [ # # ]: 0 : u[2] = 1;
208 [ # # ]: 0 : u_perp[0] = 1;
209 [ # # ]: 0 : u_perp[1] = 0;
210 [ # # ]: 0 : u_perp[2] = 0;
211 [ # # ]: 0 : M_out = MsqMatrix< 2, 2 >( 0.0 );
212 : : }
213 : : else
214 : : {
215 [ # # ]: 0 : MsqMatrix< 3, 2 > junk;
216 [ # # ]: 0 : get_u_perp( u_perp, u );
217 [ # # ]: 0 : project_to_perp_plane( M_in, u, u_perp, M_out, junk );
218 : : }
219 : : }
220 [ + - ][ + + ]: 11987566 : else if( !divide( 1.0, u_len, d_u ) )
221 : : {
222 [ + - ]: 96 : MsqMatrix< 3, 2 > junk;
223 [ + - ]: 96 : get_u_perp( u_perp, u );
224 [ + - ]: 96 : project_to_perp_plane( M_in, u, u_perp, M_out, junk );
225 : : }
226 : : else
227 : : { // the normal case (neither column is zero)
228 [ + - ]: 11987470 : u *= d_u;
229 [ + - ]: 11987470 : u_perp *= d_perp;
230 : :
231 : : // M_out = transpose(theta)*M_in
232 [ + - ]: 11987470 : M_out( 0, 0 ) = len0;
233 [ + - ][ + - ]: 11987470 : M_out( 0, 1 ) = u_perp % M_in.column( 1 );
[ + - ]
234 [ + - ]: 11987470 : M_out( 1, 0 ) = 0.0;
235 [ + - ]: 11987470 : M_out( 1, 1 ) = u_len / len0;
236 : : }
237 : 11987566 : }
238 : :
239 : 11987566 : bool TMPQualityMetric::evaluate_surface_common( PatchData& pd, Sample s, size_t e, const NodeSet& bits, size_t* indices,
240 : : size_t& num_indices, MsqVector< 2 >* derivs, MsqMatrix< 2, 2 >& W,
241 : : MsqMatrix< 2, 2 >& A, MsqMatrix< 3, 2 >& S_a_transpose_Theta,
242 : : MsqError& err )
243 : : {
244 [ + - ][ + - ]: 11987566 : EntityTopology type = pd.element_by_index( e ).get_element_type();
245 : :
246 [ + - ]: 11987566 : const MappingFunction2D* mf = pd.get_mapping_function_2D( type );
247 [ - + ]: 11987566 : if( !mf )
248 : : {
249 [ # # ][ # # ]: 0 : MSQ_SETERR( err )( "No mapping function for element type", MsqError::UNSUPPORTED_ELEMENT );
250 : 0 : return false;
251 : : }
252 : :
253 [ + - ]: 11987566 : MsqMatrix< 3, 2 > J;
254 [ + - ]: 11987566 : mf->jacobian( pd, e, bits, s, indices, derivs, num_indices, J, err );
255 : :
256 : : // If we have a 3x2 target matrix
257 [ + - ][ - + ]: 11987566 : if( targetCalc->have_surface_orient() )
258 : : {
259 [ # # ][ # # ]: 0 : MsqVector< 3 > u, u_perp;
260 [ # # ]: 0 : MsqMatrix< 3, 2 > W_hat;
261 [ # # ]: 0 : targetCalc->get_surface_target( pd, e, s, W_hat, err );
262 [ # # ][ # # ]: 0 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ # # ]
263 : : // Use the cross product of the columns of W as the normal of the
264 : : // plane to work in (i.e. u.). W should have been constructed such
265 : : // that said cross product is in the direction of (n_s)_init. And if
266 : : // for some reason it as not, then using something other than said
267 : : // cross product is likely to produce very wrong results.
268 [ # # ]: 0 : project_to_matrix_plane( W_hat, W, u, u_perp );
269 : : // Do the transforms on A to align it with W and project into the plane.
270 [ # # ][ # # ]: 0 : if( !project_to_perp_plane( J, u, u_perp, A, S_a_transpose_Theta ) ) return false;
271 : : }
272 : : // Otherwise if we have a 2x2 target matrix (i.e. the target does
273 : : // not contain orientation information), project into the plane
274 : : // tangent to J.
275 : : else
276 : : {
277 [ + - ][ + - ]: 11987566 : MsqVector< 3 > u, u_perp;
278 [ + - ]: 11987566 : targetCalc->get_2D_target( pd, e, s, W, err );
279 [ + - ][ - + ]: 11987566 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
280 [ + - ]: 11987566 : project_to_matrix_plane( J, A, u, u_perp );
281 [ + - ]: 11987566 : S_a_transpose_Theta.set_column( 0, u_perp );
282 [ + - ][ + - ]: 11987566 : S_a_transpose_Theta.set_column( 1, u * u_perp );
283 : : // If the domain is set, adjust the sign of things correctly
284 : : // for the case where the element is inverted with respect
285 : : // to the domain.
286 [ + - ][ + + ]: 11987566 : if( pd.domain_set() )
287 : : {
288 [ + - ]: 10776762 : Vector3D n;
289 [ + - ]: 10776762 : pd.get_domain_normal_at_sample( e, s, n, err );
290 [ + - ][ - + ]: 10776762 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
291 : : // if sigma == -1
292 [ + - ][ + - ]: 10776762 : if( Vector3D( u.data() ) % n < 0.0 )
[ + - ][ + + ]
293 : : {
294 : : // flip u
295 [ + - ][ + - ]: 419097 : u = -u;
296 : : // S_a_transpose_Theta == Theta, because S_a == I here.
297 : : // u_perp is unaffected by flipping u, so only the second
298 : : // column of S_a_transpose_Theta and the second row of A
299 : : // are flipped because u x u_perp will be flipped.
300 [ + - ][ + - ]: 419097 : S_a_transpose_Theta.set_column( 1, -S_a_transpose_Theta.column( 1 ) );
[ + - ]
301 [ + - ][ + - ]: 11987566 : A.set_row( 1, -A.row( 1 ) );
[ + - ]
302 : : }
303 : : }
304 : : }
305 : :
306 : 11987566 : return true;
307 : : }
308 : :
309 : 7715396 : void TMPQualityMetric::weight( PatchData& pd, Sample p_sample, size_t p_elem, int num_idx, double& value,
310 : : Vector3D* grad, SymMatrix3D* diag, Matrix3D* hess, MsqError& err )
311 : : {
312 [ + + ]: 7715396 : if( !weightCalc ) return;
313 : :
314 [ - + ][ # # ]: 236360 : double ck = weightCalc->get_weight( pd, p_elem, p_sample, err );MSQ_ERRRTN( err );
[ - + ]
315 : 236360 : value *= ck;
316 [ + - ]: 236360 : if( grad )
317 : : {
318 [ + + ]: 992740 : for( int i = 0; i < num_idx; ++i )
319 : 756380 : grad[i] *= ck;
320 : : }
321 [ - + ]: 236360 : if( diag )
322 : : {
323 [ # # ]: 0 : for( int i = 0; i < num_idx; ++i )
324 : 0 : diag[i] *= ck;
325 : : }
326 [ + - ]: 236360 : if( hess )
327 : : {
328 : 236360 : const int n = num_idx * ( num_idx + 1 ) / 2;
329 [ + + ]: 1967990 : for( int i = 0; i < n; ++i )
330 : 1731630 : hess[i] *= ck;
331 : : }
332 : : }
333 : :
334 : 60 : void TMPQualityMetric::initialize_queue( MeshDomainAssoc* mesh_and_domain, const Settings* settings, MsqError& err )
335 : : {
336 [ - + ][ # # ]: 60 : targetCalc->initialize_queue( mesh_and_domain, settings, err );MSQ_ERRRTN( err );
[ - + ]
337 [ + + ]: 60 : if( weightCalc )
338 : : {
339 [ - + ][ # # ]: 4 : weightCalc->initialize_queue( mesh_and_domain, settings, err );MSQ_ERRRTN( err );
[ - + ]
340 : : }
341 : : }
342 : :
343 [ + - ][ + - ]: 120 : } // namespace MBMesquite
|
