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 TQualityMetric.cpp
28 : : * \brief
29 : : * \author Jason Kraftcheck
30 : : */
31 : :
32 : : #undef PRINT_INFO
33 : :
34 : : #include "Mesquite.hpp"
35 : : #include "TQualityMetric.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 "TMetric.hpp"
45 : : #include "TMetricBarrier.hpp"
46 : : #include "TargetMetricUtil.hpp"
47 : : #include "TMPDerivs.hpp"
48 : :
49 : : #ifdef PRINT_INFO
50 : : #include <iostream>
51 : : #endif
52 : :
53 : : #include <functional>
54 : : #include <algorithm>
55 : :
56 : : #define NUMERICAL_2D_HESSIAN
57 : :
58 : : namespace MBMesquite
59 : : {
60 : :
61 : 46 : std::string TQualityMetric::get_name() const
62 : : {
63 : 46 : return targetMetric->get_name();
64 : : }
65 : :
66 : 9767364 : bool TQualityMetric::evaluate_internal( PatchData& pd, size_t p_handle, double& value, size_t* indices,
67 : : size_t& num_indices, MsqError& err )
68 : : {
69 [ + - ]: 9767364 : const Sample s = ElemSampleQM::sample( p_handle );
70 [ + - ]: 9767364 : const size_t e = ElemSampleQM::elem( p_handle );
71 [ + - ]: 9767364 : MsqMeshEntity& p_elem = pd.element_by_index( e );
72 [ + - ]: 9767364 : EntityTopology type = p_elem.get_element_type();
73 [ + - ]: 9767364 : unsigned edim = TopologyInfo::dimension( type );
74 [ + - ]: 9767364 : const NodeSet bits = pd.non_slave_node_set( e );
75 : :
76 : : bool rval;
77 [ + + ]: 9767364 : if( edim == 3 )
78 : : { // 3x3 or 3x2 targets ?
79 [ + - ]: 4739254 : const MappingFunction3D* mf = pd.get_mapping_function_3D( type );
80 [ - + ]: 4739254 : if( !mf )
81 : : {
82 : : MSQ_SETERR( err )
83 [ # # ][ # # ]: 0 : ( "No mapping function for element type", MsqError::UNSUPPORTED_ELEMENT );
84 : 36 : return false;
85 : : }
86 : :
87 [ + - ][ + - ]: 4739254 : MsqMatrix< 3, 3 > A, W;
88 [ + - ]: 4739254 : mf->jacobian( pd, e, bits, s, indices, mDerivs3D, num_indices, A, err );
89 [ + - ][ - + ]: 4739254 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
90 [ + - ]: 4739254 : targetCalc->get_3D_target( pd, e, s, W, err );
91 [ + - ][ - + ]: 4739254 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
92 [ + - ]: 4739254 : const MsqMatrix< 3, 3 > Winv = inverse( W );
93 [ + - ]: 4739254 : const MsqMatrix< 3, 3 > T = A * Winv;
94 [ + - ]: 4739254 : rval = targetMetric->evaluate( T, value, err );
95 [ + - ][ + + ]: 4739254 : MSQ_ERRZERO( err );
[ + - ][ + - ]
[ + + ]
96 : : #ifdef PRINT_INFO
97 : : print_info< 3 >( e, s, A, W, A * inverse( W ) );
98 : : #endif
99 : : }
100 [ + - ]: 5028110 : else if( edim == 2 )
101 : : {
102 [ + - ][ + - ]: 5028110 : MsqMatrix< 2, 2 > W, A;
103 [ + - ]: 5028110 : MsqMatrix< 3, 2 > S_a_transpose_Theta;
104 : : rval =
105 [ + - ]: 5028110 : evaluate_surface_common( pd, s, e, bits, indices, num_indices, mDerivs2D, W, A, S_a_transpose_Theta, err );
106 [ + - ][ - + ]: 5029293 : if( MSQ_CHKERR( err ) || !rval ) return false;
[ # # ][ # # ]
[ - + ][ - + ]
107 [ + - ]: 5028110 : const MsqMatrix< 2, 2 > Winv = inverse( W );
108 [ + - ]: 5028110 : const MsqMatrix< 2, 2 > T = A * Winv;
109 [ + - ]: 5028110 : rval = targetMetric->evaluate( T, value, err );
110 [ + - ][ + + ]: 5028110 : MSQ_ERRZERO( err );
[ + - ][ + - ]
[ + + ]
111 : :
112 : : #ifdef PRINT_INFO
113 : : print_info< 2 >( e, s, J, Wp, A * inverse( W ) );
114 : : #endif
115 : : }
116 : : else
117 : : {
118 : 0 : assert( false );
119 : : return false;
120 : : }
121 : :
122 : 9767364 : return rval;
123 : : }
124 : :
125 : 7479056 : bool TQualityMetric::evaluate_with_gradient( PatchData& pd, size_t p_handle, double& value,
126 : : std::vector< size_t >& indices, std::vector< Vector3D >& grad,
127 : : MsqError& err )
128 : : {
129 [ + - ]: 7479056 : const Sample s = ElemSampleQM::sample( p_handle );
130 [ + - ]: 7479056 : const size_t e = ElemSampleQM::elem( p_handle );
131 [ + - ]: 7479056 : MsqMeshEntity& p_elem = pd.element_by_index( e );
132 [ + - ]: 7479056 : EntityTopology type = p_elem.get_element_type();
133 [ + - ]: 7479056 : unsigned edim = TopologyInfo::dimension( type );
134 : 7479056 : size_t num_idx = 0;
135 [ + - ]: 7479056 : const NodeSet bits = pd.non_slave_node_set( e );
136 : :
137 : : bool rval;
138 [ + + ]: 7479056 : if( edim == 3 )
139 : : { // 3x3 or 3x2 targets ?
140 [ + - ]: 519600 : const MappingFunction3D* mf = pd.get_mapping_function_3D( type );
141 [ - + ]: 519600 : if( !mf )
142 : : {
143 : : MSQ_SETERR( err )
144 [ # # ][ # # ]: 0 : ( "No mapping function for element type", MsqError::UNSUPPORTED_ELEMENT );
145 : 0 : return false;
146 : : }
147 : :
148 [ + - ][ + - ]: 519600 : MsqMatrix< 3, 3 > A, W, dmdT;
[ + - ]
149 [ + - ]: 519600 : mf->jacobian( pd, e, bits, s, mIndices, mDerivs3D, num_idx, A, err );
150 [ + - ][ - + ]: 519600 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
151 [ + - ]: 519600 : targetCalc->get_3D_target( pd, e, s, W, err );
152 [ + - ][ - + ]: 519600 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
153 [ + - ]: 519600 : const MsqMatrix< 3, 3 > Winv = inverse( W );
154 [ + - ]: 519600 : const MsqMatrix< 3, 3 > T = A * Winv;
155 [ + - ]: 519600 : rval = targetMetric->evaluate_with_grad( T, value, dmdT, err );
156 [ + - ][ - + ]: 519600 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
157 [ + - ][ + - ]: 519600 : gradient< 3 >( num_idx, mDerivs3D, dmdT * transpose( Winv ), grad );
[ + - ]
158 : : #ifdef PRINT_INFO
159 : : print_info< 3 >( e, s, A, W, A * inverse( W ) );
160 : : #endif
161 : : }
162 [ + - ]: 6959456 : else if( edim == 2 )
163 : : {
164 [ + - ][ + - ]: 6959456 : MsqMatrix< 2, 2 > W, A, dmdT;
[ + - ]
165 [ + - ]: 6959456 : MsqMatrix< 3, 2 > S_a_transpose_Theta;
166 [ + - ]: 6959456 : rval = evaluate_surface_common( pd, s, e, bits, mIndices, num_idx, mDerivs2D, W, A, S_a_transpose_Theta, err );
167 [ + - ][ - + ]: 6959476 : if( MSQ_CHKERR( err ) || !rval ) return false;
[ # # ][ # # ]
[ - + ][ - + ]
168 [ + - ]: 6959456 : const MsqMatrix< 2, 2 > Winv = inverse( W );
169 [ + - ]: 6959456 : const MsqMatrix< 2, 2 > T = A * Winv;
170 [ + - ]: 6959456 : rval = targetMetric->evaluate_with_grad( T, value, dmdT, err );
171 [ + - ][ + + ]: 6959456 : MSQ_ERRZERO( err );
[ + - ][ + - ]
[ + + ]
172 [ + - ][ + - ]: 6959436 : gradient< 2 >( num_idx, mDerivs2D, S_a_transpose_Theta * dmdT * transpose( Winv ), grad );
[ + - ][ + - ]
173 : : #ifdef PRINT_INFO
174 : : print_info< 2 >( e, s, J, Wp, A * inverse( W ) );
175 : : #endif
176 : : }
177 : : else
178 : : {
179 : 0 : assert( false );
180 : : return false;
181 : : }
182 : :
183 : : // pass back index list
184 [ + - ]: 7479036 : indices.resize( num_idx );
185 [ + - ]: 7479036 : std::copy( mIndices, mIndices + num_idx, indices.begin() );
186 : :
187 : : // apply target weight to value
188 [ + + ][ + - ]: 7479036 : weight( pd, s, e, num_idx, value, grad.empty() ? 0 : arrptr( grad ), 0, 0, err );
[ + - ]
189 [ + - ][ - + ]: 7479036 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
190 : 7479056 : return rval;
191 : : }
192 : :
193 : 783060 : bool TQualityMetric::evaluate_with_Hessian( PatchData& pd, size_t p_handle, double& value,
194 : : std::vector< size_t >& indices, std::vector< Vector3D >& grad,
195 : : std::vector< Matrix3D >& Hessian, MsqError& err )
196 : : {
197 [ + - ]: 783060 : const Sample s = ElemSampleQM::sample( p_handle );
198 [ + - ]: 783060 : const size_t e = ElemSampleQM::elem( p_handle );
199 [ + - ]: 783060 : MsqMeshEntity& p_elem = pd.element_by_index( e );
200 [ + - ]: 783060 : EntityTopology type = p_elem.get_element_type();
201 [ + - ]: 783060 : unsigned edim = TopologyInfo::dimension( type );
202 : 783060 : size_t num_idx = 0;
203 [ + - ]: 783060 : const NodeSet bits = pd.non_slave_node_set( e );
204 : :
205 : : bool rval;
206 [ + + ]: 783060 : if( edim == 3 )
207 : : { // 3x3 or 3x2 targets ?
208 [ + - ]: 236360 : const MappingFunction3D* mf = pd.get_mapping_function_3D( type );
209 [ - + ]: 236360 : if( !mf )
210 : : {
211 : : MSQ_SETERR( err )
212 [ # # ][ # # ]: 0 : ( "No mapping function for element type", MsqError::UNSUPPORTED_ELEMENT );
213 : 0 : return false;
214 : : }
215 : :
216 [ + - ][ + - ]: 1654520 : MsqMatrix< 3, 3 > A, W, dmdT, d2mdT2[6];
[ + - ][ + - ]
[ + + ]
217 [ + - ]: 236360 : mf->jacobian( pd, e, bits, s, mIndices, mDerivs3D, num_idx, A, err );
218 [ + - ][ - + ]: 236360 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
219 [ + - ]: 236360 : targetCalc->get_3D_target( pd, e, s, W, err );
220 [ + - ][ - + ]: 236360 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
221 [ + - ]: 236360 : const MsqMatrix< 3, 3 > Winv = inverse( W );
222 [ + - ]: 236360 : const MsqMatrix< 3, 3 > T = A * Winv;
223 [ + - ]: 236360 : rval = targetMetric->evaluate_with_hess( T, value, dmdT, d2mdT2, err );
224 [ + - ][ - + ]: 236360 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
225 [ + - ][ + - ]: 236360 : gradient< 3 >( num_idx, mDerivs3D, dmdT * transpose( Winv ), grad );
[ + - ]
226 [ + - ]: 236360 : second_deriv_wrt_product_factor( d2mdT2, Winv );
227 [ + - ]: 236360 : Hessian.resize( num_idx * ( num_idx + 1 ) / 2 );
228 [ + - ][ + - ]: 236360 : if( num_idx ) hessian< 3 >( num_idx, mDerivs3D, d2mdT2, arrptr( Hessian ) );
[ + - ]
229 : :
230 : : #ifdef PRINT_INFO
231 : : print_info< 3 >( e, s, A, W, A * inverse( W ) );
232 : : #endif
233 : : }
234 [ + - ]: 546700 : else if( edim == 2 )
235 : : {
236 : : #ifdef NUMERICAL_2D_HESSIAN
237 : : // return finite difference approximation for now
238 : :
239 [ + - ]: 546700 : return QualityMetric::evaluate_with_Hessian( pd, p_handle, value, indices, grad, Hessian, err );
240 : : #else
241 : : MsqMatrix< 2, 2 > W, A, dmdT, d2mdT2[3];
242 : : MsqMatrix< 3, 2 > M;
243 : : rval = evaluate_surface_common( pd, s, e, bits, mIndices, num_idx, mDerivs2D, W, A, M, err );
244 : : if( MSQ_CHKERR( err ) || !rval ) return false;
245 : : const MsqMatrix< 2, 2 > Winv = inverse( W );
246 : : const MsqMatrix< 2, 2 > T = A * Winv;
247 : : rval = targetMetric->evaluate_with_hess( T, value, dmdT, d2mdT2, err );
248 : : MSQ_ERRZERO( err );
249 : : gradient< 2 >( num_idx, mDerivs2D, M * dmdT * transpose( Winv ), grad );
250 : : // calculate 2D hessian
251 : : second_deriv_wrt_product_factor( d2mdT2, Winv );
252 : : const size_t n = num_idx * ( num_idx + 1 ) / 2;
253 : : hess2d.resize( n );
254 : : if( n ) hessian< 2 >( num_idx, mDerivs2D, d2mdT2, arrptr( hess2d ) );
255 : : // calculate surface hessian as transform of 2D hessian
256 : : Hessian.resize( n );
257 : : for( size_t i = 0; i < n; ++i )
258 : : Hessian[i] = Matrix3D( ( M * hess2d[i] * transpose( M ) ).data() );
259 : : #ifdef PRINT_INFO
260 : : print_info< 2 >( e, s, J, Wp, A * inverse( W ) );
261 : : #endif
262 : : #endif
263 : : }
264 : : else
265 : : {
266 : 0 : assert( 0 );
267 : : return false;
268 : : }
269 : :
270 : : // pass back index list
271 [ + - ]: 236360 : indices.resize( num_idx );
272 [ + - ]: 236360 : std::copy( mIndices, mIndices + num_idx, indices.begin() );
273 : :
274 : : // apply target weight to value
275 [ - + ]: 236360 : if( !num_idx )
276 [ # # ]: 0 : weight( pd, s, e, num_idx, value, 0, 0, 0, err );
277 : : else
278 [ + - ][ + - ]: 236360 : weight( pd, s, e, num_idx, value, arrptr( grad ), 0, arrptr( Hessian ), err );
[ + - ]
279 [ + - ][ - + ]: 236360 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
280 : 783060 : return rval;
281 : : }
282 : :
283 : 0 : bool TQualityMetric::evaluate_with_Hessian_diagonal( PatchData& pd, size_t p_handle, double& value,
284 : : std::vector< size_t >& indices, std::vector< Vector3D >& grad,
285 : : std::vector< SymMatrix3D >& diagonal, MsqError& err )
286 : : {
287 [ # # ]: 0 : const Sample s = ElemSampleQM::sample( p_handle );
288 [ # # ]: 0 : const size_t e = ElemSampleQM::elem( p_handle );
289 [ # # ]: 0 : MsqMeshEntity& p_elem = pd.element_by_index( e );
290 [ # # ]: 0 : EntityTopology type = p_elem.get_element_type();
291 [ # # ]: 0 : unsigned edim = TopologyInfo::dimension( type );
292 : 0 : size_t num_idx = 0;
293 [ # # ]: 0 : const NodeSet bits = pd.non_slave_node_set( e );
294 : :
295 : : bool rval;
296 [ # # ]: 0 : if( edim == 3 )
297 : : { // 3x3 or 3x2 targets ?
298 [ # # ]: 0 : const MappingFunction3D* mf = pd.get_mapping_function_3D( type );
299 [ # # ]: 0 : if( !mf )
300 : : {
301 : : MSQ_SETERR( err )
302 [ # # ][ # # ]: 0 : ( "No mapping function for element type", MsqError::UNSUPPORTED_ELEMENT );
303 : 0 : return false;
304 : : }
305 : :
306 [ # # ][ # # ]: 0 : MsqMatrix< 3, 3 > A, W, dmdT, d2mdT2[6];
[ # # ][ # # ]
[ # # ]
307 [ # # ]: 0 : mf->jacobian( pd, e, bits, s, mIndices, mDerivs3D, num_idx, A, err );
308 [ # # ][ # # ]: 0 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ # # ]
309 [ # # ]: 0 : targetCalc->get_3D_target( pd, e, s, W, err );
310 [ # # ][ # # ]: 0 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ # # ]
311 [ # # ]: 0 : const MsqMatrix< 3, 3 > Winv = inverse( W );
312 [ # # ]: 0 : const MsqMatrix< 3, 3 > T = A * Winv;
313 [ # # ]: 0 : rval = targetMetric->evaluate_with_hess( T, value, dmdT, d2mdT2, err );
314 [ # # ][ # # ]: 0 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ # # ]
315 [ # # ][ # # ]: 0 : gradient< 3 >( num_idx, mDerivs3D, dmdT * transpose( Winv ), grad );
[ # # ]
316 [ # # ]: 0 : second_deriv_wrt_product_factor( d2mdT2, Winv );
317 : :
318 [ # # ]: 0 : diagonal.resize( num_idx );
319 [ # # ][ # # ]: 0 : hessian_diagonal< 3 >( num_idx, mDerivs3D, d2mdT2, arrptr( diagonal ) );
320 : : #ifdef PRINT_INFO
321 : : print_info< 3 >( e, s, A, W, A * inverse( W ) );
322 : : #endif
323 : : }
324 [ # # ]: 0 : else if( edim == 2 )
325 : : {
326 : : #ifdef NUMERICAL_2D_HESSIAN
327 : : // use finite diference approximation for now
328 [ # # ]: 0 : return QualityMetric::evaluate_with_Hessian_diagonal( pd, p_handle, value, indices, grad, diagonal, err );
329 : : #else
330 : : MsqMatrix< 2, 2 > W, A, dmdT, d2mdT2[3];
331 : : MsqMatrix< 3, 2 > M;
332 : : rval = evaluate_surface_common( pd, s, e, bits, mIndices, num_idx, mDerivs2D, W, A, M, err );
333 : : if( MSQ_CHKERR( err ) || !rval ) return false;
334 : : const MsqMatrix< 2, 2 > Winv = inverse( W );
335 : : const MsqMatrix< 2, 2 > T = A * Winv;
336 : : rval = targetMetric->evaluate_with_hess( T, value, dmdT, d2mdT2, err );
337 : : MSQ_ERRZERO( err );
338 : : gradient< 2 >( num_idx, mDerivs2D, M * dmdT * transpose( Winv ), grad );
339 : : second_deriv_wrt_product_factor( d2mdT2, Winv );
340 : :
341 : : diagonal.resize( num_idx );
342 : : for( size_t i = 0; i < num_idx; ++i )
343 : : {
344 : : MsqMatrix< 2, 2 > block2d;
345 : : block2d( 0, 0 ) = transpose( mDerivs2D[i] ) * d2mdT2[0] * mDerivs2D[i];
346 : : block2d( 0, 1 ) = transpose( mDerivs2D[i] ) * d2mdT2[1] * mDerivs2D[i];
347 : : block2d( 1, 0 ) = block2d( 0, 1 );
348 : : block2d( 1, 1 ) = transpose( mDerivs2D[i] ) * d2mdT2[2] * mDerivs2D[i];
349 : : MsqMatrix< 3, 2 > p = M * block2d;
350 : :
351 : : SymMatrix3D& H = diagonal[i];
352 : : H[0] = p.row( 0 ) * transpose( M.row( 0 ) );
353 : : H[1] = p.row( 0 ) * transpose( M.row( 1 ) );
354 : : H[2] = p.row( 0 ) * transpose( M.row( 2 ) );
355 : : H[3] = p.row( 1 ) * transpose( M.row( 1 ) );
356 : : H[4] = p.row( 1 ) * transpose( M.row( 2 ) );
357 : : H[5] = p.row( 2 ) * transpose( M.row( 2 ) );
358 : : }
359 : : #ifdef PRINT_INFO
360 : : print_info< 2 >( e, s, J, Wp, A * inverse( W ) );
361 : : #endif
362 : : #endif
363 : : }
364 : : else
365 : : {
366 : 0 : assert( 0 );
367 : : return false;
368 : : }
369 : :
370 : : // pass back index list
371 [ # # ]: 0 : indices.resize( num_idx );
372 [ # # ]: 0 : std::copy( mIndices, mIndices + num_idx, indices.begin() );
373 : :
374 : : // apply target weight to value
375 [ # # ]: 0 : if( !num_idx )
376 [ # # ]: 0 : weight( pd, s, e, num_idx, value, 0, 0, 0, err );
377 : : else
378 [ # # ][ # # ]: 0 : weight( pd, s, e, num_idx, value, arrptr( grad ), arrptr( diagonal ), 0, err );
[ # # ]
379 [ # # ][ # # ]: 0 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ # # ]
380 : 0 : return rval;
381 : : }
382 : :
383 [ + - ][ + - ]: 120 : } // namespace MBMesquite
|
