Branch data Line data Source code
1 : : /* *****************************************************************
2 : : MESQUITE -- The Mesh Quality Improvement Toolkit
3 : :
4 : : Copyright 2004 Sandia Corporation and Argonne National
5 : : Laboratory. Under the terms of Contract DE-AC04-94AL85000
6 : : with Sandia Corporation, the U.S. Government retains certain
7 : : rights in 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 : : [email protected], [email protected], [email protected],
24 : : [email protected], [email protected], [email protected]
25 : : [email protected]
26 : :
27 : : ***************************************************************** */
28 : : /*!
29 : : \file QualityMetric.cpp
30 : : \brief
31 : :
32 : : \author Michael Brewer
33 : : \author Thomas Leurent
34 : : \date 2002-05-14
35 : : \author Jason Kraftcheck
36 : : \date 2006-04-20
37 : : */
38 : :
39 : : #include "QualityMetric.hpp"
40 : : #include "MsqVertex.hpp"
41 : : #include "PatchData.hpp"
42 : :
43 : : namespace MBMesquite
44 : : {
45 : :
46 : 209 : void QualityMetric::initialize_queue( MeshDomainAssoc*, const Settings*, MsqError& ) {}
47 : :
48 : 343959 : void QualityMetric::get_single_pass( PatchData& pd, std::vector< size_t >& handles, bool free_vertices_only,
49 : : MsqError& err )
50 : : {
51 : 343959 : get_evaluations( pd, handles, free_vertices_only, err );
52 : 343959 : }
53 : :
54 : 986140 : static inline double get_delta_C( const PatchData& pd, const std::vector< size_t >& indices, MsqError& err )
55 : : {
56 [ - + ]: 986140 : if( indices.empty() )
57 : : {
58 [ # # ][ # # ]: 0 : MSQ_SETERR( err )( MsqError::INVALID_ARG );
59 : 0 : return 0.0;
60 : : }
61 : :
62 : 986140 : std::vector< size_t >::const_iterator beg, iter, iter2, end;
63 : :
64 [ + - ]: 986140 : std::vector< size_t > tmp_vect;
65 [ + + ]: 986140 : if( indices.size() == 1u )
66 : : {
67 [ + - ][ + - ]: 362937 : pd.get_adjacent_vertex_indices( indices.front(), tmp_vect, err );
68 [ + - ][ - + ]: 362937 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
69 [ - + ]: 362937 : assert( !tmp_vect.empty() );
70 [ + - ][ + - ]: 362937 : tmp_vect.push_back( indices.front() );
71 [ + - ]: 362937 : beg = tmp_vect.begin();
72 [ + - ]: 362937 : end = tmp_vect.end();
73 : : }
74 : : else
75 : : {
76 : 623203 : beg = indices.begin();
77 : 623203 : end = indices.end();
78 : : }
79 : :
80 : 986140 : double min_dist_sqr = HUGE_VAL, sum_dist_sqr = 0.0;
81 [ + - ][ + - ]: 7962343 : for( iter = beg; iter != end; ++iter )
[ + + ]
82 : : {
83 [ + - ][ + - ]: 45425677 : for( iter2 = iter + 1; iter2 != end; ++iter2 )
[ + - ][ + + ]
84 : : {
85 [ + - ][ + - ]: 38449474 : Vector3D diff = pd.vertex_by_index( *iter );
[ + - ]
86 [ + - ][ + - ]: 38449474 : diff -= pd.vertex_by_index( *iter2 );
[ + - ]
87 [ + - ]: 38449474 : double dist_sqr = diff % diff;
88 [ + + ]: 38449474 : if( dist_sqr < min_dist_sqr ) min_dist_sqr = dist_sqr;
89 : 38449474 : sum_dist_sqr += dist_sqr;
90 : : }
91 : : }
92 : :
93 [ + - ]: 986140 : return 3e-6 * sqrt( min_dist_sqr ) + 5e-7 * sqrt( sum_dist_sqr / ( end - beg ) );
94 : : }
95 : :
96 : 441653 : bool QualityMetric::evaluate_with_gradient( PatchData& pd, size_t handle, double& value, std::vector< size_t >& indices,
97 : : std::vector< Vector3D >& gradient, MsqError& err )
98 : : {
99 : 441653 : indices.clear();
100 : 441653 : bool valid = evaluate_with_indices( pd, handle, value, indices, err );
101 [ - + ][ # # ]: 441653 : if( MSQ_CHKERR( err ) || !valid ) return false;
[ - + ][ - + ]
102 [ + + ]: 441653 : if( indices.empty() ) return true;
103 : :
104 : : // get initial pertubation amount
105 : 441284 : double delta_C = finiteDiffEps;
106 [ + + ]: 441284 : if( !haveFiniteDiffEps )
107 : : {
108 : 440384 : delta_C = get_delta_C( pd, indices, err );
109 [ - + ][ # # ]: 440384 : MSQ_ERRZERO( err );
[ - + ]
110 [ + + ]: 440384 : if( keepFiniteDiffEps )
111 : : {
112 : 222 : finiteDiffEps = delta_C;
113 : 222 : haveFiniteDiffEps = true;
114 : : }
115 : : }
116 : 441284 : const double delta_inv_C = 1.0 / delta_C;
117 : 441284 : const int reduction_limit = 15;
118 : :
119 : 441284 : gradient.resize( indices.size() );
120 [ + + ]: 1049290 : for( size_t v = 0; v < indices.size(); ++v )
121 : : {
122 [ + - ][ + - ]: 608006 : const Vector3D pos = pd.vertex_by_index( indices[v] );
[ + - ]
123 : :
124 : : /* gradient in the x, y, z direction */
125 [ + + ]: 2432024 : for( int j = 0; j < 3; ++j )
126 : : {
127 : 1824018 : double delta = delta_C;
128 : 1824018 : double delta_inv = delta_inv_C;
129 : : double metric_value;
130 [ + - ]: 1824018 : Vector3D delta_v( 0, 0, 0 );
131 : :
132 : : // perturb the node and calculate gradient. The while loop is a
133 : : // safety net to make sure the epsilon perturbation does not take
134 : : // the element out of the feasible region.
135 : 1824018 : int counter = 0;
136 : : for( ;; )
137 : : {
138 : : // perturb the coordinates of the free vertex in the j direction
139 : : // by delta
140 [ + - ]: 1824018 : delta_v[j] = delta;
141 [ + - ][ + - ]: 1824018 : pd.set_vertex_coordinates( pos + delta_v, indices[v], err );
[ + - ]
142 [ + - ][ - + ]: 1824018 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
143 : :
144 : : // compute the function at the perturbed point location
145 [ + - ]: 1824018 : valid = evaluate( pd, handle, metric_value, err );
146 [ + - ][ - + ]: 1824018 : if( !MSQ_CHKERR( err ) && valid ) break;
[ # # ][ # # ]
[ + - ][ + - ]
147 : :
148 [ # # ]: 0 : if( ++counter >= reduction_limit )
149 : : {
150 : : MSQ_SETERR( err )
151 [ # # ][ # # ]: 0 : ( "Perturbing vertex by delta caused an inverted element.", MsqError::INTERNAL_ERROR );
152 : 0 : return false;
153 : : }
154 : :
155 : 0 : delta *= 0.1;
156 : 0 : delta_inv *= 10.;
157 : : }
158 : : // put the coordinates back where they belong
159 [ + - ][ + - ]: 1824018 : pd.set_vertex_coordinates( pos, indices[v], err );
160 : : // compute the numerical gradient
161 [ + - ][ + - ]: 1824018 : gradient[v][j] = ( metric_value - value ) * delta_inv;
162 : 0 : } // for(j)
163 : : } // for(indices)
164 : 441653 : return true;
165 : : }
166 : :
167 : 546922 : bool QualityMetric::evaluate_with_Hessian( PatchData& pd, size_t handle, double& value, std::vector< size_t >& indices,
168 : : std::vector< Vector3D >& gradient, std::vector< Matrix3D >& Hessian,
169 : : MsqError& err )
170 : : {
171 : 546922 : indices.clear();
172 : 546922 : gradient.clear();
173 : 546922 : keepFiniteDiffEps = true;
174 [ + - ]: 546922 : bool valid = evaluate_with_gradient( pd, handle, value, indices, gradient, err );
175 : 546922 : keepFiniteDiffEps = false;
176 [ + - ][ - + ]: 546922 : if( MSQ_CHKERR( err ) || !valid )
[ # # ][ # # ]
[ - + ][ - + ]
177 : : {
178 : 0 : haveFiniteDiffEps = false;
179 : 0 : return false;
180 : : }
181 [ + + ]: 546922 : if( indices.empty() )
182 : : {
183 : 944 : haveFiniteDiffEps = false;
184 : 944 : return true;
185 : : }
186 : :
187 : : // get initial pertubation amount
188 : : double delta_C;
189 [ + + ]: 545978 : if( haveFiniteDiffEps ) { delta_C = finiteDiffEps; }
190 : : else
191 : : {
192 [ + - ]: 545756 : delta_C = get_delta_C( pd, indices, err );
193 [ + - ][ - + ]: 545756 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
194 : : }
195 [ - + ]: 545978 : assert( delta_C < 1e30 );
196 : 545978 : const double delta_inv_C = 1.0 / delta_C;
197 : 545978 : const int reduction_limit = 15;
198 : :
199 [ + - ]: 545978 : std::vector< Vector3D > temp_gradient( indices.size() );
200 : 545978 : const int num_hess = indices.size() * ( indices.size() + 1 ) / 2;
201 [ + - ]: 545978 : Hessian.resize( num_hess );
202 : :
203 [ + + ]: 2141134 : for( unsigned v = 0; v < indices.size(); ++v )
204 : : {
205 [ + - ][ + - ]: 1595156 : const Vector3D pos = pd.vertex_by_index( indices[v] );
[ + - ]
206 [ + + ]: 6380624 : for( int j = 0; j < 3; ++j )
207 : : { // x, y, and z
208 : 4785468 : double delta = delta_C;
209 : 4785468 : double delta_inv = delta_inv_C;
210 : : double metric_value;
211 [ + - ]: 4785468 : Vector3D delta_v( 0, 0, 0 );
212 : :
213 : : // find finite difference for gradient
214 : 4785468 : int counter = 0;
215 : : for( ;; )
216 : : {
217 [ + - ]: 4785468 : delta_v[j] = delta;
218 [ + - ][ + - ]: 4785468 : pd.set_vertex_coordinates( pos + delta_v, indices[v], err );
[ + - ]
219 [ + - ][ - + ]: 4785468 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
220 [ + - ]: 4785468 : valid = evaluate_with_gradient( pd, handle, metric_value, indices, temp_gradient, err );
221 [ + - ][ - + ]: 4785468 : if( !MSQ_CHKERR( err ) && valid ) break;
[ # # ][ # # ]
[ + - ][ + - ]
222 : :
223 [ # # ]: 0 : if( ++counter >= reduction_limit )
224 : : {
225 : : MSQ_SETERR( err )
226 : : ( "Algorithm did not successfully compute element's "
227 : : "Hessian.\n",
228 [ # # ][ # # ]: 0 : MsqError::INTERNAL_ERROR );
229 : 0 : haveFiniteDiffEps = false;
230 : 0 : return false;
231 : : }
232 : :
233 : 0 : delta *= 0.1;
234 : 0 : delta_inv *= 10.0;
235 : : }
236 [ + - ][ + - ]: 4785468 : pd.set_vertex_coordinates( pos, indices[v], err );
237 [ + - ][ - + ]: 4785468 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
238 : :
239 : : // compute the numerical Hessian
240 [ + + ]: 14273523 : for( unsigned w = 0; w <= v; ++w )
241 : : {
242 : : // finite difference to get some entries of the Hessian
243 [ + - ][ + - ]: 9488055 : Vector3D fd( temp_gradient[w] );
244 [ + - ][ + - ]: 9488055 : fd -= gradient[w];
245 [ + - ]: 9488055 : fd *= delta_inv;
246 : : // For the block at position w,v in a matrix, we need the corresponding index
247 : : // (mat_index) in a 1D array containing only upper triangular blocks.
248 : 9488055 : unsigned sum_w = w * ( w + 1 ) / 2; // 1+2+3+...+w
249 : 9488055 : unsigned mat_index = w * indices.size() + v - sum_w;
250 [ + - ][ + - ]: 9488055 : Hessian[mat_index][0][j] = fd[0];
[ + - ]
251 [ + - ][ + - ]: 9488055 : Hessian[mat_index][1][j] = fd[1];
[ + - ]
252 [ + - ][ + - ]: 9488055 : Hessian[mat_index][2][j] = fd[2];
[ + - ]
253 : : }
254 : 0 : } // for(j)
255 : : } // for(indices)
256 : 545978 : haveFiniteDiffEps = false;
257 : 546922 : return true;
258 : : }
259 : :
260 : 0 : bool QualityMetric::evaluate_with_Hessian_diagonal( PatchData& pd, size_t handle, double& value,
261 : : std::vector< size_t >& indices, std::vector< Vector3D >& gradient,
262 : : std::vector< SymMatrix3D >& Hessian_diagonal, MsqError& err )
263 : : {
264 [ # # ]: 0 : bool rval = evaluate_with_Hessian( pd, handle, value, indices, gradient, tmpHess, err );
265 [ # # ][ # # ]: 0 : if( MSQ_CHKERR( err ) || !rval ) return rval;
[ # # ][ # # ]
[ # # ][ # # ]
266 : 0 : size_t s = indices.size();
267 [ # # ]: 0 : Hessian_diagonal.resize( s );
268 [ # # ]: 0 : std::vector< Matrix3D >::const_iterator h = tmpHess.begin();
269 [ # # ]: 0 : for( size_t i = 0; i < indices.size(); ++i )
270 : : {
271 [ # # ][ # # ]: 0 : Hessian_diagonal[i] = h->upper();
[ # # ]
272 [ # # ]: 0 : h += s--;
273 : : }
274 : 0 : return rval;
275 : : }
276 : :
277 : 292591 : uint32_t QualityMetric::fixed_vertex_bitmap( PatchData& pd, const MsqMeshEntity* elem, std::vector< size_t >& indices )
278 : : {
279 : 292591 : indices.clear();
280 : 292591 : uint32_t result = ~(uint32_t)0;
281 : 292591 : unsigned num_vtx = elem->vertex_count();
282 : 292591 : const size_t* vertices = elem->get_vertex_index_array();
283 : 292591 : indices.clear();
284 [ + + ]: 1571814 : for( unsigned i = 0; i < num_vtx; ++i )
285 : : {
286 [ + + ]: 1279223 : if( vertices[i] < pd.num_free_vertices() )
287 : : {
288 : 865973 : indices.push_back( vertices[i] );
289 : 865973 : result &= ~( uint32_t )( 1 << i );
290 : : }
291 : : }
292 : 292591 : return result;
293 : : }
294 : :
295 : 292588 : void QualityMetric::remove_fixed_gradients( EntityTopology elem_type, uint32_t fixed, std::vector< Vector3D >& grads )
296 : : {
297 : 292588 : const unsigned num_vertex = TopologyInfo::corners( elem_type );
298 : : unsigned r, w;
299 [ + + ][ + + ]: 912233 : for( r = 0; r < num_vertex && !( fixed & ( 1 << r ) ); ++r )
300 : : ;
301 [ + + ]: 766921 : for( w = r++; r < num_vertex; ++r )
302 : : {
303 [ + + ]: 474333 : if( !( fixed & ( 1 << r ) ) )
304 : : {
305 : 246316 : grads[w] = grads[r];
306 : 246316 : ++w;
307 : : }
308 : : }
309 : 292588 : grads.resize( w );
310 : 292588 : }
311 : :
312 : 0 : void QualityMetric::remove_fixed_diagonals( EntityTopology type, uint32_t fixed, std::vector< Vector3D >& grads,
313 : : std::vector< SymMatrix3D >& diags )
314 : : {
315 : 0 : const unsigned num_vertex = TopologyInfo::corners( type );
316 : : unsigned r, w;
317 [ # # ][ # # ]: 0 : for( r = 0; r < num_vertex && !( fixed & ( 1 << r ) ); ++r )
318 : : ;
319 [ # # ]: 0 : for( w = r++; r < num_vertex; ++r )
320 : : {
321 [ # # ]: 0 : if( !( fixed & ( 1 << r ) ) )
322 : : {
323 : 0 : grads[w] = grads[r];
324 : 0 : diags[w] = diags[r];
325 : 0 : ++w;
326 : : }
327 : : }
328 : 0 : grads.resize( w );
329 : 0 : diags.resize( w );
330 : 0 : }
331 : :
332 : 139808 : void QualityMetric::remove_fixed_hessians( EntityTopology elem_type, uint32_t fixed, std::vector< Matrix3D >& hessians )
333 : : {
334 : 139808 : const unsigned num_vertex = TopologyInfo::corners( elem_type );
335 : 139808 : unsigned r, c, i = 0, w = 0;
336 [ + + ]: 751030 : for( r = 0; r < num_vertex; ++r )
337 : : {
338 [ + + ]: 611222 : if( fixed & ( 1 << r ) )
339 : : {
340 : 197674 : i += num_vertex - r;
341 : 197674 : continue;
342 : : }
343 [ + + ]: 1562124 : for( c = r; c < num_vertex; ++c )
344 : : {
345 [ + + ]: 1148576 : if( !( fixed & ( 1 << c ) ) )
346 : : {
347 : 977442 : hessians[w] = hessians[i];
348 : 977442 : ++w;
349 : : }
350 : 1148576 : ++i;
351 : : }
352 : : }
353 : 139808 : hessians.resize( w );
354 : 139808 : }
355 : :
356 : 0 : double QualityMetric::weighted_average_metrics( const double coef[], const double metric_values[],
357 : : const int& num_values, MsqError& /*err*/ )
358 : : {
359 : : // MSQ_MAX needs to be made global?
360 : : // double MSQ_MAX=1e10;
361 : 0 : double total_value = 0.0;
362 : 0 : int i = 0;
363 : : // if no values, return zero
364 [ # # ]: 0 : if( num_values <= 0 ) { return 0.0; }
365 : :
366 [ # # ]: 0 : for( i = 0; i < num_values; ++i )
367 : : {
368 : 0 : total_value += coef[i] * metric_values[i];
369 : : }
370 : 0 : total_value /= (double)num_values;
371 : :
372 : 0 : return total_value;
373 : : }
374 : :
375 [ + - ][ + - ]: 120 : } // namespace MBMesquite
|
