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 : :
26 : : ***************************************************************** */
27 : : /*!
28 : : \file IdealWeightInverseMeanRatio.cpp
29 : : \brief
30 : :
31 : : \author Michael Brewer
32 : : \author Todd Munson
33 : : \author Thomas Leurent
34 : :
35 : : \date 2002-06-9
36 : : */
37 : : #include "IdealWeightInverseMeanRatio.hpp"
38 : : #include "MeanRatioFunctions.hpp"
39 : : #include "MsqTimer.hpp"
40 : : #include "MsqDebug.hpp"
41 : : #include "MsqError.hpp"
42 : : #include "PatchData.hpp"
43 : :
44 : : #include <vector>
45 : : using std::vector;
46 : : #include <math.h>
47 : :
48 : : namespace MBMesquite
49 : : {
50 : :
51 : 9 : IdealWeightInverseMeanRatio::IdealWeightInverseMeanRatio( MsqError& err, double pow_dbl )
52 [ + - ][ + - ]: 1053 : : AveragingQM( QualityMetric::LINEAR )
[ + + ][ + - ]
[ + + ][ + - ]
[ + + ][ + - ]
[ + - ][ + - ]
[ + - # #
# # # # ]
53 : : {
54 [ + - ][ + - ]: 9 : set_metric_power( pow_dbl, err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
[ # # ]
55 : : }
56 : :
57 [ + - ][ + - ]: 1062 : IdealWeightInverseMeanRatio::IdealWeightInverseMeanRatio() : AveragingQM( QualityMetric::LINEAR )
[ + + ][ + - ]
[ + + ][ + - ]
[ + + ][ + - ]
[ + - ][ + - ]
[ + - # #
# # # # ]
58 : : {
59 [ + - ]: 9 : MsqError err;
60 [ + - ]: 9 : set_metric_power( 1.0, err );
61 [ + - ][ - + ]: 9 : assert( !err );
62 [ # # ]: 9 : }
63 : :
64 : 28 : std::string IdealWeightInverseMeanRatio::get_name() const
65 : : {
66 [ + - ]: 28 : return "Inverse Mean Ratio";
67 : : }
68 : :
69 : 152739 : int IdealWeightInverseMeanRatio::get_negate_flag() const
70 : : {
71 [ - + ]: 152739 : return b2Con.value() < 0 ? -1 : 1;
72 : : }
73 : :
74 : : //! Sets the power value in the metric computation.
75 : 18 : void IdealWeightInverseMeanRatio::set_metric_power( double pow_dbl, MsqError& err )
76 : : {
77 [ - + ]: 18 : if( fabs( pow_dbl ) <= MSQ_MIN )
78 : : {
79 [ # # ]: 0 : MSQ_SETERR( err )( MsqError::INVALID_ARG );
80 : 18 : return;
81 : : }
82 : 18 : a2Con = pow( .5, pow_dbl );
83 : 18 : b2Con = pow_dbl;
84 : 18 : c2Con = -pow_dbl;
85 : 18 : a3Con = pow( 1.0 / 3.0, pow_dbl );
86 : 18 : b3Con = pow_dbl;
87 : 18 : c3Con = -2.0 * pow_dbl / 3.0;
88 : : }
89 : :
90 : 84464 : bool IdealWeightInverseMeanRatio::evaluate( PatchData& pd, size_t handle, double& m, MsqError& err )
91 : : {
92 [ + - ]: 84464 : const MsqMeshEntity* e = &pd.element_by_index( handle );
93 [ + - ]: 84464 : EntityTopology topo = e->get_element_type();
94 : :
95 [ + - ]: 84464 : const MsqVertex* vertices = pd.get_vertex_array( err );
96 [ + - ][ - + ]: 84464 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
97 [ + - ]: 84464 : const size_t* v_i = e->get_vertex_index_array();
98 : :
99 [ + - ]: 84464 : Vector3D n; // Surface normal for 2D objects
100 : :
101 : : // Prism and Hex element descriptions
102 : : static const int locs_pri[6][4] = { { 0, 1, 2, 3 }, { 1, 2, 0, 4 }, { 2, 0, 1, 5 },
103 : : { 3, 5, 4, 0 }, { 4, 3, 5, 1 }, { 5, 4, 3, 2 } };
104 : : static const int locs_hex[8][4] = { { 0, 1, 3, 4 }, { 1, 2, 0, 5 }, { 2, 3, 1, 6 }, { 3, 0, 2, 7 },
105 : : { 4, 7, 5, 0 }, { 5, 4, 6, 1 }, { 6, 5, 7, 2 }, { 7, 6, 4, 3 } };
106 : :
107 [ + - ]: 84464 : const Vector3D d_con( 1.0, 1.0, 1.0 );
108 : :
109 : : int i;
110 : :
111 : 84464 : m = 0.0;
112 : 84464 : bool metric_valid = false;
113 [ + + + + : 84464 : switch( topo )
- + - ]
114 : : {
115 : : case TRIANGLE:
116 [ + - ]: 367 : pd.get_domain_normal_at_element( e, n, err );
117 [ + - ][ - + ]: 367 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
118 [ + - ][ + - ]: 367 : n = n / n.length(); // Need unit normal
[ + - ]
119 [ + - ]: 367 : mCoords[0] = vertices[v_i[0]];
120 [ + - ]: 367 : mCoords[1] = vertices[v_i[1]];
121 [ + - ]: 367 : mCoords[2] = vertices[v_i[2]];
122 [ + - ]: 367 : metric_valid = m_fcn_2e( m, mCoords, n, a2Con, b2Con, c2Con );
123 [ + + ]: 367 : if( !metric_valid ) return false;
124 : 364 : break;
125 : :
126 : : case QUADRILATERAL:
127 [ + - ]: 17784 : pd.get_domain_normal_at_element( e, n, err );
128 [ + - ][ - + ]: 17784 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
129 [ + - ][ + - ]: 17784 : n = n / n.length(); // Need unit normal
[ + - ]
130 [ + + ]: 88808 : for( i = 0; i < 4; ++i )
131 : : {
132 [ + - ]: 71052 : mCoords[0] = vertices[v_i[locs_hex[i][0]]];
133 [ + - ]: 71052 : mCoords[1] = vertices[v_i[locs_hex[i][1]]];
134 [ + - ]: 71052 : mCoords[2] = vertices[v_i[locs_hex[i][2]]];
135 [ + - ]: 71052 : metric_valid = m_fcn_2i( mMetrics[i], mCoords, n, a2Con, b2Con, c2Con, d_con );
136 [ + + ]: 71052 : if( !metric_valid ) return false;
137 : : }
138 [ + - ]: 17756 : m = average_metrics( mMetrics, 4, err );
139 : 17756 : break;
140 : :
141 : : case TETRAHEDRON:
142 [ + - ]: 57270 : mCoords[0] = vertices[v_i[0]];
143 [ + - ]: 57270 : mCoords[1] = vertices[v_i[1]];
144 [ + - ]: 57270 : mCoords[2] = vertices[v_i[2]];
145 [ + - ]: 57270 : mCoords[3] = vertices[v_i[3]];
146 [ + - ]: 57270 : metric_valid = m_fcn_3e( m, mCoords, a3Con, b3Con, c3Con );
147 [ - + ]: 57270 : if( !metric_valid ) return false;
148 : 57270 : break;
149 : :
150 : : case PYRAMID:
151 [ + + ]: 5690 : for( i = 0; i < 4; ++i )
152 : : {
153 [ + - ]: 4553 : mCoords[0] = vertices[v_i[i]];
154 [ + - ]: 4553 : mCoords[1] = vertices[v_i[( i + 1 ) % 4]];
155 [ + - ]: 4553 : mCoords[2] = vertices[v_i[( i + 3 ) % 4]];
156 [ + - ]: 4553 : mCoords[3] = vertices[v_i[4]];
157 [ + - ]: 4553 : metric_valid = m_fcn_3p( mMetrics[i], mCoords, a3Con, b3Con, c3Con );
158 [ + + ]: 4553 : if( !metric_valid ) return false;
159 : : }
160 [ + - ]: 1137 : m = average_metrics( mMetrics, 4, err );
161 [ + - ][ - + ]: 1137 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
162 : 1137 : break;
163 : :
164 : : case PRISM:
165 [ # # ]: 0 : for( i = 0; i < 6; ++i )
166 : : {
167 [ # # ]: 0 : mCoords[0] = vertices[v_i[locs_pri[i][0]]];
168 [ # # ]: 0 : mCoords[1] = vertices[v_i[locs_pri[i][1]]];
169 [ # # ]: 0 : mCoords[2] = vertices[v_i[locs_pri[i][2]]];
170 [ # # ]: 0 : mCoords[3] = vertices[v_i[locs_pri[i][3]]];
171 [ # # ]: 0 : metric_valid = m_fcn_3w( mMetrics[i], mCoords, a3Con, b3Con, c3Con );
172 [ # # ]: 0 : if( !metric_valid ) return false;
173 : : }
174 [ # # ]: 0 : m = average_metrics( mMetrics, 6, err );
175 [ # # ][ # # ]: 0 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ # # ]
176 : 0 : break;
177 : :
178 : : case HEXAHEDRON:
179 [ + + ]: 71121 : for( i = 0; i < 8; ++i )
180 : : {
181 [ + - ]: 63219 : mCoords[0] = vertices[v_i[locs_hex[i][0]]];
182 [ + - ]: 63219 : mCoords[1] = vertices[v_i[locs_hex[i][1]]];
183 [ + - ]: 63219 : mCoords[2] = vertices[v_i[locs_hex[i][2]]];
184 [ + - ]: 63219 : mCoords[3] = vertices[v_i[locs_hex[i][3]]];
185 [ + - ]: 63219 : metric_valid = m_fcn_3i( mMetrics[i], mCoords, a3Con, b3Con, c3Con, d_con );
186 [ + + ]: 63219 : if( !metric_valid ) return false;
187 : : }
188 [ + - ]: 7902 : m = average_metrics( mMetrics, 8, err );
189 [ + - ][ - + ]: 7902 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
190 : 7902 : break;
191 : :
192 : : default:
193 : : MSQ_SETERR( err )
194 : : ( MsqError::UNSUPPORTED_ELEMENT, "Element type (%d) not supported in IdealWeightInverseMeanRatio",
195 [ # # ][ # # ]: 0 : (int)topo );
196 : 0 : return false;
197 : : } // end switch over element type
198 : 84464 : return true;
199 : : }
200 : :
201 : 152617 : bool IdealWeightInverseMeanRatio::evaluate_with_gradient( PatchData& pd, size_t handle, double& m,
202 : : std::vector< size_t >& indices, std::vector< Vector3D >& g,
203 : : MsqError& err )
204 : : {
205 [ + - ]: 152617 : const MsqMeshEntity* e = &pd.element_by_index( handle );
206 [ + - ]: 152617 : EntityTopology topo = e->get_element_type();
207 : :
208 [ + - ][ - + ]: 152617 : if( !analytical_average_gradient() && topo != TRIANGLE && topo != TETRAHEDRON )
[ # # ][ # # ]
[ - + ]
209 : : {
210 : : static bool print = true;
211 [ # # ]: 0 : if( print )
212 : : {
213 : : MSQ_DBGOUT( 1 ) << "Analyical gradient not available for selected averaging scheme. "
214 : : << "Using (possibly much slower) numerical approximation of gradient"
215 : : << " of quality metric. " << std::endl;
216 : 0 : print = false;
217 : : }
218 [ # # ]: 0 : return QualityMetric::evaluate_with_gradient( pd, handle, m, indices, g, err );
219 : : }
220 : :
221 [ + - ]: 152617 : const MsqVertex* vertices = pd.get_vertex_array( err );
222 [ + - ][ - + ]: 152617 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
223 [ + - ]: 152617 : const size_t* v_i = e->get_vertex_index_array();
224 : :
225 [ + - ]: 152617 : Vector3D n; // Surface normal for 2D objects
226 : :
227 : : // Prism and Hex element descriptions
228 : : static const int locs_pri[6][4] = { { 0, 1, 2, 3 }, { 1, 2, 0, 4 }, { 2, 0, 1, 5 },
229 : : { 3, 5, 4, 0 }, { 4, 3, 5, 1 }, { 5, 4, 3, 2 } };
230 : : static const int locs_hex[8][4] = { { 0, 1, 3, 4 }, { 1, 2, 0, 5 }, { 2, 3, 1, 6 }, { 3, 0, 2, 7 },
231 : : { 4, 7, 5, 0 }, { 5, 4, 6, 1 }, { 6, 5, 7, 2 }, { 7, 6, 4, 3 } };
232 : :
233 [ + - ]: 152617 : const Vector3D d_con( 1.0, 1.0, 1.0 );
234 : :
235 : : int i;
236 : :
237 : 152617 : bool metric_valid = false;
238 [ + - ]: 152617 : const uint32_t fm = fixed_vertex_bitmap( pd, e, indices );
239 : :
240 : 152617 : m = 0.0;
241 : :
242 [ + + + + : 152617 : switch( topo )
+ + - ]
243 : : {
244 : : case TRIANGLE:
245 [ + - ]: 101 : pd.get_domain_normal_at_element( e, n, err );
246 [ + - ][ - + ]: 101 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
247 [ + - ][ + - ]: 101 : n /= n.length(); // Need unit normal
248 [ + - ]: 101 : mCoords[0] = vertices[v_i[0]];
249 [ + - ]: 101 : mCoords[1] = vertices[v_i[1]];
250 [ + - ]: 101 : mCoords[2] = vertices[v_i[2]];
251 [ + - ]: 101 : g.resize( 3 );
252 [ + - ][ + - ]: 101 : if( !g_fcn_2e( m, arrptr( g ), mCoords, n, a2Con, b2Con, c2Con ) ) return false;
[ - + ]
253 : 101 : break;
254 : :
255 : : case QUADRILATERAL:
256 [ + - ]: 532 : pd.get_domain_normal_at_element( e, n, err );
257 [ + - ][ - + ]: 532 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
258 [ + - ][ + - ]: 532 : n /= n.length(); // Need unit normal
259 [ + + ]: 2660 : for( i = 0; i < 4; ++i )
260 : : {
261 [ + - ]: 2128 : mCoords[0] = vertices[v_i[locs_hex[i][0]]];
262 [ + - ]: 2128 : mCoords[1] = vertices[v_i[locs_hex[i][1]]];
263 [ + - ]: 2128 : mCoords[2] = vertices[v_i[locs_hex[i][2]]];
264 [ + - ][ - + ]: 2128 : if( !g_fcn_2i( mMetrics[i], mGradients + 3 * i, mCoords, n, a2Con, b2Con, c2Con, d_con ) ) return false;
265 : : }
266 : :
267 [ + - ]: 532 : g.resize( 4 );
268 [ + - ][ + - ]: 532 : m = average_corner_gradients( QUADRILATERAL, fm, 4, mMetrics, mGradients, arrptr( g ), err );
269 [ + - ][ - + ]: 532 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
270 : 532 : break;
271 : :
272 : : case TETRAHEDRON:
273 [ + - ]: 136059 : mCoords[0] = vertices[v_i[0]];
274 [ + - ]: 136059 : mCoords[1] = vertices[v_i[1]];
275 [ + - ]: 136059 : mCoords[2] = vertices[v_i[2]];
276 [ + - ]: 136059 : mCoords[3] = vertices[v_i[3]];
277 [ + - ]: 136059 : g.resize( 4 );
278 [ + - ][ + - ]: 136059 : metric_valid = g_fcn_3e( m, arrptr( g ), mCoords, a3Con, b3Con, c3Con );
279 [ - + ]: 136059 : if( !metric_valid ) return false;
280 : 136059 : break;
281 : :
282 : : case PYRAMID:
283 [ + + ]: 11382 : for( i = 0; i < 4; ++i )
284 : : {
285 [ + - ]: 9106 : mCoords[0] = vertices[v_i[i]];
286 [ + - ]: 9106 : mCoords[1] = vertices[v_i[( i + 1 ) % 4]];
287 [ + - ]: 9106 : mCoords[2] = vertices[v_i[( i + 3 ) % 4]];
288 [ + - ]: 9106 : mCoords[3] = vertices[v_i[4]];
289 [ + - ]: 9106 : metric_valid = g_fcn_3p( mMetrics[i], mGradients + 4 * i, mCoords, a3Con, b3Con, c3Con );
290 [ + + ]: 9106 : if( !metric_valid ) return false;
291 : : }
292 : :
293 [ + - ]: 2276 : g.resize( 5 );
294 [ + - ][ + - ]: 2276 : m = average_corner_gradients( PYRAMID, fm, 4, mMetrics, mGradients, arrptr( g ), err );
295 [ + - ][ - + ]: 2276 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
296 : 2276 : break;
297 : :
298 : : case PRISM:
299 [ + + ]: 462 : for( i = 0; i < 6; ++i )
300 : : {
301 [ + - ]: 396 : mCoords[0] = vertices[v_i[locs_pri[i][0]]];
302 [ + - ]: 396 : mCoords[1] = vertices[v_i[locs_pri[i][1]]];
303 [ + - ]: 396 : mCoords[2] = vertices[v_i[locs_pri[i][2]]];
304 [ + - ]: 396 : mCoords[3] = vertices[v_i[locs_pri[i][3]]];
305 [ + - ][ - + ]: 396 : if( !g_fcn_3w( mMetrics[i], mGradients + 4 * i, mCoords, a3Con, b3Con, c3Con ) ) return false;
306 : : }
307 : :
308 [ + - ]: 66 : g.resize( 6 );
309 [ + - ][ + - ]: 66 : m = average_corner_gradients( PRISM, fm, 6, mMetrics, mGradients, arrptr( g ), err );
310 [ + - ][ - + ]: 66 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
311 : 66 : break;
312 : :
313 : : case HEXAHEDRON:
314 [ + + ]: 122229 : for( i = 0; i < 8; ++i )
315 : : {
316 [ + - ]: 108648 : mCoords[0] = vertices[v_i[locs_hex[i][0]]];
317 [ + - ]: 108648 : mCoords[1] = vertices[v_i[locs_hex[i][1]]];
318 [ + - ]: 108648 : mCoords[2] = vertices[v_i[locs_hex[i][2]]];
319 [ + - ]: 108648 : mCoords[3] = vertices[v_i[locs_hex[i][3]]];
320 [ + - ][ - + ]: 108648 : if( !g_fcn_3i( mMetrics[i], mGradients + 4 * i, mCoords, a3Con, b3Con, c3Con, d_con ) ) return false;
321 : : }
322 : :
323 [ + - ]: 13581 : g.resize( 8 );
324 [ + - ][ + - ]: 13581 : m = average_corner_gradients( HEXAHEDRON, fm, 8, mMetrics, mGradients, arrptr( g ), err );
325 [ + - ][ - + ]: 13581 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
326 : 13581 : break;
327 : :
328 : : default:
329 : : MSQ_SETERR( err )
330 : : ( MsqError::UNSUPPORTED_ELEMENT, "Element type (%d) not supported in IdealWeightInverseMeanRatio",
331 [ # # ][ # # ]: 0 : (int)topo );
332 : 0 : return false;
333 : : } // end switch over element type
334 : :
335 [ + - ]: 152615 : remove_fixed_gradients( topo, fm, g );
336 : 152617 : return true;
337 : : }
338 : :
339 : 0 : bool IdealWeightInverseMeanRatio::evaluate_with_Hessian_diagonal( PatchData& pd, size_t handle, double& m,
340 : : std::vector< size_t >& indices,
341 : : std::vector< Vector3D >& g,
342 : : std::vector< SymMatrix3D >& h, MsqError& err )
343 : : {
344 [ # # ]: 0 : const MsqMeshEntity* e = &pd.element_by_index( handle );
345 [ # # ]: 0 : EntityTopology topo = e->get_element_type();
346 : :
347 [ # # ][ # # ]: 0 : if( !analytical_average_hessian() && topo != TRIANGLE && topo != TETRAHEDRON )
[ # # ][ # # ]
[ # # ]
348 : : {
349 : : static bool print = true;
350 [ # # ]: 0 : if( print )
351 : : {
352 : : MSQ_DBGOUT( 1 ) << "Analyical gradient not available for selected averaging scheme. "
353 : : << "Using (possibly much slower) numerical approximation of gradient"
354 : : << " of quality metric. " << std::endl;
355 : 0 : print = false;
356 : : }
357 [ # # ]: 0 : return QualityMetric::evaluate_with_Hessian_diagonal( pd, handle, m, indices, g, h, err );
358 : : }
359 : :
360 [ # # ]: 0 : const MsqVertex* vertices = pd.get_vertex_array( err );
361 [ # # ][ # # ]: 0 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ # # ]
362 [ # # ]: 0 : const size_t* v_i = e->get_vertex_index_array();
363 : :
364 [ # # ]: 0 : Vector3D n; // Surface normal for 2D objects
365 : :
366 : : // Prism and Hex element descriptions
367 : : static const int locs_pri[6][4] = { { 0, 1, 2, 3 }, { 1, 2, 0, 4 }, { 2, 0, 1, 5 },
368 : : { 3, 5, 4, 0 }, { 4, 3, 5, 1 }, { 5, 4, 3, 2 } };
369 : : static const int locs_hex[8][4] = { { 0, 1, 3, 4 }, { 1, 2, 0, 5 }, { 2, 3, 1, 6 }, { 3, 0, 2, 7 },
370 : : { 4, 7, 5, 0 }, { 5, 4, 6, 1 }, { 6, 5, 7, 2 }, { 7, 6, 4, 3 } };
371 : :
372 [ # # ]: 0 : const Vector3D d_con( 1.0, 1.0, 1.0 );
373 : :
374 : : int i;
375 : :
376 : 0 : bool metric_valid = false;
377 [ # # ]: 0 : const uint32_t fm = fixed_vertex_bitmap( pd, e, indices );
378 : :
379 : 0 : m = 0.0;
380 : :
381 [ # # # # : 0 : switch( topo )
# # # ]
382 : : {
383 : : case TRIANGLE:
384 [ # # ]: 0 : pd.get_domain_normal_at_element( e, n, err );
385 [ # # ][ # # ]: 0 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ # # ]
386 [ # # ][ # # ]: 0 : n /= n.length(); // Need unit normal
387 [ # # ]: 0 : mCoords[0] = vertices[v_i[0]];
388 [ # # ]: 0 : mCoords[1] = vertices[v_i[1]];
389 [ # # ]: 0 : mCoords[2] = vertices[v_i[2]];
390 [ # # ]: 0 : g.resize( 3 );
391 [ # # ]: 0 : h.resize( 3 );
392 [ # # ][ # # ]: 0 : if( !h_fcn_2e( m, arrptr( g ), mHessians, mCoords, n, a2Con, b2Con, c2Con ) ) return false;
[ # # ]
393 [ # # ][ # # ]: 0 : h[0] = mHessians[0].upper();
394 [ # # ][ # # ]: 0 : h[1] = mHessians[3].upper();
395 [ # # ][ # # ]: 0 : h[2] = mHessians[5].upper();
396 : 0 : break;
397 : :
398 : : case QUADRILATERAL:
399 [ # # ]: 0 : pd.get_domain_normal_at_element( e, n, err );
400 [ # # ][ # # ]: 0 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ # # ]
401 [ # # ][ # # ]: 0 : n /= n.length(); // Need unit normal
402 [ # # ]: 0 : for( i = 0; i < 4; ++i )
403 : : {
404 [ # # ]: 0 : mCoords[0] = vertices[v_i[locs_hex[i][0]]];
405 [ # # ]: 0 : mCoords[1] = vertices[v_i[locs_hex[i][1]]];
406 [ # # ]: 0 : mCoords[2] = vertices[v_i[locs_hex[i][2]]];
407 [ # # ]: 0 : if( !h_fcn_2i( mMetrics[i], mGradients + 3 * i, mHessians + 6 * i, mCoords, n, a2Con, b2Con, c2Con,
408 [ # # ]: 0 : d_con ) )
409 : 0 : return false;
410 : : }
411 : :
412 [ # # ]: 0 : g.resize( 4 );
413 [ # # ]: 0 : h.resize( 4 );
414 : : m = average_corner_hessian_diagonals( QUADRILATERAL, fm, 4, mMetrics, mGradients, mHessians, arrptr( g ),
415 [ # # ][ # # ]: 0 : arrptr( h ), err );
[ # # ]
416 [ # # ][ # # ]: 0 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ # # ]
417 : 0 : break;
418 : :
419 : : case TETRAHEDRON:
420 [ # # ]: 0 : mCoords[0] = vertices[v_i[0]];
421 [ # # ]: 0 : mCoords[1] = vertices[v_i[1]];
422 [ # # ]: 0 : mCoords[2] = vertices[v_i[2]];
423 [ # # ]: 0 : mCoords[3] = vertices[v_i[3]];
424 [ # # ]: 0 : g.resize( 4 );
425 [ # # ]: 0 : h.resize( 4 );
426 [ # # ][ # # ]: 0 : metric_valid = h_fcn_3e( m, arrptr( g ), mHessians, mCoords, a3Con, b3Con, c3Con );
427 [ # # ]: 0 : if( !metric_valid ) return false;
428 [ # # ][ # # ]: 0 : h[0] = mHessians[0].upper();
429 [ # # ][ # # ]: 0 : h[1] = mHessians[4].upper();
430 [ # # ][ # # ]: 0 : h[2] = mHessians[7].upper();
431 [ # # ][ # # ]: 0 : h[3] = mHessians[9].upper();
432 : 0 : break;
433 : :
434 : : case PYRAMID:
435 [ # # ]: 0 : for( i = 0; i < 4; ++i )
436 : : {
437 [ # # ]: 0 : mCoords[0] = vertices[v_i[i]];
438 [ # # ]: 0 : mCoords[1] = vertices[v_i[( i + 1 ) % 4]];
439 [ # # ]: 0 : mCoords[2] = vertices[v_i[( i + 3 ) % 4]];
440 [ # # ]: 0 : mCoords[3] = vertices[v_i[4]];
441 : : metric_valid =
442 [ # # ]: 0 : h_fcn_3p( mMetrics[i], mGradients + 4 * i, mHessians + 10 * i, mCoords, a3Con, b3Con, c3Con );
443 : :
444 [ # # ]: 0 : if( !metric_valid ) return false;
445 : : }
446 : :
447 [ # # ]: 0 : g.resize( 5 );
448 [ # # ]: 0 : h.resize( 5 );
449 : : m = average_corner_hessian_diagonals( PYRAMID, fm, 4, mMetrics, mGradients, mHessians, arrptr( g ),
450 [ # # ][ # # ]: 0 : arrptr( h ), err );
[ # # ]
451 [ # # ][ # # ]: 0 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ # # ]
452 : 0 : break;
453 : :
454 : : case PRISM:
455 [ # # ]: 0 : for( i = 0; i < 6; ++i )
456 : : {
457 [ # # ]: 0 : mCoords[0] = vertices[v_i[locs_pri[i][0]]];
458 [ # # ]: 0 : mCoords[1] = vertices[v_i[locs_pri[i][1]]];
459 [ # # ]: 0 : mCoords[2] = vertices[v_i[locs_pri[i][2]]];
460 [ # # ]: 0 : mCoords[3] = vertices[v_i[locs_pri[i][3]]];
461 [ # # ][ # # ]: 0 : if( !h_fcn_3w( mMetrics[i], mGradients + 4 * i, mHessians + 10 * i, mCoords, a3Con, b3Con, c3Con ) )
462 : 0 : return false;
463 : : }
464 : :
465 [ # # ]: 0 : g.resize( 6 );
466 [ # # ]: 0 : h.resize( 6 );
467 : : m = average_corner_hessian_diagonals( PRISM, fm, 6, mMetrics, mGradients, mHessians, arrptr( g ),
468 [ # # ][ # # ]: 0 : arrptr( h ), err );
[ # # ]
469 [ # # ][ # # ]: 0 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ # # ]
470 : 0 : break;
471 : :
472 : : case HEXAHEDRON:
473 [ # # ]: 0 : for( i = 0; i < 8; ++i )
474 : : {
475 [ # # ]: 0 : mCoords[0] = vertices[v_i[locs_hex[i][0]]];
476 [ # # ]: 0 : mCoords[1] = vertices[v_i[locs_hex[i][1]]];
477 [ # # ]: 0 : mCoords[2] = vertices[v_i[locs_hex[i][2]]];
478 [ # # ]: 0 : mCoords[3] = vertices[v_i[locs_hex[i][3]]];
479 [ # # ]: 0 : if( !h_fcn_3i( mMetrics[i], mGradients + 4 * i, mHessians + 10 * i, mCoords, a3Con, b3Con, c3Con,
480 [ # # ]: 0 : d_con ) )
481 : 0 : return false;
482 : : }
483 : :
484 [ # # ]: 0 : g.resize( 8 );
485 [ # # ]: 0 : h.resize( 8 );
486 : : m = average_corner_hessian_diagonals( HEXAHEDRON, fm, 8, mMetrics, mGradients, mHessians, arrptr( g ),
487 [ # # ][ # # ]: 0 : arrptr( h ), err );
[ # # ]
488 [ # # ][ # # ]: 0 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ # # ]
489 : 0 : break;
490 : :
491 : : default:
492 : : MSQ_SETERR( err )
493 : : ( MsqError::UNSUPPORTED_ELEMENT, "Element type (%d) not supported in IdealWeightInverseMeanRatio",
494 [ # # ][ # # ]: 0 : (int)topo );
495 : 0 : return false;
496 : : } // end switch over element type
497 : :
498 [ # # ]: 0 : remove_fixed_diagonals( topo, fm, g, h );
499 : 0 : return true;
500 : : }
501 : :
502 : 139646 : bool IdealWeightInverseMeanRatio::evaluate_with_Hessian( PatchData& pd, size_t handle, double& m,
503 : : std::vector< size_t >& indices, std::vector< Vector3D >& g,
504 : : std::vector< Matrix3D >& h, MsqError& err )
505 : : {
506 [ + - ]: 139646 : const MsqMeshEntity* e = &pd.element_by_index( handle );
507 [ + - ]: 139646 : EntityTopology topo = e->get_element_type();
508 : :
509 [ + - ][ - + ]: 139646 : if( !analytical_average_hessian() && topo != TRIANGLE && topo != TETRAHEDRON )
[ # # ][ # # ]
[ - + ]
510 : : {
511 : : static bool print = true;
512 [ # # ]: 0 : if( print )
513 : : {
514 : : MSQ_DBGOUT( 1 ) << "Analyical gradient not available for selected averaging scheme. "
515 : : << "Using (possibly much slower) numerical approximation of gradient"
516 : : << " of quality metric. " << std::endl;
517 : 0 : print = false;
518 : : }
519 [ # # ]: 0 : return QualityMetric::evaluate_with_Hessian( pd, handle, m, indices, g, h, err );
520 : : }
521 : :
522 [ + - ]: 139646 : const MsqVertex* vertices = pd.get_vertex_array( err );
523 [ + - ][ - + ]: 139646 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
524 [ + - ]: 139646 : const size_t* v_i = e->get_vertex_index_array();
525 : :
526 [ + - ]: 139646 : Vector3D n; // Surface normal for 2D objects
527 : :
528 : : // Prism and Hex element descriptions
529 : : static const int locs_pri[6][4] = { { 0, 1, 2, 3 }, { 1, 2, 0, 4 }, { 2, 0, 1, 5 },
530 : : { 3, 5, 4, 0 }, { 4, 3, 5, 1 }, { 5, 4, 3, 2 } };
531 : : static const int locs_hex[8][4] = { { 0, 1, 3, 4 }, { 1, 2, 0, 5 }, { 2, 3, 1, 6 }, { 3, 0, 2, 7 },
532 : : { 4, 7, 5, 0 }, { 5, 4, 6, 1 }, { 6, 5, 7, 2 }, { 7, 6, 4, 3 } };
533 : :
534 [ + - ]: 139646 : const Vector3D d_con( 1.0, 1.0, 1.0 );
535 : :
536 : : int i;
537 : :
538 : 139646 : bool metric_valid = false;
539 [ + - ]: 139646 : const uint32_t fm = fixed_vertex_bitmap( pd, e, indices );
540 : :
541 : 139646 : m = 0.0;
542 : :
543 [ - - + + : 139646 : switch( topo )
+ + - ]
544 : : {
545 : : case TRIANGLE:
546 [ # # ]: 0 : pd.get_domain_normal_at_element( e, n, err );
547 [ # # ][ # # ]: 0 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ # # ]
548 [ # # ][ # # ]: 0 : n /= n.length(); // Need unit normal
549 [ # # ]: 0 : mCoords[0] = vertices[v_i[0]];
550 [ # # ]: 0 : mCoords[1] = vertices[v_i[1]];
551 [ # # ]: 0 : mCoords[2] = vertices[v_i[2]];
552 [ # # ]: 0 : g.resize( 3 );
553 [ # # ]: 0 : h.resize( 6 );
554 [ # # ][ # # ]: 0 : if( !h_fcn_2e( m, arrptr( g ), arrptr( h ), mCoords, n, a2Con, b2Con, c2Con ) ) return false;
[ # # ][ # # ]
555 : 0 : break;
556 : :
557 : : case QUADRILATERAL:
558 [ # # ]: 0 : pd.get_domain_normal_at_element( e, n, err );
559 [ # # ][ # # ]: 0 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ # # ]
560 [ # # ][ # # ]: 0 : n /= n.length(); // Need unit normal
561 [ # # ]: 0 : for( i = 0; i < 4; ++i )
562 : : {
563 [ # # ]: 0 : mCoords[0] = vertices[v_i[locs_hex[i][0]]];
564 [ # # ]: 0 : mCoords[1] = vertices[v_i[locs_hex[i][1]]];
565 [ # # ]: 0 : mCoords[2] = vertices[v_i[locs_hex[i][2]]];
566 [ # # ]: 0 : if( !h_fcn_2i( mMetrics[i], mGradients + 3 * i, mHessians + 6 * i, mCoords, n, a2Con, b2Con, c2Con,
567 [ # # ]: 0 : d_con ) )
568 : 0 : return false;
569 : : }
570 : :
571 [ # # ]: 0 : g.resize( 4 );
572 [ # # ]: 0 : h.resize( 10 );
573 : : m = average_corner_hessians( QUADRILATERAL, fm, 4, mMetrics, mGradients, mHessians, arrptr( g ),
574 [ # # ][ # # ]: 0 : arrptr( h ), err );
[ # # ]
575 [ # # ][ # # ]: 0 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ # # ]
576 : 0 : break;
577 : :
578 : : case TETRAHEDRON:
579 [ + - ]: 125015 : mCoords[0] = vertices[v_i[0]];
580 [ + - ]: 125015 : mCoords[1] = vertices[v_i[1]];
581 [ + - ]: 125015 : mCoords[2] = vertices[v_i[2]];
582 [ + - ]: 125015 : mCoords[3] = vertices[v_i[3]];
583 [ + - ]: 125015 : g.resize( 4 );
584 [ + - ]: 125015 : h.resize( 10 );
585 [ + - ][ + - ]: 125015 : metric_valid = h_fcn_3e( m, arrptr( g ), arrptr( h ), mCoords, a3Con, b3Con, c3Con );
[ + - ]
586 [ - + ]: 125015 : if( !metric_valid ) return false;
587 : 125015 : break;
588 : :
589 : : case PYRAMID:
590 [ + + ]: 11050 : for( i = 0; i < 4; ++i )
591 : : {
592 [ + - ]: 8840 : mCoords[0] = vertices[v_i[i]];
593 [ + - ]: 8840 : mCoords[1] = vertices[v_i[( i + 1 ) % 4]];
594 [ + - ]: 8840 : mCoords[2] = vertices[v_i[( i + 3 ) % 4]];
595 [ + - ]: 8840 : mCoords[3] = vertices[v_i[4]];
596 : : metric_valid =
597 [ + - ]: 8840 : h_fcn_3p( mMetrics[i], mGradients + 4 * i, mHessians + 10 * i, mCoords, a3Con, b3Con, c3Con );
598 : :
599 [ - + ]: 8840 : if( !metric_valid ) return false;
600 : : }
601 : :
602 [ + - ]: 2210 : g.resize( 5 );
603 [ + - ]: 2210 : h.resize( 15 );
604 : : m = average_corner_hessians( PYRAMID, fm, 4, mMetrics, mGradients, mHessians, arrptr( g ), arrptr( h ),
605 [ + - ][ + - ]: 2210 : err );
[ + - ]
606 [ + - ][ - + ]: 2210 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
607 : 2210 : break;
608 : :
609 : : case PRISM:
610 [ + + ]: 462 : for( i = 0; i < 6; ++i )
611 : : {
612 [ + - ]: 396 : mCoords[0] = vertices[v_i[locs_pri[i][0]]];
613 [ + - ]: 396 : mCoords[1] = vertices[v_i[locs_pri[i][1]]];
614 [ + - ]: 396 : mCoords[2] = vertices[v_i[locs_pri[i][2]]];
615 [ + - ]: 396 : mCoords[3] = vertices[v_i[locs_pri[i][3]]];
616 [ + - ][ - + ]: 396 : if( !h_fcn_3w( mMetrics[i], mGradients + 4 * i, mHessians + 10 * i, mCoords, a3Con, b3Con, c3Con ) )
617 : 0 : return false;
618 : : }
619 : :
620 [ + - ]: 66 : g.resize( 6 );
621 [ + - ]: 66 : h.resize( 21 );
622 [ + - ][ + - ]: 66 : m = average_corner_hessians( PRISM, fm, 6, mMetrics, mGradients, mHessians, arrptr( g ), arrptr( h ), err );
[ + - ]
623 [ + - ][ - + ]: 66 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
624 : 66 : break;
625 : :
626 : : case HEXAHEDRON:
627 [ + + ]: 111195 : for( i = 0; i < 8; ++i )
628 : : {
629 [ + - ]: 98840 : mCoords[0] = vertices[v_i[locs_hex[i][0]]];
630 [ + - ]: 98840 : mCoords[1] = vertices[v_i[locs_hex[i][1]]];
631 [ + - ]: 98840 : mCoords[2] = vertices[v_i[locs_hex[i][2]]];
632 [ + - ]: 98840 : mCoords[3] = vertices[v_i[locs_hex[i][3]]];
633 [ - + ]: 98840 : if( !h_fcn_3i( mMetrics[i], mGradients + 4 * i, mHessians + 10 * i, mCoords, a3Con, b3Con, c3Con,
634 [ + - ]: 98840 : d_con ) )
635 : 0 : return false;
636 : : }
637 : :
638 [ + - ]: 12355 : g.resize( 8 );
639 [ + - ]: 12355 : h.resize( 36 );
640 : : m = average_corner_hessians( HEXAHEDRON, fm, 8, mMetrics, mGradients, mHessians, arrptr( g ), arrptr( h ),
641 [ + - ][ + - ]: 12355 : err );
[ + - ]
642 [ + - ][ - + ]: 12355 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
643 : 12355 : break;
644 : :
645 : : default:
646 : : MSQ_SETERR( err )
647 : : ( MsqError::UNSUPPORTED_ELEMENT, "Element type (%d) not supported in IdealWeightInverseMeanRatio",
648 [ # # ][ # # ]: 0 : (int)topo );
649 : 0 : return false;
650 : : } // end switch over element type
651 : :
652 [ + - ]: 139646 : remove_fixed_gradients( topo, fm, g );
653 [ + - ]: 139646 : remove_fixed_hessians( topo, fm, h );
654 : 139646 : return true;
655 : : }
656 : :
657 [ + - ][ + - ]: 56 : } // namespace MBMesquite
|
