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