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) kraftche@cae.wisc.edu
24
25 ***************************************************************** */
26
27 #include "QualityMetricTester.hpp"
28 #include "PatchData.hpp"
29 #include "QualityMetric.hpp"
30 #include "IdealElements.hpp"
31 #include "UnitUtil.hpp"
32 #include "ElemSampleQM.hpp"
33 #include "TopologyInfo.hpp"
34 #include "EdgeQM.hpp"
35
36 #include <cppunit/extensions/HelperMacros.h>
37
types_in_group(QualityMetricTester::ElemTypeGroup group)38 static std::vector<EntityTopology> types_in_group( QualityMetricTester::ElemTypeGroup group )
39 {
40 std::vector<EntityTopology> types;
41 switch (group) {
42 default:
43 case QualityMetricTester::ALL:
44 types.push_back( POLYHEDRON );
45 types.push_back( POLYGON );
46 case QualityMetricTester::ALL_FE:
47 types.push_back( SEPTAHEDRON );
48 case QualityMetricTester::ALL_FE_EXCEPT_SEPTAHEDRON:
49 types.push_back( PYRAMID );
50 types.push_back( PRISM );
51 case QualityMetricTester::NON_MIXED_FE:
52 types.push_back( HEXAHEDRON );
53 types.push_back( QUADRILATERAL );
54 case QualityMetricTester::SIMPLICIES:
55 types.push_back(TETRAHEDRON);
56 types.push_back(TRIANGLE);
57 break;
58
59 case QualityMetricTester::THREE_D:
60 types.push_back( POLYHEDRON );
61 case QualityMetricTester::THREE_D_FE:
62 types.push_back( SEPTAHEDRON );
63 case QualityMetricTester::THREE_D_FE_EXCEPT_SEPTAHEDRON:
64 types.push_back( PYRAMID );
65 types.push_back( PRISM );
66 case QualityMetricTester::THREE_D_NON_MIXED_FE:
67 types.push_back( HEXAHEDRON );
68 types.push_back(TETRAHEDRON);
69 break;
70
71 case QualityMetricTester::TWO_D:
72 types.push_back( POLYGON );
73 case QualityMetricTester::TWO_D_FE:
74 types.push_back( TRIANGLE );
75 types.push_back( QUADRILATERAL );
76 break;
77 }
78 std::reverse( types.begin(), types.end() );
79 return types;
80 }
81
QualityMetricTester(const EntityTopology * supported_elem_types,size_t len,const Settings * settings)82 QualityMetricTester::QualityMetricTester(
83 const EntityTopology* supported_elem_types,
84 size_t len,
85 const Settings* settings )
86 : degenHexPyramid(false),
87 types( len ),
88 mSettings( settings ),
89 geomPlane( Vector3D(0,0,1), Vector3D(0,0,0) )
90 {
91 std::copy( supported_elem_types, supported_elem_types+len, types.begin() );
92 }
93
QualityMetricTester(ElemTypeGroup group,const Settings * settings)94 QualityMetricTester::QualityMetricTester(
95 ElemTypeGroup group,
96 const Settings* settings )
97 : degenHexPyramid(false),
98 types( types_in_group(group) ),
99 mSettings( settings ),
100 geomPlane( Vector3D(0,0,1), Vector3D(0,0,0) )
101 {
102 }
103
get_ideal_tris(PatchData & pd,bool unit_area)104 void QualityMetricTester::get_ideal_tris( PatchData& pd, bool unit_area )
105 {
106 // 6 ------- 5 .
107 // /\ /\ .
108 // / \ / \ .
109 // / \ / \ .
110 // 1 <------X------> 4 .
111 // \ /0\ / .
112 // \ / \ / .
113 // \/ \/ .
114 // 2 ------- 3 .
115 static const double y3 = MSQ_SQRT_THREE_DIV_TWO;
116 double ideal_tri_verts[] = { 0.0, 0.0, 0.0,
117 -1.0, 0.0, 0.0,
118 -0.5, -y3, 0.0,
119 0.5, -y3, 0.0,
120 1.0, 0.0, 0.0,
121 0.5, y3, 0.0,
122 -0.5, y3, 0.0 };
123 const size_t ideal_tri_elems[] = { 0, 1, 2,
124 0, 2, 3,
125 0, 4, 5,
126 0, 4, 5,
127 0, 5, 6,
128 0, 6, 1 };
129 const bool fixed[] = { false, true, true, true, true, true, true };
130
131 if (unit_area) {
132 const double unit_tri_scale = 2.0/std::sqrt(6.0);
133 for (size_t i = 0; i < sizeof(ideal_tri_verts)/sizeof(double); ++i)
134 ideal_tri_verts[i] *= unit_tri_scale;
135 }
136
137 pd.set_domain( &geomPlane );
138
139 MsqPrintError err( std::cout );
140 pd.fill( 7, ideal_tri_verts, 6, TRIANGLE, ideal_tri_elems, fixed, err );
141 CPPUNIT_ASSERT(!MSQ_CHKERR(err));
142 }
143
get_ideal_quads(PatchData & pd)144 void QualityMetricTester::get_ideal_quads( PatchData& pd )
145 {
146 // 7
147 // 8 +------+------+ 6
148 // | | |
149 // | | |
150 // | |0 |
151 // 1 +------+------+ 5
152 // | | |
153 // | | |
154 // | | |
155 // 2 +------+------+ 4
156 // 3
157 static const double ideal_quad_verts[] = { 0.0, 0.0, 0.0,
158 -1.0, 0.0, 0.0,
159 -1.0,-1.0, 0.0,
160 0.0,-1.0, 0.0,
161 1.0,-1.0, 0.0,
162 1.0, 0.0, 0.0,
163 1.0, 1.0, 0.0,
164 0.0, 1.0, 0.0,
165 -1.0, 1.0, 0.0 };
166 static const size_t ideal_quad_elems[] = { 0, 1, 2, 3,
167 0, 3, 4, 5,
168 0, 5, 6, 7,
169 0, 7, 8, 1 };
170 const bool fixed[] = { false, true, true, true, true, true, true, true, true };
171
172 pd.set_domain( &geomPlane );
173
174 MsqPrintError err( std::cout );
175 pd.fill( 9, ideal_quad_verts, 4, QUADRILATERAL, ideal_quad_elems, fixed, err );
176 CPPUNIT_ASSERT(!MSQ_CHKERR(err));
177 }
178
179
get_ideal_hexes(PatchData & pd)180 void QualityMetricTester::get_ideal_hexes( PatchData& pd )
181 {
182 static const double ideal_hex_verts[] = { 0.0, 0.0, 0.0,
183 -1.0, 0.0, 0.0,
184 -1.0,-1.0, 0.0,
185 0.0,-1.0, 0.0,
186 1.0,-1.0, 0.0,
187 1.0, 0.0, 0.0,
188 1.0, 1.0, 0.0,
189 0.0, 1.0, 0.0,
190 -1.0, 1.0, 0.0,
191 0.0, 0.0,-1.0,
192 -1.0, 0.0,-1.0,
193 -1.0,-1.0,-1.0,
194 0.0,-1.0,-1.0,
195 1.0,-1.0,-1.0,
196 1.0, 0.0,-1.0,
197 1.0, 1.0,-1.0,
198 0.0, 1.0,-1.0,
199 -1.0, 1.0,-1.0,
200 0.0, 0.0, 1.0,
201 -1.0, 0.0, 1.0,
202 -1.0,-1.0, 1.0,
203 0.0,-1.0, 1.0,
204 1.0,-1.0, 1.0,
205 1.0, 0.0, 1.0,
206 1.0, 1.0, 1.0,
207 0.0, 1.0, 1.0,
208 -1.0, 1.0, 1.0 };
209 static const size_t ideal_hex_elems[] = { 9, 10, 11, 12, 0, 1, 2, 3,
210 9, 12, 13, 14, 0, 3, 4, 5,
211 9, 14, 15, 16, 0, 5, 6, 7,
212 9, 16, 17, 10, 0, 7, 8, 1,
213 0, 1, 2, 3, 18, 19, 20, 21,
214 0, 3, 4, 5, 18, 21, 22, 23,
215 0, 5, 6, 7, 18, 23, 24, 25,
216 0, 7, 8, 1, 18, 25, 26, 19 };
217 bool fixed[27];
218 fixed[0] = false;
219 for (size_t f = 1; f < 27; ++f)
220 fixed[f] = true;
221
222 MsqPrintError err( std::cout );
223 pd.fill( 27, ideal_hex_verts, 8, HEXAHEDRON, ideal_hex_elems, fixed, err );
224 CPPUNIT_ASSERT(!MSQ_CHKERR(err));
225 }
226
227 // free_vertex_index:
228 // >= 0 : all vertice except specified one are fixed
229 // -1 : all vertices are free
230 // -2 : first vertex is fixed, all others are free
get_ideal_element(EntityTopology type,bool unit_area,PatchData & pd,int free_vertex_index)231 void QualityMetricTester::get_ideal_element( EntityTopology type,
232 bool unit_area,
233 PatchData& pd,
234 int free_vertex_index )
235 {
236 if (TopologyInfo::dimension(type) == 2)
237 pd.set_domain( &geomPlane );
238
239 const size_t n = TopologyInfo::corners(type);
240 const Vector3D* coords = unit_area ? unit_element(type,degenHexPyramid) : unit_edge_element(type,degenHexPyramid);
241 CPPUNIT_ASSERT(coords != 0);
242
243 CPPUNIT_ASSERT(sizeof(double)*3 == sizeof(Vector3D));
244 const double* elem_verts = reinterpret_cast<const double*>(coords);
245
246 bool* fixed = new bool[n];
247 std::vector<size_t> conn(n);
248 for (size_t i = 0; i < n; ++i) {
249 conn[i] = i;
250 fixed[i] = (free_vertex_index >= 0);
251 }
252 if (free_vertex_index == -2)
253 fixed[0] = true;
254 else if (free_vertex_index >= 0)
255 fixed[free_vertex_index] = false;
256
257 MsqPrintError err( std::cout );
258 pd.fill( n, elem_verts, 1, type, arrptr(conn), fixed, err );
259 delete [] fixed;
260 CPPUNIT_ASSERT(!MSQ_CHKERR(err));
261 }
get_ideal_element(EntityTopology type,bool unit_area,PatchData & pd,bool first_vertex_fixed)262 void QualityMetricTester::get_ideal_element( EntityTopology type,
263 bool unit_area,
264 PatchData& pd,
265 bool first_vertex_fixed )
266 {
267 get_ideal_element( type, unit_area, pd, first_vertex_fixed ? -2 : -1 );
268 }
269
270 /** Begin with an ideal element, and move the first vertex one half
271 * of the distnace to the centroid */
get_nonideal_element(EntityTopology type,PatchData & pd,int free_vtx_type)272 void QualityMetricTester::get_nonideal_element( EntityTopology type,
273 PatchData& pd,
274 int free_vtx_type )
275 {
276 MsqError err;
277 get_ideal_element( type, false, pd, free_vtx_type );
278 const size_t* verts = pd.element_by_index(0).get_vertex_index_array();
279 const MsqVertex* coords = pd.get_vertex_array(err);
280 size_t n = pd.element_by_index(0).vertex_count();
281
282 Vector3D sum(0,0,0);
283 for (unsigned i = 0; i < n; ++i)
284 sum += coords[verts[i]];
285
286 sum /= n;
287 sum += coords[verts[0]];
288 sum *= 0.5;
289 pd.set_vertex_coordinates( sum, verts[0], err );
290 }
291
get_nonideal_element(EntityTopology type,PatchData & pd,bool first_vertex_fixed)292 void QualityMetricTester::get_nonideal_element( EntityTopology type,
293 PatchData& pd,
294 bool first_vertex_fixed )
295 {
296 get_nonideal_element( type, pd, first_vertex_fixed ? -2 : -1 );
297 }
298
299 /** Collapse first and second vertices of an ideal element */
get_degenerate_element(EntityTopology type,PatchData & pd)300 void QualityMetricTester::get_degenerate_element( EntityTopology type, PatchData& pd )
301 {
302 MsqError err;
303 get_ideal_element( type, false, pd );
304 const size_t* verts = pd.element_by_index(0).get_vertex_index_array();
305 const MsqVertex* coords = pd.get_vertex_array(err);
306 pd.set_vertex_coordinates( coords[verts[1]], verts[0], err );
307 }
308
309 /** Create element with zero area/volume
310 * For quads and hexes, the results are also inverted elements.
311 */
get_zero_element(EntityTopology type,PatchData & pd)312 void QualityMetricTester::get_zero_element( EntityTopology type, PatchData& pd )
313 {
314 MsqError err;
315 get_ideal_element( type, false, pd );
316 MsqMeshEntity& elem = pd.element_by_index(0);
317 const size_t* verts = elem.get_vertex_index_array();
318 const MsqVertex* coords = pd.get_vertex_array(err);
319 unsigned i;
320 Vector3D sum(0,0,0);
321
322 switch (type) {
323 case TRIANGLE:
324 pd.set_vertex_coordinates( 0.5*(coords[verts[0]]+coords[verts[1]]), verts[2], err );
325 break;
326 case QUADRILATERAL:
327 pd.set_vertex_coordinates( coords[verts[0]], verts[1], err );
328 pd.set_vertex_coordinates( coords[verts[3]], verts[2], err );
329 break;
330 case TETRAHEDRON:
331 for (i = 0; i < 3; ++i)
332 sum += coords[verts[i]];
333 pd.set_vertex_coordinates( sum/3.0, verts[3], err );
334 break;
335 case PYRAMID:
336 for (i = 0; i < 4; ++i)
337 sum += coords[verts[i]];
338 pd.set_vertex_coordinates( 0.25*sum, verts[4], err );
339 break;
340 case PRISM:
341 pd.set_vertex_coordinates( 0.5*(coords[verts[0]]+coords[verts[1]]), verts[2], err );
342 pd.set_vertex_coordinates( 0.5*(coords[verts[3]]+coords[verts[4]]), verts[5], err );
343 break;
344 case HEXAHEDRON:
345 pd.set_vertex_coordinates( coords[verts[0]], verts[4], err );
346 pd.set_vertex_coordinates( coords[verts[1]], verts[5], err );
347 pd.set_vertex_coordinates( coords[verts[2]], verts[6], err );
348 pd.set_vertex_coordinates( coords[verts[3]], verts[7], err );
349 break;
350 default:
351 CPPUNIT_ASSERT(false);
352 break;
353 }
354 }
355
356 /** Create inverted elements.
357 * For tri and tet elements, reflect one vertex about the opposite side.
358 * For all other elements, introduce a concave vertex in one of the
359 * quadrilateral faces.
360 */
get_inverted_element(EntityTopology type,PatchData & pd)361 void QualityMetricTester::get_inverted_element( EntityTopology type, PatchData& pd )
362 {
363 MsqError err;
364 get_ideal_element( type, false, pd );
365 MsqMeshEntity& elem = pd.element_by_index(0);
366 const size_t* verts = elem.get_vertex_index_array();
367 const MsqVertex* coords = pd.get_vertex_array(err);
368 unsigned i;
369 Vector3D sum(0,0,0);
370
371 // switch (type) {
372 // case TRIANGLE:
373 // coords[verts[2]] = coords[verts[0]] + coords[verts[1]] - coords[verts[2]];
374 // break;
375 // case QUADRILATERAL:
376 // case PYRAMID:
377 // case HEXAHEDRON:
378 // coords[verts[2]] += 0.75 * (coords[verts[0]] - coords[verts[2]]);
379 // break;
380 // case TETRAHEDRON:
381 // for (i = 0; i < 3; ++i)
382 // sum += coords[verts[i]];
383 // coords[verts[3]] += 0.5*sum - 1.5*coords[verts[3]];
384 // break;
385 // case PRISM:
386 // coords[verts[4]] += 0.75 * (coords[verts[0]] - coords[verts[4]]);
387 // break;
388 // default:
389 // CPPUNIT_ASSERT(false);
390 // break;
391 // }
392 if (type == TRIANGLE)
393 pd.set_vertex_coordinates( coords[verts[0]] + coords[verts[1]] - coords[verts[2]], verts[2], err );
394 else if (type == QUADRILATERAL || type == PYRAMID || type == HEXAHEDRON)
395 pd.set_vertex_coordinates( 0.75 * (coords[verts[0]] - coords[verts[2]]), verts[2], err );
396 else if (type == TETRAHEDRON) {
397 for (i = 0; i < 3; ++i)
398 sum += coords[verts[i]];
399 pd.set_vertex_coordinates( 0.5*sum - 1.5*coords[verts[3]], verts[3], err );
400 }
401 else if (type == PRISM)
402 pd.set_vertex_coordinates( 0.75 * (coords[verts[0]] - coords[verts[4]]), verts[4], err );
403 else
404 CPPUNIT_ASSERT(false);
405
406 if (TopologyInfo::dimension(type) == 3 || pd.domain_set()) {
407 int inverted, total;
408 elem.check_element_orientation(pd,inverted,total,err);
409 CPPUNIT_ASSERT(!err);
410 CPPUNIT_ASSERT(total);
411 }
412 }
413
test_type_is_supported(EntityTopology type,QualityMetric * qm)414 void QualityMetricTester::test_type_is_supported( EntityTopology type, QualityMetric* qm )
415 {
416 MsqPrintError err(std::cout);
417 bool rval;
418 // create patch containing only elements of sepcified type
419 PatchData pd;
420 if (mSettings)
421 pd.attach_settings(mSettings);
422 get_ideal_element( type, false, pd );
423 // get list of evaluation locations in patch
424 std::vector<size_t> handles;
425 qm->get_evaluations( pd, handles, false, err );
426 CPPUNIT_ASSERT(!err);
427 CPPUNIT_ASSERT(!handles.empty());
428 // evaluate the metric at one location
429 double value = -HUGE_VAL;
430 rval = qm->evaluate( pd, handles[0], value, err );
431 CPPUNIT_ASSERT(!err);
432 CPPUNIT_ASSERT(rval);
433 // if metric returned input value for 'value' try again
434 // with a different initial value to verify that the metric
435 // is setting 'value' to something.
436 if (value == -HUGE_VAL) {
437 value = 0.0;
438 rval = qm->evaluate( pd, handles[0], value, err );
439 CPPUNIT_ASSERT(!err && rval);
440 CPPUNIT_ASSERT( value != 0.0 );
441 }
442 std::vector<size_t> indices;
443 rval = qm->evaluate_with_indices( pd, handles[0], value, indices, err );
444 CPPUNIT_ASSERT(!MSQ_CHKERR(err)); CPPUNIT_ASSERT(rval);
445 std::vector<Vector3D> grad;
446 rval = qm->evaluate_with_gradient( pd, handles[0], value, indices, grad, err );
447 CPPUNIT_ASSERT(!MSQ_CHKERR(err)); CPPUNIT_ASSERT(rval);
448 std::vector<Matrix3D> hess;
449 rval = qm->evaluate_with_Hessian( pd, handles[0], value, indices, grad, hess, err );
450 CPPUNIT_ASSERT(!MSQ_CHKERR(err)); CPPUNIT_ASSERT(rval);
451 }
452
test_type_is_not_supported(EntityTopology type,QualityMetric * qm)453 void QualityMetricTester::test_type_is_not_supported( EntityTopology type, QualityMetric* qm )
454 {
455 MsqError err;
456 bool rval;
457 // create patch containing only elements of sepcified type
458 PatchData pd;
459 if (mSettings)
460 pd.attach_settings(mSettings);
461 get_ideal_element( type, false, pd );
462 // get list of evaluation locations in patch
463 std::vector<size_t> handles;
464 qm->get_evaluations( pd, handles, false, err );
465 if (err.error_code() == MsqError::UNSUPPORTED_ELEMENT)
466 return;
467 CPPUNIT_ASSERT(!err);
468 if (handles.empty())
469 return;
470 // evaluate the metric at one location
471 double value;
472 rval = qm->evaluate( pd, handles[0], value, err );
473 CPPUNIT_ASSERT_EQUAL(MsqError::UNSUPPORTED_ELEMENT, err.error_code());
474 std::vector<size_t> indices;
475 rval = qm->evaluate_with_indices( pd, handles[0], value, indices, err );
476 CPPUNIT_ASSERT_EQUAL(MsqError::UNSUPPORTED_ELEMENT, err.error_code());
477 std::vector<Vector3D> grad;
478 rval = qm->evaluate_with_gradient( pd, handles[0], value, indices, grad, err );
479 CPPUNIT_ASSERT_EQUAL(MsqError::UNSUPPORTED_ELEMENT, err.error_code());
480 std::vector<Matrix3D> hess;
481 rval = qm->evaluate_with_Hessian( pd, handles[0], value, indices, grad, hess, err );
482 CPPUNIT_ASSERT_EQUAL(MsqError::UNSUPPORTED_ELEMENT, err.error_code());
483 }
484
test_supported_element_types(QualityMetric * qm)485 void QualityMetricTester::test_supported_element_types( QualityMetric* qm )
486 {
487 for (int i = TRIANGLE; i < MIXED; ++i) {
488 if (i == POLYGON || i == POLYHEDRON || i == SEPTAHEDRON)
489 continue;
490 else if (type_is_supported( (EntityTopology)i ))
491 test_type_is_supported( (EntityTopology)i, qm );
492 else
493 test_type_is_not_supported( (EntityTopology)i, qm );
494 }
495 }
496
test_evaluate(PatchData & pd,QualityMetric * qm,double value)497 static void test_evaluate( PatchData& pd, QualityMetric* qm, double value )
498 {
499 MsqPrintError err( std::cout );
500 std::vector<size_t> handles;
501 qm->get_evaluations( pd, handles, false, err );
502 CPPUNIT_ASSERT(!MSQ_CHKERR(err));
503 CPPUNIT_ASSERT(!handles.empty());
504 for (unsigned i = 0; i < handles.size(); ++i) {
505 double qmval;
506 bool rval = qm->evaluate( pd, handles[i], qmval, err );
507 CPPUNIT_ASSERT(!MSQ_CHKERR(err));
508 CPPUNIT_ASSERT(rval);
509 CPPUNIT_ASSERT_DOUBLES_EQUAL( value, qmval, 1e-6 );
510 }
511 }
512
test_evaluate_unit_element(QualityMetric * qm,EntityTopology type,double value)513 void QualityMetricTester::test_evaluate_unit_element( QualityMetric* qm, EntityTopology type, double value )
514 {
515 PatchData pd;
516 if (mSettings)
517 pd.attach_settings(mSettings);
518 get_ideal_element( type, true, pd );
519 test_evaluate( pd, qm, value );
520 }
521
test_evaluate_unit_edge_element(QualityMetric * qm,EntityTopology type,double value)522 void QualityMetricTester::test_evaluate_unit_edge_element( QualityMetric* qm, EntityTopology type, double value )
523 {
524 PatchData pd;
525 if (mSettings)
526 pd.attach_settings(mSettings);
527 get_ideal_element( type, false, pd );
528 test_evaluate( pd, qm, value );
529 }
530
test_evaluate_unit_tris_about_vertex(QualityMetric * qm,double value)531 void QualityMetricTester::test_evaluate_unit_tris_about_vertex( QualityMetric* qm, double value )
532 {
533 PatchData pd;
534 if (mSettings)
535 pd.attach_settings(mSettings);
536 get_ideal_tris( pd, true );
537 test_evaluate( pd, qm, value );
538 }
539
test_evaluate_unit_quads_about_vertex(QualityMetric * qm,double value)540 void QualityMetricTester::test_evaluate_unit_quads_about_vertex( QualityMetric* qm, double value )
541 {
542 PatchData pd;
543 if (mSettings)
544 pd.attach_settings(mSettings);
545 get_ideal_quads( pd );
546 test_evaluate( pd, qm, value );
547 }
548
test_evaluate_unit_hexes_about_vertex(QualityMetric * qm,double value)549 void QualityMetricTester::test_evaluate_unit_hexes_about_vertex( QualityMetric* qm, double value )
550 {
551 PatchData pd;
552 if (mSettings)
553 pd.attach_settings(mSettings);
554 get_ideal_hexes( pd );
555 test_evaluate( pd, qm, value );
556 }
557
test_evaluate_unit_edge_tris_about_vertex(QualityMetric * qm,double value)558 void QualityMetricTester::test_evaluate_unit_edge_tris_about_vertex( QualityMetric* qm, double value )
559 {
560 PatchData pd;
561 if (mSettings)
562 pd.attach_settings(mSettings);
563 get_ideal_tris( pd, false );
564 test_evaluate( pd, qm, value );
565 }
566
test_get_element_evaluations(QualityMetric * qm)567 void QualityMetricTester::test_get_element_evaluations( QualityMetric* qm )
568 {
569 MsqPrintError err( std::cout );
570 PatchData pd;
571 if (mSettings)
572 pd.attach_settings(mSettings);
573 unsigned count = 0;
574 std::vector<size_t> handles;
575
576 if (type_is_supported(HEXAHEDRON)) {
577 get_ideal_hexes( pd );
578 qm->get_evaluations( pd, handles, false, err );
579 CPPUNIT_ASSERT(!MSQ_CHKERR(err));
580 CPPUNIT_ASSERT_EQUAL( pd.num_elements(), handles.size() );
581 std::sort( handles.begin(), handles.end() );
582 handles.erase( std::unique( handles.begin(), handles.end() ), handles.end() );
583 CPPUNIT_ASSERT_EQUAL( pd.num_elements(), handles.size() );
584 ++count;
585 }
586
587 if (type_is_supported(QUADRILATERAL)) {
588 get_ideal_quads( pd );
589 qm->get_evaluations( pd, handles, false, err );
590 CPPUNIT_ASSERT(!MSQ_CHKERR(err));
591 CPPUNIT_ASSERT_EQUAL( pd.num_elements(), handles.size() );
592 std::sort( handles.begin(), handles.end() );
593 handles.erase( std::unique( handles.begin(), handles.end() ), handles.end() );
594 CPPUNIT_ASSERT_EQUAL( pd.num_elements(), handles.size() );
595 ++count;
596 }
597
598 if (type_is_supported(TRIANGLE)) {
599 get_ideal_tris( pd, false );
600 qm->get_evaluations( pd, handles, false, err );
601 CPPUNIT_ASSERT(!MSQ_CHKERR(err));
602 CPPUNIT_ASSERT_EQUAL( pd.num_elements(), handles.size() );
603 std::sort( handles.begin(), handles.end() );
604 handles.erase( std::unique( handles.begin(), handles.end() ), handles.end() );
605 CPPUNIT_ASSERT_EQUAL( pd.num_elements(), handles.size() );
606 ++count;
607 }
608
609 CPPUNIT_ASSERT(count > 0);
610 }
611
test_get_vertex_evaluations(QualityMetric * qm)612 void QualityMetricTester::test_get_vertex_evaluations( QualityMetric* qm )
613 {
614 MsqPrintError err( std::cout );
615 PatchData pd;
616 if (mSettings)
617 pd.attach_settings(mSettings);
618 unsigned count = 0;
619 std::vector<size_t> handles;
620
621 if (type_is_supported(HEXAHEDRON)) {
622 get_ideal_hexes( pd );
623 qm->get_evaluations( pd, handles, false, err );
624 CPPUNIT_ASSERT(!MSQ_CHKERR(err));
625 CPPUNIT_ASSERT_EQUAL( (size_t)1, handles.size() );
626 ++count;
627 }
628
629 if (type_is_supported(TRIANGLE)) {
630 get_ideal_tris( pd, false );
631 qm->get_evaluations( pd, handles, false, err );
632 CPPUNIT_ASSERT(!MSQ_CHKERR(err));
633 CPPUNIT_ASSERT_EQUAL( (size_t)1, handles.size() );
634 ++count;
635 }
636
637 if (type_is_supported(QUADRILATERAL)) {
638 get_ideal_quads( pd );
639 qm->get_evaluations( pd, handles, false, err );
640 CPPUNIT_ASSERT(!MSQ_CHKERR(err));
641 CPPUNIT_ASSERT_EQUAL( (size_t)1, handles.size() );
642 ++count;
643 }
644
645 CPPUNIT_ASSERT(count > 0);
646 }
647
test_get_sample_evaluations(QualityMetric * qm)648 void QualityMetricTester::test_get_sample_evaluations( QualityMetric* qm )
649 {
650 MsqPrintError err( std::cout );
651 PatchData pd;
652 if (mSettings)
653 pd.attach_settings(mSettings);
654 unsigned count = 0;
655 std::vector<size_t> handles;
656
657 if (type_is_supported(HEXAHEDRON)) {
658 get_ideal_hexes( pd );
659 size_t expected_evals = 0;
660 for (size_t i= 0; i < pd.num_elements(); ++i)
661 expected_evals += pd.get_samples(i).num_nodes();
662
663 qm->get_evaluations( pd, handles, false, err );
664 CPPUNIT_ASSERT(!MSQ_CHKERR(err));
665 CPPUNIT_ASSERT_EQUAL( expected_evals, handles.size() );
666
667 std::sort( handles.begin(), handles.end() );
668 handles.erase( std::unique( handles.begin(), handles.end() ), handles.end() );
669 CPPUNIT_ASSERT_EQUAL(expected_evals, handles.size() );
670
671 ++count;
672 }
673
674 if (type_is_supported(TRIANGLE)) {
675 get_ideal_tris( pd, false );
676 size_t expected_evals = 0;
677 for (size_t i= 0; i < pd.num_elements(); ++i)
678 expected_evals += pd.get_samples(i).num_nodes();
679
680 qm->get_evaluations( pd, handles, false, err );
681 CPPUNIT_ASSERT(!MSQ_CHKERR(err));
682 CPPUNIT_ASSERT_EQUAL( expected_evals, handles.size() );
683
684 std::sort( handles.begin(), handles.end() );
685 handles.erase( std::unique( handles.begin(), handles.end() ), handles.end() );
686 CPPUNIT_ASSERT_EQUAL(expected_evals, handles.size() );
687
688 ++count;
689 }
690
691 if (type_is_supported(QUADRILATERAL)) {
692 get_ideal_quads( pd );
693 size_t expected_evals = 0;
694 for (size_t i= 0; i < pd.num_elements(); ++i)
695 expected_evals += pd.get_samples(i).num_nodes();
696
697 qm->get_evaluations( pd, handles, false, err );
698 CPPUNIT_ASSERT(!MSQ_CHKERR(err));
699 CPPUNIT_ASSERT_EQUAL( expected_evals, handles.size() );
700
701 std::sort( handles.begin(), handles.end() );
702 handles.erase( std::unique( handles.begin(), handles.end() ), handles.end() );
703 CPPUNIT_ASSERT_EQUAL(expected_evals, handles.size() );
704
705 ++count;
706 }
707
708 CPPUNIT_ASSERT(count > 0);
709 }
710 void get_ideal_element( EntityTopology type,
711 bool unit_area,
712 PatchData& pd,
713 int free_vertex_index );
714
test_get_edge_evaluations(EdgeQM * qm)715 void QualityMetricTester::test_get_edge_evaluations( EdgeQM* qm )
716 {
717 MsqPrintError err( std::cout );
718 PatchData pd;
719 if (mSettings)
720 pd.attach_settings(mSettings);
721 std::vector<size_t> handles;
722
723 CPPUNIT_ASSERT(!types.empty());
724 for (size_t i = 0; i < types.size(); ++i) {
725 get_ideal_element( types[i], true, pd );
726 handles.clear();
727 qm->get_evaluations( pd, handles, false, err );
728 ASSERT_NO_ERROR(err);
729 CPPUNIT_ASSERT_EQUAL( TopologyInfo::edges(types[i]), (unsigned)handles.size() );
730 }
731 }
732
test_get_in_element_evaluations(ElemSampleQM * qm)733 void QualityMetricTester::test_get_in_element_evaluations( ElemSampleQM* qm )
734 {
735 MsqPrintError err( std::cout );
736 PatchData pd;
737 if (mSettings)
738 pd.attach_settings(mSettings);
739 std::vector<size_t> handles1, handles2;
740
741 CPPUNIT_ASSERT(!types.empty());
742
743 get_ideal_element( types[0], false, pd );
744 CPPUNIT_ASSERT(pd.num_elements() == 1);
745 qm->get_evaluations( pd, handles1, false, err );
746 CPPUNIT_ASSERT(!MSQ_CHKERR(err));
747 CPPUNIT_ASSERT(!handles1.empty());
748 qm->get_element_evaluations( pd, 0, handles2, err );
749 CPPUNIT_ASSERT(!MSQ_CHKERR(err));
750 CPPUNIT_ASSERT(!handles2.empty());
751 std::sort( handles1.begin(), handles1.end() );
752 std::sort( handles2.begin(), handles2.end() );
753 CPPUNIT_ASSERT( handles1 == handles2 );
754
755 get_ideal_element( types.back(), false, pd );
756 CPPUNIT_ASSERT(pd.num_elements() == 1);
757 qm->get_evaluations( pd, handles1, false, err );
758 CPPUNIT_ASSERT(!MSQ_CHKERR(err));
759 CPPUNIT_ASSERT(!handles1.empty());
760 qm->get_element_evaluations( pd, 0, handles2, err );
761 CPPUNIT_ASSERT(!MSQ_CHKERR(err));
762 CPPUNIT_ASSERT(!handles2.empty());
763 std::sort( handles1.begin(), handles1.end() );
764 std::sort( handles2.begin(), handles2.end() );
765 CPPUNIT_ASSERT( handles1 == handles2 );
766 }
767
test_get_indices_fixed(QualityMetric * qm)768 void QualityMetricTester::test_get_indices_fixed( QualityMetric* qm )
769 {
770 MsqPrintError err( std::cout );
771 PatchData pd1, pd2;
772 if (mSettings) {
773 pd1.attach_settings(mSettings);
774 pd2.attach_settings(mSettings);
775 }
776 std::vector<size_t> handles1, handles2, indices1, indices2;
777 double qm_val1, qm_val2;
778
779 CPPUNIT_ASSERT(!types.empty());
780 // get element with no fixed vertices
781 get_ideal_element( types.back(), false, pd1, false );
782 qm->get_evaluations( pd1, handles1, false, err );
783 CPPUNIT_ASSERT(!MSQ_CHKERR(err));
784 CPPUNIT_ASSERT(!handles1.empty());
785 // call evaluate with indices
786 size_t handle = handles1.back();
787 qm->evaluate_with_indices( pd1, handle, qm_val1, indices1, err );
788 CPPUNIT_ASSERT(!MSQ_CHKERR(err));
789 CPPUNIT_ASSERT(!indices1.empty());
790 // get element with one fixed vertex
791 get_ideal_element( types.back(), false, pd2, true );
792 qm->get_evaluations( pd2, handles2, false, err );
793 CPPUNIT_ASSERT(!MSQ_CHKERR(err));
794 CPPUNIT_ASSERT(!handles2.empty());
795 // If vertex-based metric, need to find corresponding handle
796 // in second patch. The vertex order changed as a result of
797 // one of the vertices beging fixed in the second patch.
798 if (qm->get_metric_type() == QualityMetric::VERTEX_BASED )
799 {
800 handle = handles1.back();
801 const size_t* conn1 = pd1.element_by_index(0).get_vertex_index_array();
802 const size_t* conn2 = pd2.element_by_index(0).get_vertex_index_array();
803 const size_t len = TopologyInfo::corners( types.back() );
804 const size_t* ptr = std::find( conn1, conn1+len, handle );
805 handle = conn2[ptr - conn1];
806 }
807 // call evaluate with indices
808 qm->evaluate_with_indices( pd2, handle, qm_val2, indices2, err );
809 CPPUNIT_ASSERT(!MSQ_CHKERR(err));
810 CPPUNIT_ASSERT(!indices2.empty());
811 // make sure we got the same QM value
812 // (shoudn't be affected by fixed/non-fixed)
813 CPPUNIT_ASSERT_DOUBLES_EQUAL( qm_val1, qm_val2, 1e-6 );
814 // should have gotten back one less index (for element-based metrics)
815 if (handles1.size() == 1) {
816 CPPUNIT_ASSERT_EQUAL( indices1.size(), indices2.size()+1 );
817 }
818 // indices2 shouldn't contain any fixed vertices
819 std::sort( indices2.begin(), indices2.end() );
820 CPPUNIT_ASSERT( indices2.back() < pd2.num_free_vertices() );
821 // indices2 shouldn/t contain any duplicates
822 CPPUNIT_ASSERT( std::unique(indices2.begin(), indices2.end()) == indices2.end() );
823 }
824
test_get_element_indices(QualityMetric * qm)825 void QualityMetricTester::test_get_element_indices( QualityMetric* qm )
826 {
827 MsqPrintError err( std::cout );
828 PatchData pd;
829 if (mSettings)
830 pd.attach_settings(mSettings);
831 std::vector<size_t> handles, indices, verts;
832 double qm_val;
833
834 CPPUNIT_ASSERT(qm->get_metric_type() == QualityMetric::ELEMENT_BASED);
835 CPPUNIT_ASSERT( !types.empty() );
836
837 for (size_t i = 0; i < types.size(); ++i) {
838 // construct patch w/ one fixed vertex
839 get_ideal_element( types[i], false, pd, true );
840 qm->get_evaluations( pd, handles, false, err );
841 CPPUNIT_ASSERT( !MSQ_CHKERR(err) );
842 CPPUNIT_ASSERT_EQUAL( (size_t)1, handles.size() );
843 // get vertex list
844 verts.clear();
845 pd.element_by_index(0).get_vertex_indices( verts );
846 // remove fixed vertex
847 for (size_t j = 0; j <verts.size(); ++j)
848 if (verts[j] >= pd.num_free_vertices()) {
849 verts.erase( verts.begin() + j);
850 break;
851 }
852 // evaluate metric
853 qm->evaluate_with_indices( pd, handles[0], qm_val, indices, err );
854 CPPUNIT_ASSERT( !MSQ_CHKERR(err) );
855 // index list should match vertex list (element metric
856 // should depend on all free vertices in element).
857 CPPUNIT_ASSERT_EQUAL( verts.size(), indices.size() );
858 std::sort( verts.begin(), verts.end() );
859 std::sort( indices.begin(), indices.end() );
860 CPPUNIT_ASSERT( verts == indices );
861 }
862 }
863
test_get_vertex_indices(QualityMetric * qm)864 void QualityMetricTester::test_get_vertex_indices( QualityMetric* qm )
865 {
866 MsqPrintError err( std::cout );
867 PatchData pd;
868 if (mSettings)
869 pd.attach_settings(mSettings);
870 std::vector<size_t> indices, verts;
871 double qm_val;
872
873 CPPUNIT_ASSERT(qm->get_metric_type() == QualityMetric::VERTEX_BASED);
874 CPPUNIT_ASSERT( !types.empty() );
875
876 for (size_t i = 0; i < types.size(); ++i) {
877 // construct patch w/ one fixed vertex
878 get_ideal_element( types[i], false, pd, true );
879 // get adjacent vertices
880 verts.clear();
881 pd.get_adjacent_vertex_indices( 0, verts, err );
882 CPPUNIT_ASSERT( !MSQ_CHKERR(err) );
883 // remove fixed vertex
884 for (size_t j = 0; j <verts.size(); ++j)
885 if (verts[j] >= pd.num_free_vertices()) {
886 verts.erase( verts.begin() + j);
887 break;
888 }
889 verts.push_back( 0 );
890 // evaluate metric
891 qm->evaluate_with_indices( pd, 0, qm_val, indices, err );
892 CPPUNIT_ASSERT( !MSQ_CHKERR(err) );
893 // vertex metric should depend on all adjacent free vertices
894 CPPUNIT_ASSERT_EQUAL( verts.size(), indices.size() );
895 std::sort( verts.begin(), verts.end() );
896 std::sort( indices.begin(), indices.end() );
897 CPPUNIT_ASSERT( verts == indices );
898 }
899 }
900
test_get_edge_indices(EdgeQM * qm)901 void QualityMetricTester::test_get_edge_indices( EdgeQM* qm )
902 {
903 MsqPrintError err( std::cout );
904 PatchData pd, pd2;
905 if (mSettings) {
906 pd.attach_settings(mSettings);
907 pd2.attach_settings(mSettings);
908 }
909 std::vector<size_t> handles, indices, indices2, verts;
910 std::vector<size_t>::iterator it;
911 double qm_val;
912
913 CPPUNIT_ASSERT(qm->get_metric_type() == QualityMetric::VERTEX_BASED);
914 CPPUNIT_ASSERT( !types.empty() );
915
916 for (size_t i = 0; i < types.size(); ++i) {
917 // construct patch w/ no free vertices
918 get_ideal_element( types[i], false, pd, false );
919
920 qm->get_evaluations( pd, handles, false, err );
921 ASSERT_NO_ERROR(err);
922 CPPUNIT_ASSERT_EQUAL( (size_t)TopologyInfo::edges(types[i]), handles.size() );
923
924 for (size_t j = 0; j < handles.size(); ++j) {
925 // evaluate metric
926 qm->evaluate_with_indices( pd, handles[j], qm_val, indices, err );
927 ASSERT_NO_ERROR(err);
928 // evaluation at each edge should depend on at least the two end vertices
929 CPPUNIT_ASSERT( indices.size() >= (size_t)2 );
930 }
931 }
932
933 for (size_t i = 0; i < types.size(); ++i) {
934 // construct patch w/ one free vertex
935 get_ideal_element( types[i], false, pd, true );
936 const size_t fixed_vertex = pd.num_free_vertices();
937
938 qm->get_evaluations( pd, handles, false, err );
939 ASSERT_NO_ERROR(err);
940 CPPUNIT_ASSERT_EQUAL( (size_t)TopologyInfo::edges(types[i]), handles.size() );
941
942 for (size_t j = 0; j < handles.size(); ++j) {
943 // evaluate metric
944 qm->evaluate_with_indices( pd, handles[j], qm_val, indices, err );
945 ASSERT_NO_ERROR(err);
946 // evaluation at each edge should depend on at least the two end vertices
947 CPPUNIT_ASSERT( !indices.empty() );
948
949 // indices should never contain the index of a fixed vertex
950 it = std::find(indices.begin(), indices.end(), fixed_vertex );
951 CPPUNIT_ASSERT( it == indices.end() );
952 }
953 }
954 }
955
test_get_sample_indices(QualityMetric * qm)956 void QualityMetricTester::test_get_sample_indices( QualityMetric* qm )
957 {
958 MsqPrintError err( std::cout );
959 PatchData pd;
960 if (mSettings)
961 pd.attach_settings(mSettings);
962 std::vector<size_t> handles, indices, verts;
963 double qm_val;
964
965 CPPUNIT_ASSERT(qm->get_metric_type() == QualityMetric::ELEMENT_BASED);
966 CPPUNIT_ASSERT( !types.empty() );
967
968 for (size_t i = 0; i < types.size(); ++i) {
969 // get evaluation locations
970 get_ideal_element( types[i], false, pd, true );
971 const size_t count = pd.get_samples(0).num_nodes();
972 CPPUNIT_ASSERT( count > 0 );
973 qm->get_evaluations( pd, handles, false, err );
974 CPPUNIT_ASSERT( !MSQ_CHKERR(err) );
975 // should have one evaluation for each sample point
976 CPPUNIT_ASSERT_EQUAL( count, handles.size() );
977 // get vertex list
978 verts.clear();
979 pd.element_by_index(0).get_vertex_indices( verts );
980 // remove fixed vertex
981 for (size_t j = 0; j <verts.size(); ++j)
982 if (verts[j] >= pd.num_free_vertices()) {
983 verts.erase( verts.begin() + j);
984 break;
985 }
986 // evaluate metric
987 qm->evaluate_with_indices( pd, handles[0], qm_val, indices, err );
988 CPPUNIT_ASSERT( !MSQ_CHKERR(err) );
989 // only one fixed vertex, and presumably every possible evaluation
990 // depends on more than one vertex, so should be at least one
991 // other, non-fixed, vertex.
992 CPPUNIT_ASSERT( indices.size() > 0 );
993 // make sure no duplicates
994 std::sort( indices.begin(), indices.end() );
995 CPPUNIT_ASSERT( std::unique(indices.begin(), indices.end()) == indices.end() );
996 // check that all indices are in vertex list
997 CPPUNIT_ASSERT( indices.size() <= verts.size() );
998 for (std::vector<size_t>::iterator j = indices.begin(); j != indices.end(); ++j)
999 CPPUNIT_ASSERT( std::find( verts.begin(), verts.end(), *j ) != verts.end() );
1000 }
1001 }
1002
compare_eval_and_eval_with_indices(QualityMetric * qm)1003 void QualityMetricTester::compare_eval_and_eval_with_indices( QualityMetric* qm )
1004 {
1005 MsqPrintError err( std::cout );
1006 PatchData pd;
1007 if (mSettings)
1008 pd.attach_settings(mSettings);
1009 std::vector<size_t> handles, indices;
1010 double qm_val1, qm_val2;
1011 bool rval;
1012
1013 CPPUNIT_ASSERT( !types.empty() );
1014 for (size_t i = 0; i < types.size(); ++i) {
1015 get_ideal_element( types[i], false, pd );
1016 qm->get_evaluations( pd, handles, false, err );
1017 CPPUNIT_ASSERT( !MSQ_CHKERR(err) );
1018 CPPUNIT_ASSERT( !handles.empty() );
1019 for (size_t j = 0; j < handles.size(); ++j) {
1020 rval = qm->evaluate( pd, handles[j], qm_val1, err );
1021 CPPUNIT_ASSERT( !MSQ_CHKERR(err) );
1022 CPPUNIT_ASSERT( rval );
1023 rval = qm->evaluate_with_indices( pd, handles[j], qm_val2, indices, err );
1024 CPPUNIT_ASSERT( !MSQ_CHKERR(err) );
1025 CPPUNIT_ASSERT( rval );
1026
1027 CPPUNIT_ASSERT_DOUBLES_EQUAL( qm_val1, qm_val2, 1e-6 );
1028 }
1029 }
1030 }
1031
compare_eval_with_indices_and_eval_with_gradient(QualityMetric * qm,PatchData & pd)1032 void QualityMetricTester::compare_eval_with_indices_and_eval_with_gradient( QualityMetric* qm, PatchData& pd )
1033 {
1034 MsqPrintError err( std::cout );
1035 std::vector<size_t> handles, indices1, indices2;
1036 std::vector<Vector3D> grad;
1037 double qm_val1, qm_val2;
1038 bool rval;
1039
1040 qm->get_evaluations( pd, handles, false, err );
1041 CPPUNIT_ASSERT( !MSQ_CHKERR(err) );
1042 CPPUNIT_ASSERT( !handles.empty() );
1043 for (size_t j = 0; j < handles.size(); ++j) {
1044 rval = qm->evaluate_with_indices( pd, handles[j], qm_val1, indices1, err );
1045 CPPUNIT_ASSERT( !MSQ_CHKERR(err) );
1046 CPPUNIT_ASSERT( rval );
1047 rval = qm->evaluate_with_gradient( pd, handles[j], qm_val2, indices2, grad, err );
1048 CPPUNIT_ASSERT( !MSQ_CHKERR(err) );
1049 CPPUNIT_ASSERT( rval );
1050
1051 CPPUNIT_ASSERT_DOUBLES_EQUAL( qm_val1, qm_val2, 1e-6 );
1052 CPPUNIT_ASSERT_EQUAL( indices1.size(), indices2.size() );
1053
1054 std::sort( indices1.begin(), indices1.end() );
1055 std::sort( indices2.begin(), indices2.end() );
1056 CPPUNIT_ASSERT( indices1 == indices2 );
1057 CPPUNIT_ASSERT( !indices1.empty() );
1058 }
1059 }
1060
compare_eval_with_indices_and_eval_with_gradient(QualityMetric * qm)1061 void QualityMetricTester::compare_eval_with_indices_and_eval_with_gradient( QualityMetric* qm )
1062 {
1063 PatchData pd;
1064 if (mSettings)
1065 pd.attach_settings(mSettings);
1066
1067 CPPUNIT_ASSERT( !types.empty() );
1068 for (size_t i = 0; i < types.size(); ++i) {
1069 get_ideal_element( types[i], false, pd, true );
1070 compare_eval_with_indices_and_eval_with_gradient( qm, pd );
1071 }
1072 }
1073
compare_eval_with_indices_and_eval_with_hessian(QualityMetric * qm,PatchData & pd)1074 void QualityMetricTester::compare_eval_with_indices_and_eval_with_hessian(
1075 QualityMetric* qm, PatchData& pd )
1076 {
1077 MsqPrintError err( std::cout );
1078 std::vector<size_t> handles, indices1, indices2;
1079 std::vector<Vector3D> grad;
1080 std::vector<Matrix3D> Hess;
1081 double qm_val1, qm_val2;
1082 bool rval;
1083
1084 qm->get_evaluations( pd, handles, false, err );
1085 CPPUNIT_ASSERT( !MSQ_CHKERR(err) );
1086 CPPUNIT_ASSERT( !handles.empty() );
1087 for (size_t j = 0; j < handles.size(); ++j) {
1088 rval = qm->evaluate_with_indices( pd, handles[j], qm_val1, indices1, err );
1089 CPPUNIT_ASSERT( !MSQ_CHKERR(err) );
1090 CPPUNIT_ASSERT( rval );
1091 rval = qm->evaluate_with_Hessian( pd, handles[j], qm_val2, indices2, grad, Hess, err );
1092 CPPUNIT_ASSERT( !MSQ_CHKERR(err) );
1093 CPPUNIT_ASSERT( rval );
1094
1095 CPPUNIT_ASSERT_DOUBLES_EQUAL( qm_val1, qm_val2, 1e-6 );
1096 CPPUNIT_ASSERT_EQUAL( indices1.size(), indices2.size() );
1097
1098 std::sort( indices1.begin(), indices1.end() );
1099 std::sort( indices2.begin(), indices2.end() );
1100 CPPUNIT_ASSERT( indices1 == indices2 );
1101 CPPUNIT_ASSERT( !indices1.empty() );
1102 }
1103 }
1104
compare_eval_with_indices_and_eval_with_hessian(QualityMetric * qm)1105 void QualityMetricTester::compare_eval_with_indices_and_eval_with_hessian( QualityMetric* qm )
1106 {
1107 PatchData pd;
1108 if (mSettings)
1109 pd.attach_settings(mSettings);
1110
1111 CPPUNIT_ASSERT( !types.empty() );
1112 for (size_t i = 0; i < types.size(); ++i) {
1113 get_nonideal_element( types[i], pd, true );
1114 compare_eval_with_indices_and_eval_with_hessian( qm, pd );
1115 }
1116 }
1117
compare_eval_with_indices_and_eval_with_diagonal(QualityMetric * qm,PatchData & pd)1118 void QualityMetricTester::compare_eval_with_indices_and_eval_with_diagonal(
1119 QualityMetric* qm, PatchData& pd )
1120 {
1121 MsqPrintError err( std::cout );
1122 std::vector<size_t> handles, indices1, indices2;
1123 std::vector<Vector3D> grad;
1124 std::vector<SymMatrix3D> Hess;
1125 double qm_val1, qm_val2;
1126 bool rval;
1127
1128 qm->get_evaluations( pd, handles, false, err );
1129 CPPUNIT_ASSERT( !MSQ_CHKERR(err) );
1130 CPPUNIT_ASSERT( !handles.empty() );
1131 for (size_t j = 0; j < handles.size(); ++j) {
1132 rval = qm->evaluate_with_indices( pd, handles[j], qm_val1, indices1, err );
1133 CPPUNIT_ASSERT( !MSQ_CHKERR(err) );
1134 CPPUNIT_ASSERT( rval );
1135 rval = qm->evaluate_with_Hessian_diagonal( pd, handles[j], qm_val2, indices2, grad, Hess, err );
1136 CPPUNIT_ASSERT( !MSQ_CHKERR(err) );
1137 CPPUNIT_ASSERT( rval );
1138
1139 CPPUNIT_ASSERT_DOUBLES_EQUAL( qm_val1, qm_val2, 1e-6 );
1140 CPPUNIT_ASSERT_EQUAL( indices1.size(), indices2.size() );
1141
1142 std::sort( indices1.begin(), indices1.end() );
1143 std::sort( indices2.begin(), indices2.end() );
1144 CPPUNIT_ASSERT( indices1 == indices2 );
1145 CPPUNIT_ASSERT( !indices1.empty() );
1146 }
1147 }
1148
compare_eval_with_indices_and_eval_with_diagonal(QualityMetric * qm)1149 void QualityMetricTester::compare_eval_with_indices_and_eval_with_diagonal( QualityMetric* qm )
1150 {
1151 PatchData pd;
1152 if (mSettings)
1153 pd.attach_settings(mSettings);
1154
1155 CPPUNIT_ASSERT( !types.empty() );
1156 for (size_t i = 0; i < types.size(); ++i) {
1157 get_nonideal_element( types[i], pd, true );
1158 compare_eval_with_indices_and_eval_with_diagonal( qm, pd );
1159 }
1160 }
1161
compare_eval_with_grad_and_eval_with_hessian(QualityMetric * qm,PatchData & pd)1162 void QualityMetricTester::compare_eval_with_grad_and_eval_with_hessian(
1163 QualityMetric* qm, PatchData& pd )
1164 {
1165 MsqPrintError err( std::cout );
1166 std::vector<size_t> handles, indices1, indices2;
1167 std::vector<Vector3D> grad1, grad2;
1168 std::vector<Matrix3D> Hess;
1169 double qm_val1, qm_val2;
1170 bool rval;
1171
1172 qm->get_evaluations( pd, handles, false, err );
1173 CPPUNIT_ASSERT( !MSQ_CHKERR(err) );
1174 CPPUNIT_ASSERT( !handles.empty() );
1175 for (size_t j = 0; j < handles.size(); ++j) {
1176 rval = qm->evaluate_with_gradient( pd, handles[j], qm_val1, indices1, grad1, err );
1177 CPPUNIT_ASSERT( !MSQ_CHKERR(err) );
1178 CPPUNIT_ASSERT( rval );
1179 rval = qm->evaluate_with_Hessian( pd, handles[j], qm_val2, indices2, grad2, Hess, err );
1180 CPPUNIT_ASSERT( !MSQ_CHKERR(err) );
1181 CPPUNIT_ASSERT( rval );
1182
1183 CPPUNIT_ASSERT_DOUBLES_EQUAL( qm_val1, qm_val2, 1e-6 );
1184 CPPUNIT_ASSERT_EQUAL( indices1.size(), indices2.size() );
1185
1186 for (size_t k = 0; k < indices1.size(); ++k) {
1187 std::vector<size_t>::iterator it =
1188 std::find( indices2.begin(), indices2.end(), indices1[k] );
1189 CPPUNIT_ASSERT( it != indices2.end() );
1190 size_t m = it - indices2.begin();
1191 CPPUNIT_ASSERT_VECTORS_EQUAL( grad1[k], grad2[m], 1e-6 );
1192 }
1193 }
1194 }
1195
compare_eval_with_grad_and_eval_with_hessian(QualityMetric * qm)1196 void QualityMetricTester::compare_eval_with_grad_and_eval_with_hessian( QualityMetric* qm )
1197 {
1198 PatchData pd;
1199 if (mSettings)
1200 pd.attach_settings(mSettings);
1201
1202 CPPUNIT_ASSERT( !types.empty() );
1203 for (size_t i = 0; i < types.size(); ++i) {
1204 get_nonideal_element( types[i], pd, true );
1205 compare_eval_with_grad_and_eval_with_hessian( qm, pd );
1206 }
1207 }
1208
compare_eval_with_grad_and_eval_with_diagonal(QualityMetric * qm,PatchData & pd)1209 void QualityMetricTester::compare_eval_with_grad_and_eval_with_diagonal(
1210 QualityMetric* qm, PatchData& pd )
1211 {
1212 MsqPrintError err( std::cout );
1213 std::vector<size_t> handles, indices1, indices2;
1214 std::vector<Vector3D> grad1, grad2;
1215 std::vector<SymMatrix3D> Hess;
1216 double qm_val1, qm_val2;
1217 bool rval;
1218
1219 qm->get_evaluations( pd, handles, false, err );
1220 CPPUNIT_ASSERT( !MSQ_CHKERR(err) );
1221 CPPUNIT_ASSERT( !handles.empty() );
1222 for (size_t j = 0; j < handles.size(); ++j) {
1223 rval = qm->evaluate_with_gradient( pd, handles[j], qm_val1, indices1, grad1, err );
1224 CPPUNIT_ASSERT( !MSQ_CHKERR(err) );
1225 CPPUNIT_ASSERT( rval );
1226 rval = qm->evaluate_with_Hessian_diagonal( pd, handles[j], qm_val2, indices2, grad2, Hess, err );
1227 CPPUNIT_ASSERT( !MSQ_CHKERR(err) );
1228 CPPUNIT_ASSERT( rval );
1229
1230 CPPUNIT_ASSERT_DOUBLES_EQUAL( qm_val1, qm_val2, 1e-6 );
1231 CPPUNIT_ASSERT_EQUAL( indices1.size(), indices2.size() );
1232
1233 for (size_t k = 0; k < indices1.size(); ++k) {
1234 std::vector<size_t>::iterator it =
1235 std::find( indices2.begin(), indices2.end(), indices1[k] );
1236 CPPUNIT_ASSERT( it != indices2.end() );
1237 size_t m = it - indices2.begin();
1238 CPPUNIT_ASSERT_VECTORS_EQUAL( grad1[k], grad2[m], 1e-6 );
1239 }
1240 }
1241 }
1242
compare_eval_with_grad_and_eval_with_diagonal(QualityMetric * qm)1243 void QualityMetricTester::compare_eval_with_grad_and_eval_with_diagonal( QualityMetric* qm )
1244 {
1245 PatchData pd;
1246 if (mSettings)
1247 pd.attach_settings(mSettings);
1248
1249 CPPUNIT_ASSERT( !types.empty() );
1250 for (size_t i = 0; i < types.size(); ++i) {
1251 get_nonideal_element( types[i], pd, true );
1252 compare_eval_with_grad_and_eval_with_diagonal( qm, pd );
1253 }
1254 }
1255
compare_eval_with_diag_and_eval_with_hessian(QualityMetric * qm,PatchData & pd)1256 void QualityMetricTester::compare_eval_with_diag_and_eval_with_hessian(
1257 QualityMetric* qm, PatchData& pd )
1258 {
1259 MsqPrintError err( std::cout );
1260 std::vector<size_t> handles, indices1, indices2;
1261 std::vector<Vector3D> grad1, grad2;
1262 std::vector<SymMatrix3D> hess1;
1263 std::vector<Matrix3D> hess2;
1264 double qm_val1, qm_val2;
1265 bool rval;
1266
1267 qm->get_evaluations( pd, handles, false, err );
1268 CPPUNIT_ASSERT( !MSQ_CHKERR(err) );
1269 CPPUNIT_ASSERT( !handles.empty() );
1270 for (size_t j = 0; j < handles.size(); ++j) {
1271 rval = qm->evaluate_with_Hessian_diagonal( pd, handles[j], qm_val1, indices1, grad1, hess1, err );
1272 CPPUNIT_ASSERT( !MSQ_CHKERR(err) );
1273 CPPUNIT_ASSERT( rval );
1274 rval = qm->evaluate_with_Hessian( pd, handles[j], qm_val2, indices2, grad2, hess2, err );
1275 CPPUNIT_ASSERT( !MSQ_CHKERR(err) );
1276 CPPUNIT_ASSERT( rval );
1277
1278 CPPUNIT_ASSERT_DOUBLES_EQUAL( qm_val1, qm_val2, 1e-6 );
1279
1280 CPPUNIT_ASSERT_EQUAL( indices1.size(), indices2.size() );
1281 CPPUNIT_ASSERT_EQUAL( grad1.size(), grad2.size() );
1282 CPPUNIT_ASSERT_EQUAL( hess2.size(), hess1.size() * (hess1.size()+1) / 2 );
1283 unsigned h2step = indices2.size();
1284 std::vector<Matrix3D>::const_iterator h2i = hess2.begin();
1285
1286 for (size_t k = 0; k < indices2.size(); ++k) {
1287 std::vector<size_t>::iterator it =
1288 std::find( indices1.begin(), indices1.end(), indices2[k] );
1289 CPPUNIT_ASSERT( it != indices1.end() );
1290 size_t m = it - indices1.begin();
1291 CPPUNIT_ASSERT_VECTORS_EQUAL( grad1[m], grad2[k], 5e-5 );
1292 CPPUNIT_ASSERT_MATRICES_EQUAL( Matrix3D(hess1[m]), *h2i, 5e-5 );
1293 h2i += h2step--;
1294 }
1295 }
1296 }
1297
compare_eval_with_diag_and_eval_with_hessian(QualityMetric * qm)1298 void QualityMetricTester::compare_eval_with_diag_and_eval_with_hessian( QualityMetric* qm )
1299 {
1300 PatchData pd;
1301 if (mSettings)
1302 pd.attach_settings(mSettings);
1303
1304 CPPUNIT_ASSERT( !types.empty() );
1305 for (size_t i = 0; i < types.size(); ++i) {
1306 get_nonideal_element( types[i], pd, true );
1307 compare_eval_with_diag_and_eval_with_hessian( qm, pd );
1308 }
1309 }
1310
compare_analytical_and_numerical_gradients(QualityMetric * qm,PatchData & pd)1311 void QualityMetricTester::compare_analytical_and_numerical_gradients(
1312 QualityMetric* qm, PatchData& pd )
1313 {
1314 MsqPrintError err( std::cout );
1315 std::vector<size_t> handles, indices1, indices2;
1316 std::vector<Vector3D> grad1, grad2;
1317 double qm_val1, qm_val2;
1318 bool rval;
1319
1320 qm->get_evaluations( pd, handles, false, err );
1321 CPPUNIT_ASSERT( !MSQ_CHKERR(err) );
1322 CPPUNIT_ASSERT( !handles.empty() );
1323 for (size_t j = 0; j < handles.size(); ++j) {
1324 rval = qm->QualityMetric::evaluate_with_gradient( pd, handles[j], qm_val1, indices1, grad1, err );
1325 CPPUNIT_ASSERT( !MSQ_CHKERR(err) );
1326 CPPUNIT_ASSERT( rval );
1327 rval = qm->evaluate_with_gradient( pd, handles[j], qm_val2, indices2, grad2, err );
1328 CPPUNIT_ASSERT( !MSQ_CHKERR(err) );
1329 CPPUNIT_ASSERT( rval );
1330
1331 CPPUNIT_ASSERT_DOUBLES_EQUAL( qm_val1, qm_val2, 1e-6 );
1332 CPPUNIT_ASSERT_EQUAL( indices1.size(), indices2.size() );
1333 CPPUNIT_ASSERT( !indices1.empty() );
1334
1335 std::vector<size_t>::iterator it1, it2;
1336 for (it1 = indices1.begin(); it1 != indices1.end(); ++it1) {
1337 it2 = std::find( indices2.begin(), indices2.end(), *it1 );
1338 CPPUNIT_ASSERT( it2 != indices2.end() );
1339
1340 size_t idx1 = it1 - indices1.begin();
1341 size_t idx2 = it2 - indices2.begin();
1342 CPPUNIT_ASSERT_VECTORS_EQUAL( grad1[idx1], grad2[idx2], 0.01 );
1343 }
1344 }
1345 }
1346
compare_analytical_and_numerical_gradients(QualityMetric * qm)1347 void QualityMetricTester::compare_analytical_and_numerical_gradients( QualityMetric* qm )
1348 {
1349 PatchData pd;
1350 if (mSettings)
1351 pd.attach_settings(mSettings);
1352
1353 CPPUNIT_ASSERT( !types.empty() );
1354 for (size_t i = 0; i < types.size(); ++i) {
1355 get_nonideal_element( types[i], pd, true );
1356 compare_analytical_and_numerical_gradients( qm, pd );
1357 }
1358 }
1359
1360
compare_analytical_and_numerical_hessians(QualityMetric * qm,PatchData & pd)1361 void QualityMetricTester::compare_analytical_and_numerical_hessians(
1362 QualityMetric* qm, PatchData& pd )
1363 {
1364 MsqPrintError err( std::cout );
1365 std::vector<size_t> handles, indices1, indices2;
1366 std::vector<Vector3D> grad1, grad2;
1367 std::vector<Matrix3D> Hess1, Hess2;
1368 double qm_val1, qm_val2;
1369 bool rval;
1370
1371 qm->get_evaluations( pd, handles, false, err );
1372 CPPUNIT_ASSERT( !MSQ_CHKERR(err) );
1373 CPPUNIT_ASSERT( !handles.empty() );
1374 for (size_t j = 0; j < handles.size(); ++j) {
1375 rval = qm->QualityMetric::evaluate_with_Hessian( pd, handles[j], qm_val1, indices1, grad1, Hess1, err );
1376 CPPUNIT_ASSERT( !MSQ_CHKERR(err) );
1377 CPPUNIT_ASSERT( rval );
1378 rval = qm->evaluate_with_Hessian( pd, handles[j], qm_val2, indices2, grad2, Hess2, err );
1379 CPPUNIT_ASSERT( !MSQ_CHKERR(err) );
1380 CPPUNIT_ASSERT( rval );
1381
1382 CPPUNIT_ASSERT_DOUBLES_EQUAL( qm_val1, qm_val2, 1e-6 );
1383 CPPUNIT_ASSERT_EQUAL( indices1.size(), indices2.size() );
1384 CPPUNIT_ASSERT( !indices1.empty() );
1385
1386 std::vector<size_t>::iterator it;
1387 unsigned h = 0;
1388 for (unsigned r = 0; r < indices1.size(); ++r) {
1389 it = std::find( indices2.begin(), indices2.end(), indices1[r] );
1390 CPPUNIT_ASSERT( it != indices2.end() );
1391 unsigned r2 = it - indices2.begin();
1392
1393 for (unsigned c = r; c < indices1.size(); ++c, ++h) {
1394 it = std::find( indices2.begin(), indices2.end(), indices1[c] );
1395 CPPUNIT_ASSERT( it != indices2.end() );
1396 unsigned c2 = it - indices2.begin();
1397
1398 unsigned h2;
1399 if (r2 <= c2)
1400 h2 = indices2.size()*r - r*(r+1)/2 + c;
1401 else
1402 h2 = indices2.size()*c - c*(c+1)/2 + r;
1403
1404 //if (!utest_mat_equal(Hess1[h],Hess2[h2],0.001))
1405 // assert(false);
1406 CPPUNIT_ASSERT_MATRICES_EQUAL( Hess1[h], Hess2[h2], 0.001 );
1407 }
1408 }
1409 }
1410 }
1411
compare_analytical_and_numerical_hessians(QualityMetric * qm)1412 void QualityMetricTester::compare_analytical_and_numerical_hessians( QualityMetric* qm )
1413 {
1414 PatchData pd;
1415 if (mSettings)
1416 pd.attach_settings(mSettings);
1417
1418 CPPUNIT_ASSERT( !types.empty() );
1419 for (size_t i = 0; i < types.size(); ++i) {
1420 get_nonideal_element( types[i], pd, true );
1421 compare_analytical_and_numerical_hessians( qm, pd );
1422 }
1423 }
1424
compare_analytical_and_numerical_diagonals(QualityMetric * qm,PatchData & pd)1425 void QualityMetricTester::compare_analytical_and_numerical_diagonals(
1426 QualityMetric* qm, PatchData& pd )
1427 {
1428 MsqPrintError err( std::cout );
1429 std::vector<size_t> handles, indices1, indices2;
1430 std::vector<Vector3D> grad1, grad2;
1431 std::vector<Matrix3D> Hess1;
1432 std::vector<SymMatrix3D> Hess2;
1433 double qm_val1, qm_val2;
1434 bool rval;
1435
1436 qm->get_evaluations( pd, handles, false, err );
1437 CPPUNIT_ASSERT( !MSQ_CHKERR(err) );
1438 CPPUNIT_ASSERT( !handles.empty() );
1439 for (size_t j = 0; j < handles.size(); ++j) {
1440 rval = qm->QualityMetric::evaluate_with_Hessian( pd, handles[j], qm_val1, indices1, grad1, Hess1, err );
1441 CPPUNIT_ASSERT( !MSQ_CHKERR(err) );
1442 CPPUNIT_ASSERT( rval );
1443 rval = qm->evaluate_with_Hessian_diagonal( pd, handles[j], qm_val2, indices2, grad2, Hess2, err );
1444 CPPUNIT_ASSERT( !MSQ_CHKERR(err) );
1445 CPPUNIT_ASSERT( rval );
1446
1447 CPPUNIT_ASSERT_DOUBLES_EQUAL( qm_val1, qm_val2, 1e-6 );
1448 CPPUNIT_ASSERT_EQUAL( indices1.size(), indices2.size() );
1449 CPPUNIT_ASSERT( !indices1.empty() );
1450 CPPUNIT_ASSERT_EQUAL( indices1.size() * (indices1.size()+1) / 2, Hess1.size() );
1451 CPPUNIT_ASSERT_EQUAL( indices2.size(), Hess2.size() );
1452
1453 size_t h = 0;
1454 std::vector<size_t>::iterator it;
1455 for (unsigned r = 0; r < indices1.size(); ++r) {
1456 it = std::find( indices2.begin(), indices2.end(), indices1[r] );
1457 CPPUNIT_ASSERT( it != indices2.end() );
1458 unsigned r2 = it - indices2.begin();
1459 //if (!utest_mat_equal(Hess1[h],Hess2[r2],0.001))
1460 // assert(false);
1461 CPPUNIT_ASSERT_MATRICES_EQUAL( Hess1[h], Hess2[r2], 0.001 );
1462 h += indices1.size() - r;
1463 }
1464 }
1465 }
1466
compare_analytical_and_numerical_diagonals(QualityMetric * qm)1467 void QualityMetricTester::compare_analytical_and_numerical_diagonals( QualityMetric* qm )
1468 {
1469 PatchData pd;
1470 if (mSettings)
1471 pd.attach_settings(mSettings);
1472 MsqError err;
1473
1474 CPPUNIT_ASSERT( !types.empty() );
1475 for (size_t i = 0; i < types.size(); ++i) {
1476 get_nonideal_element( types[i], pd, true );
1477 compare_analytical_and_numerical_diagonals( qm, pd );
1478 }
1479 }
1480
1481 class PatchTranslate : public QualityMetricTester::PatchXform
1482 {
1483 Vector3D offset;
1484 public:
PatchTranslate(Vector3D s)1485 inline PatchTranslate( Vector3D s ) : offset(s) {}
1486 virtual ~PatchTranslate();
1487 virtual void xform( PatchData& pd, PlanarDomain* dom ) ;
1488 virtual void xform_grad( std::vector<Vector3D>& );
1489 };
1490
1491
~PatchTranslate()1492 PatchTranslate::~PatchTranslate() {}
xform(PatchData & pd,PlanarDomain * dom)1493 void PatchTranslate::xform( PatchData& pd, PlanarDomain* dom )
1494 {
1495 // If patch has associated planar geometry, translate that too
1496 MeshDomain* geom = pd.get_domain();
1497 if (geom) {
1498 PlanarDomain* plane = dynamic_cast<PlanarDomain*>(geom);
1499 CPPUNIT_ASSERT(plane);
1500 CPPUNIT_ASSERT(dom);
1501
1502 Vector3D norm = plane->get_normal();
1503 Vector3D point = plane->get_origin() + offset;
1504 dom->set_plane( norm, point );
1505 pd.set_domain( dom );
1506 }
1507
1508 // Translate mesh vertices
1509 MsqError err;
1510 for (size_t i = 0; i < pd.num_nodes(); ++i)
1511 pd.move_vertex( offset, i, err );
1512 }
xform_grad(std::vector<Vector3D> &)1513 void PatchTranslate::xform_grad( std::vector<Vector3D>& ) {}
1514
1515
1516 class PatchScale : public QualityMetricTester::PatchXform
1517 {
1518 double scaleFactor;
1519 public:
PatchScale(double s)1520 inline PatchScale( double s ) : scaleFactor(s) {}
1521 virtual ~PatchScale();
1522 virtual void xform( PatchData& pd, PlanarDomain* dom ) ;
1523 virtual void xform_grad( std::vector<Vector3D>& );
1524 };
1525
~PatchScale()1526 PatchScale::~PatchScale() {}
xform(PatchData & pd,PlanarDomain * dom)1527 void PatchScale::xform( PatchData& pd, PlanarDomain* dom )
1528 {
1529 // If patch has associated planar geometry, scale distance from origin
1530 MeshDomain* geom = pd.get_domain();
1531 if (geom) {
1532 PlanarDomain* plane = dynamic_cast<PlanarDomain*>(geom);
1533 CPPUNIT_ASSERT(plane);
1534 CPPUNIT_ASSERT(dom);
1535
1536 Vector3D norm = plane->get_normal();
1537 Vector3D point = scaleFactor * plane->get_origin();
1538 dom->set_plane( norm, point );
1539 pd.set_domain( dom );
1540 }
1541
1542 // Scale about origin
1543 MsqError err;
1544 for (size_t i = 0; i < pd.num_nodes(); ++i)
1545 pd.set_vertex_coordinates( pd.vertex_by_index(i) * scaleFactor, i, err );
1546 }
xform_grad(std::vector<Vector3D> &)1547 void PatchScale::xform_grad( std::vector<Vector3D>& ) {}
1548
1549 class PatchRotate : public QualityMetricTester::PatchXform
1550 {
1551 Matrix3D rotation;
1552 public:
PatchRotate(Matrix3D s)1553 inline PatchRotate( Matrix3D s ) : rotation(s) {}
1554 virtual ~PatchRotate();
1555 virtual void xform( PatchData& pd, PlanarDomain* dom ) ;
1556 virtual void xform_grad( std::vector<Vector3D>& grad );
1557 };
1558
~PatchRotate()1559 PatchRotate::~PatchRotate() {}
xform(PatchData & pd,PlanarDomain * dom)1560 void PatchRotate::xform( PatchData& pd, PlanarDomain* dom )
1561 {
1562 // If patch has associated planar geometry, rotate that also
1563 MeshDomain* geom = pd.get_domain();
1564 if (geom) {
1565 PlanarDomain* plane = dynamic_cast<PlanarDomain*>(geom);
1566 CPPUNIT_ASSERT(plane);
1567 CPPUNIT_ASSERT(dom);
1568
1569 Vector3D norm = rotation * plane->get_normal();
1570 Vector3D point = rotation * plane->get_origin();
1571 dom->set_plane( norm, point );
1572 pd.set_domain( dom );
1573 }
1574
1575 // Scale about origin
1576 MsqError err;
1577 for (size_t i = 0; i < pd.num_nodes(); ++i)
1578 pd.set_vertex_coordinates( rotation * pd.vertex_by_index(i), i, err );
1579 }
xform_grad(std::vector<Vector3D> & grad)1580 void PatchRotate::xform_grad( std::vector<Vector3D>& grad )
1581 {
1582 for (unsigned i = 0; i < grad.size(); ++i)
1583 grad[i] = rotation * grad[i];
1584 }
1585
test_transform_invariant(QualityMetric * qm,QualityMetricTester::PatchXform & xform,bool untangle)1586 void QualityMetricTester::test_transform_invariant( QualityMetric* qm,
1587 QualityMetricTester::PatchXform& xform,
1588 bool untangle )
1589 {
1590 MsqPrintError err( std::cout );
1591 PatchData pd;
1592 if (mSettings)
1593 pd.attach_settings(mSettings);
1594 std::vector<size_t> handles;
1595 double qmval;
1596 bool rval;
1597 PlanarDomain dom(Vector3D(0,0,1),Vector3D(0,0,0));
1598 size_t i, j;
1599
1600 CPPUNIT_ASSERT( !types.empty() );
1601 for (i = 0; i < types.size(); ++i) {
1602 if (untangle)
1603 get_inverted_element( types[i], pd );
1604 else
1605 get_nonideal_element( types[i], pd );
1606 qm->get_evaluations( pd, handles, false, err );
1607 CPPUNIT_ASSERT( !MSQ_CHKERR(err) );
1608 CPPUNIT_ASSERT( !handles.empty() );
1609
1610 // get value(s) for ideal element
1611 std::vector<double> values(handles.size());
1612 for (j = 0; j < handles.size(); ++j) {
1613 rval = qm->evaluate( pd, handles[j], values[j], err );
1614 CPPUNIT_ASSERT(!MSQ_CHKERR(err));
1615 CPPUNIT_ASSERT(rval);
1616 }
1617
1618 // compare values for transformed patch
1619 xform.xform( pd, &dom );
1620 for (j = 0; j < handles.size(); ++j) {
1621 rval = qm->evaluate( pd, handles[j], qmval, err );
1622 CPPUNIT_ASSERT(!MSQ_CHKERR(err));
1623 CPPUNIT_ASSERT(rval);
1624 CPPUNIT_ASSERT_DOUBLES_EQUAL( values[j], qmval, 1e-3 );
1625 }
1626 }
1627 }
1628
test_grad_transform_invariant(QualityMetric * qm,QualityMetricTester::PatchXform & xform,bool untangle)1629 void QualityMetricTester::test_grad_transform_invariant( QualityMetric* qm,
1630 QualityMetricTester::PatchXform& xform,
1631 bool untangle )
1632 {
1633 MsqPrintError err( std::cout );
1634 PatchData pd;
1635 if (mSettings)
1636 pd.attach_settings(mSettings);
1637 std::vector<size_t> handles, indices2;
1638 std::vector<Vector3D> grad2;
1639 double qmval;
1640 bool rval;
1641 PlanarDomain dom(Vector3D(0,0,1),Vector3D(0,0,0));
1642 size_t i, j, k;
1643
1644 CPPUNIT_ASSERT( !types.empty() );
1645 for (i = 0; i < types.size(); ++i) {
1646 if (untangle)
1647 get_inverted_element( types[i], pd );
1648 else
1649 get_nonideal_element( types[i], pd );
1650 qm->get_evaluations( pd, handles, false, err );
1651 CPPUNIT_ASSERT( !MSQ_CHKERR(err) );
1652 CPPUNIT_ASSERT( !handles.empty() );
1653
1654 // get value(s) for orignal patch
1655 std::vector< std::vector<Vector3D> > grad1(handles.size());
1656 std::vector< std::vector<size_t> > indices1(handles.size());
1657 for (j = 0; j < handles.size(); ++j) {
1658 rval = qm->evaluate_with_gradient( pd, handles[j], qmval, indices1[j], grad1[j], err );
1659 CPPUNIT_ASSERT(!MSQ_CHKERR(err));
1660 CPPUNIT_ASSERT(rval);
1661 CPPUNIT_ASSERT( !grad1[j].empty() );
1662 }
1663
1664 // transform patch
1665 xform.xform( pd, &dom );
1666
1667 // compare values
1668 for (j = 0; j < handles.size(); ++j) {
1669 rval = qm->evaluate_with_gradient( pd, handles[j], qmval, indices2, grad2, err );
1670 CPPUNIT_ASSERT(!MSQ_CHKERR(err));
1671 CPPUNIT_ASSERT(rval);
1672
1673 CPPUNIT_ASSERT_EQUAL( indices1[j].size(), indices2.size() );
1674 CPPUNIT_ASSERT_EQUAL( grad1[j].size(), grad2.size() );
1675 CPPUNIT_ASSERT( indices1[j] == indices2 );
1676
1677 xform.xform_grad( grad1[j] );
1678 for (k = 0; k < grad2.size(); ++k)
1679 CPPUNIT_ASSERT_VECTORS_EQUAL( grad1[j][k], grad2[k], 1e-3 );
1680 }
1681 }
1682 }
1683
test_hessian_transform_invariant(QualityMetric * qm,QualityMetricTester::PatchXform & xform,bool untangle)1684 void QualityMetricTester::test_hessian_transform_invariant( QualityMetric* qm,
1685 QualityMetricTester::PatchXform& xform,
1686 bool untangle )
1687 {
1688 MsqPrintError err( std::cout );
1689 PatchData pd;
1690 if (mSettings)
1691 pd.attach_settings(mSettings);
1692 std::vector<size_t> handles, indices2;
1693 std::vector<Vector3D> grad1, grad2;
1694 std::vector<Matrix3D> hess2;
1695 double qmval;
1696 bool rval;
1697 PlanarDomain dom(Vector3D(0,0,1),Vector3D(0,0,0));
1698 size_t i, j, k;
1699
1700 CPPUNIT_ASSERT( !types.empty() );
1701 for (i = 0; i < types.size(); ++i) {
1702 if (untangle)
1703 get_inverted_element( types[i], pd );
1704 else
1705 get_nonideal_element( types[i], pd );
1706 qm->get_evaluations( pd, handles, false, err );
1707 CPPUNIT_ASSERT( !MSQ_CHKERR(err) );
1708 CPPUNIT_ASSERT( !handles.empty() );
1709
1710 std::vector< std::vector<size_t> > indices1(handles.size());
1711 std::vector< std::vector<Matrix3D> > hess1(handles.size());
1712 for (j = 0; j < handles.size(); ++j) {
1713 rval = qm->evaluate_with_Hessian( pd, handles[j], qmval, indices1[j], grad1, hess1[j], err );
1714 CPPUNIT_ASSERT(!MSQ_CHKERR(err));
1715 CPPUNIT_ASSERT(rval);
1716 CPPUNIT_ASSERT(!hess1[j].empty());
1717 }
1718
1719 xform.xform( pd, &dom );
1720
1721 for (j = 0; j < handles.size(); ++j) {
1722 rval = qm->evaluate_with_Hessian( pd, handles[j], qmval, indices2, grad2, hess2, err );
1723 CPPUNIT_ASSERT(!MSQ_CHKERR(err));
1724 CPPUNIT_ASSERT(rval);
1725
1726 CPPUNIT_ASSERT_EQUAL( indices1[j].size(), indices2.size() );
1727 CPPUNIT_ASSERT_EQUAL( hess1[j].size(), hess2.size() );
1728 CPPUNIT_ASSERT( indices1[j] == indices2 );
1729
1730 for (k = 0; k < hess2.size(); ++k)
1731 CPPUNIT_ASSERT_MATRICES_EQUAL( hess1[j][k], hess2[k], 1e-3 );
1732 }
1733 }
1734 }
1735
test_location_invariant(QualityMetric * qm,bool untangle)1736 void QualityMetricTester::test_location_invariant( QualityMetric* qm, bool untangle )
1737 {
1738 PatchTranslate xform( Vector3D(1,1,-1) );
1739 test_transform_invariant( qm, xform, untangle );
1740 }
1741
test_scale_invariant(QualityMetric * qm,bool untangle)1742 void QualityMetricTester::test_scale_invariant( QualityMetric* qm, bool untangle )
1743 {
1744 PatchScale scale1( 0.5 ), scale2( 2.0 );
1745 test_transform_invariant( qm, scale1, untangle );
1746 test_transform_invariant( qm, scale2, untangle );
1747 }
1748
test_orient_invariant(QualityMetric * qm,bool untangle)1749 void QualityMetricTester::test_orient_invariant( QualityMetric* qm, bool untangle )
1750 {
1751 // rotate 45 degrees about x-axis
1752 const double f = std::sqrt(2.0)/2;
1753 Matrix3D rotation( 1, 0, 0,
1754 0, f, f,
1755 0,-f, f );
1756 PatchRotate xform(rotation);
1757 test_transform_invariant( qm, xform, untangle );
1758 }
1759
test_grad_location_invariant(QualityMetric * qm,bool untangle)1760 void QualityMetricTester::test_grad_location_invariant( QualityMetric* qm, bool untangle )
1761 {
1762 PatchTranslate xform( Vector3D(1,1,-1) );
1763 test_grad_transform_invariant( qm, xform, untangle );
1764 }
1765
test_hessian_location_invariant(QualityMetric * qm,bool untangle)1766 void QualityMetricTester::test_hessian_location_invariant( QualityMetric* qm, bool untangle )
1767 {
1768 PatchTranslate xform( Vector3D(1,1,-1) );
1769 test_hessian_transform_invariant( qm, xform, untangle );
1770 }
1771
test_grad_orient_invariant(QualityMetric * qm,bool untangle)1772 void QualityMetricTester::test_grad_orient_invariant( QualityMetric* qm, bool untangle )
1773 {
1774 // rotate 45 degrees about x-axis
1775 const double f = std::sqrt(2.0)/2;
1776 Matrix3D rotation( 1, 0, 0,
1777 0, f, f,
1778 0,-f, f );
1779 PatchRotate xform(rotation);
1780 test_grad_transform_invariant( qm, xform, untangle );
1781 }
1782
test_measures_transform(QualityMetric * qm,QualityMetricTester::PatchXform & xform,bool unit_area)1783 void QualityMetricTester::test_measures_transform( QualityMetric* qm,
1784 QualityMetricTester::PatchXform& xform,
1785 bool unit_area )
1786 {
1787 MsqPrintError err( std::cout );
1788 PatchData pd;
1789 if (mSettings)
1790 pd.attach_settings(mSettings);
1791 std::vector<size_t> handles;
1792 double qmval1, qmval2;
1793 bool rval;
1794 PlanarDomain dom(Vector3D(0,0,1),Vector3D(0,0,0));
1795 bool decreasing = false;
1796 if (qm->get_negate_flag() != 1) {
1797 CPPUNIT_ASSERT_EQUAL( qm->get_negate_flag(), -1 );
1798 decreasing = true;
1799 }
1800
1801 CPPUNIT_ASSERT( !types.empty() );
1802 for (size_t i = 0; i < types.size(); ++i) {
1803 get_ideal_element( types[i], unit_area, pd );
1804 qm->get_evaluations( pd, handles, false, err );
1805 CPPUNIT_ASSERT( !MSQ_CHKERR(err) );
1806 CPPUNIT_ASSERT( !handles.empty() );
1807
1808 rval = qm->evaluate( pd, handles[0], qmval1, err );
1809 CPPUNIT_ASSERT(!MSQ_CHKERR(err));
1810 CPPUNIT_ASSERT(rval);
1811
1812 xform.xform( pd, &dom );
1813 rval = qm->evaluate( pd, handles[0], qmval2, err );
1814 CPPUNIT_ASSERT(!MSQ_CHKERR(err));
1815 CPPUNIT_ASSERT(rval);
1816 if (decreasing)
1817 CPPUNIT_ASSERT( qmval2 < qmval1 );
1818 else
1819 CPPUNIT_ASSERT( qmval2 > qmval1 );
1820 }
1821 }
1822
test_measures_size(QualityMetric * qm,bool unit_area)1823 void QualityMetricTester::test_measures_size( QualityMetric* qm, bool unit_area )
1824 {
1825 PatchScale s1(0.5), s2(2.0);
1826 test_measures_transform( qm, s1, unit_area );
1827 test_measures_transform( qm, s2, unit_area );
1828 }
1829
test_measures_in_plane_orientation(QualityMetric * qm)1830 void QualityMetricTester::test_measures_in_plane_orientation( QualityMetric* qm )
1831 {
1832 // rotate 45 degrees about z-axis
1833 const double f = std::sqrt(2.0)/2;
1834 Matrix3D rotation( f,-f, 0,
1835 f, f, 0,
1836 0, 0, 1 );
1837 PatchRotate xform( rotation );
1838 test_measures_transform( qm, xform, false );
1839 }
1840
test_measures_out_of_plane_orientation(QualityMetric * qm)1841 void QualityMetricTester::test_measures_out_of_plane_orientation( QualityMetric* qm )
1842 {
1843 // rotate 45 degrees about z-axis
1844 const double f = std::sqrt(2.0)/2;
1845 Matrix3D rotation( 1, 0, 0,
1846 0, f, f,
1847 0,-f, f );
1848 PatchRotate xform( rotation );
1849 test_measures_transform( qm, xform, false );
1850 }
1851
1852 typedef void (QualityMetricTester::*patch_func_t)( EntityTopology, PatchData& );
test_bad_element(QualityMetricTester * instance,QualityMetric * qm,patch_func_t type_func,bool should_succeed,const std::vector<EntityTopology> & types,const Settings * mSettings)1853 static void test_bad_element( QualityMetricTester* instance,
1854 QualityMetric* qm,
1855 patch_func_t type_func,
1856 bool should_succeed,
1857 const std::vector<EntityTopology>& types,
1858 const Settings* mSettings )
1859 {
1860 MsqPrintError err( std::cout );
1861 PatchData pd;
1862 if (mSettings)
1863 pd.attach_settings(mSettings);
1864 std::vector<size_t> handles;
1865 double qmval;
1866 bool rval;
1867
1868 CPPUNIT_ASSERT( !types.empty() );
1869 for (size_t i = 0; i < types.size(); ++i) {
1870 (instance->*type_func)( types[i], pd );
1871 qm->get_evaluations( pd, handles, false, err );
1872 CPPUNIT_ASSERT( !MSQ_CHKERR(err) );
1873 CPPUNIT_ASSERT( !handles.empty() );
1874
1875 rval = true;
1876 for (size_t j = 0; j < handles.size(); ++j) {
1877 rval = rval && qm->evaluate( pd, handles[j], qmval, err );
1878 CPPUNIT_ASSERT(!MSQ_CHKERR(err));
1879 }
1880 CPPUNIT_ASSERT_EQUAL( should_succeed, rval );
1881 }
1882 }
1883
test_evaluate_inverted_element(QualityMetric * qm,bool should_succeed)1884 void QualityMetricTester::test_evaluate_inverted_element( QualityMetric* qm, bool should_succeed )
1885 { test_bad_element( this, qm, &QualityMetricTester::get_inverted_element, should_succeed, types, mSettings ); }
test_evaluate_degenerate_element(QualityMetric * qm,bool should_succeed)1886 void QualityMetricTester::test_evaluate_degenerate_element( QualityMetric* qm, bool should_succeed )
1887 { test_bad_element( this, qm, &QualityMetricTester::get_degenerate_element, should_succeed, types, mSettings ); }
test_evaluate_zero_element(QualityMetric * qm,bool should_succeed)1888 void QualityMetricTester::test_evaluate_zero_element( QualityMetric* qm, bool should_succeed )
1889 { test_bad_element( this, qm, &QualityMetricTester::get_zero_element, should_succeed, types, mSettings ); }
1890
test_measures_quality(QualityMetric * qm)1891 void QualityMetricTester::test_measures_quality( QualityMetric* qm )
1892 {
1893 MsqPrintError err( std::cout );
1894 PatchData pd;
1895 if (mSettings)
1896 pd.attach_settings(mSettings);
1897 std::vector<size_t> handles;
1898 double qmval1, qmval2;
1899 bool rval;
1900 bool decreasing = false;
1901 if (qm->get_negate_flag() != 1) {
1902 CPPUNIT_ASSERT_EQUAL( qm->get_negate_flag(), -1 );
1903 decreasing = true;
1904 }
1905
1906 CPPUNIT_ASSERT( !types.empty() );
1907 for (size_t i = 0; i < types.size(); ++i) {
1908 // get a patch
1909 get_ideal_element( types[i], false, pd );
1910 qm->get_evaluations( pd, handles, false, err );
1911 CPPUNIT_ASSERT( !MSQ_CHKERR(err) );
1912 CPPUNIT_ASSERT( !handles.empty() );
1913 // evaluate for ideal element
1914 rval = qm->evaluate( pd, handles[0], qmval1, err );
1915 CPPUNIT_ASSERT(!MSQ_CHKERR(err));
1916 CPPUNIT_ASSERT( rval );
1917 // calculate element centroid
1918 const size_t* verts = pd.element_by_index(0).get_vertex_index_array();
1919 const MsqVertex* coords = pd.get_vertex_array(err);
1920 size_t n = pd.element_by_index(0).vertex_count();
1921 Vector3D cent(0,0,0);
1922 Vector3D coords0 = coords[verts[0]];
1923 for (unsigned j = 0; j < n; ++j)
1924 cent += coords[verts[j]];
1925 cent /= n;
1926 // evaluate for non-ideal element (decreased area)
1927 pd.set_vertex_coordinates( 0.5*(cent + coords0), verts[0], err );
1928 rval = qm->evaluate( pd, handles[0], qmval2, err );
1929 CPPUNIT_ASSERT(!MSQ_CHKERR(err));
1930 CPPUNIT_ASSERT( rval );
1931 // check values
1932 if (decreasing) {
1933 CPPUNIT_ASSERT( qmval2 < qmval1 );
1934 }
1935 else {
1936 CPPUNIT_ASSERT( qmval2 > qmval1 );
1937 }
1938 // evaluate for non-ideal element (increased area)
1939 pd.set_vertex_coordinates( 3*coords0 - 2*cent, verts[0], err );
1940 rval = qm->evaluate( pd, handles[0], qmval2, err );
1941 CPPUNIT_ASSERT(!MSQ_CHKERR(err));
1942 CPPUNIT_ASSERT( rval );
1943 // check values
1944 if (decreasing) {
1945 CPPUNIT_ASSERT( qmval2 < qmval1 );
1946 }
1947 else {
1948 CPPUNIT_ASSERT( qmval2 > qmval1 );
1949 }
1950 }
1951 }
1952
test_domain_deviation_quality(QualityMetric * qm)1953 void QualityMetricTester::test_domain_deviation_quality( QualityMetric* qm )
1954 {
1955 MsqPrintError err( std::cout );
1956 PatchData pd;
1957 if (mSettings)
1958 pd.attach_settings(mSettings);
1959 std::vector<size_t> handles, indices;
1960 size_t i;
1961 double qmval1, qmval2;
1962 bool rval;
1963 bool decreasing = false;
1964 if (qm->get_negate_flag() != 1) {
1965 CPPUNIT_ASSERT_EQUAL( qm->get_negate_flag(), -1 );
1966 decreasing = true;
1967 }
1968
1969 CPPUNIT_ASSERT( type_is_supported(TRIANGLE)||type_is_supported(QUADRILATERAL) );
1970
1971 if (type_is_supported(TRIANGLE)) {
1972 // get a patch
1973 get_ideal_tris( pd, false );
1974 qm->get_evaluations( pd, handles, false, err );
1975 CPPUNIT_ASSERT( !MSQ_CHKERR(err) );
1976 CPPUNIT_ASSERT( !handles.empty() );
1977 // find an evaluation that depends on the free vertex
1978 for (i = 0; i < handles.size(); ++i) {
1979 rval = qm->evaluate_with_indices( pd, handles[i], qmval1, indices, err );
1980 CPPUNIT_ASSERT(!MSQ_CHKERR(err));
1981 CPPUNIT_ASSERT(rval);
1982 if (!indices.empty())
1983 break;
1984 }
1985 CPPUNIT_ASSERT_EQUAL( (size_t)1, indices.size() ); // one free vertex in patch
1986 CPPUNIT_ASSERT_EQUAL( (size_t)0, indices.front() ); // one free vertex in patch
1987 const size_t handle = handles[i];
1988 // rotate patch about origin such that elements are no longer in the plane
1989 const Matrix3D rot( M_PI/4.0, Vector3D( 1, 1, 0 ) );
1990 for (i = 0; i < pd.num_nodes(); ++i)
1991 pd.set_vertex_coordinates( rot * pd.vertex_by_index(i), i, err );
1992 // evaluate for rotated patch
1993 rval = qm->evaluate( pd, handle, qmval2, err );
1994 CPPUNIT_ASSERT(!MSQ_CHKERR(err));
1995 CPPUNIT_ASSERT(rval);
1996 // check quality
1997 if (decreasing) {
1998 CPPUNIT_ASSERT( qmval1 > qmval2 );
1999 }
2000 else {
2001 CPPUNIT_ASSERT( qmval1 < qmval2 );
2002 }
2003 }
2004
2005 if (type_is_supported(QUADRILATERAL)) {
2006 // get a patch
2007 get_ideal_quads( pd );
2008 qm->get_evaluations( pd, handles, false, err );
2009 CPPUNIT_ASSERT( !MSQ_CHKERR(err) );
2010 CPPUNIT_ASSERT( !handles.empty() );
2011 // find an evaluation that depends on the free vertex
2012 for (i = 0; i < handles.size(); ++i) {
2013 rval = qm->evaluate_with_indices( pd, handles[i], qmval1, indices, err );
2014 CPPUNIT_ASSERT(!MSQ_CHKERR(err));
2015 CPPUNIT_ASSERT(rval);
2016 if (!indices.empty())
2017 break;
2018 }
2019 CPPUNIT_ASSERT_EQUAL( (size_t)1, indices.size() ); // one free vertex in patch
2020 CPPUNIT_ASSERT_EQUAL( (size_t)0, indices.front() ); // one free vertex in patch
2021 const size_t handle = handles[i];
2022 // rotate patch about origin such that elements are no longer in the plane
2023 const Matrix3D rot( M_PI/3.0, Vector3D( 1, 0, 0 ) );
2024 for (i = 0; i < pd.num_nodes(); ++i)
2025 pd.set_vertex_coordinates( rot * pd.vertex_by_index(i), i, err );
2026 // evaluate for rotated patch
2027 rval = qm->evaluate( pd, handle, qmval2, err );
2028 CPPUNIT_ASSERT(!MSQ_CHKERR(err));
2029 CPPUNIT_ASSERT(rval);
2030 // check quality
2031 if (decreasing) {
2032 CPPUNIT_ASSERT( qmval1 > qmval2 );
2033 }
2034 else {
2035 CPPUNIT_ASSERT( qmval1 < qmval2 );
2036 }
2037 }
2038 }
2039
2040
test_domain_deviation_gradient(QualityMetric * qm)2041 void QualityMetricTester::test_domain_deviation_gradient( QualityMetric* qm )
2042 {
2043 MsqPrintError err( std::cout );
2044 PatchData pd;
2045 if (mSettings)
2046 pd.attach_settings(mSettings);
2047 std::vector<size_t> handles, indices;
2048 std::vector<Vector3D> grad;
2049 size_t i;
2050 double qmval;
2051 bool rval;
2052
2053 CPPUNIT_ASSERT( type_is_supported(TRIANGLE)||type_is_supported(QUADRILATERAL) );
2054
2055 if (type_is_supported(TRIANGLE)) {
2056 // get a patch
2057 get_ideal_tris( pd, false );
2058 qm->get_evaluations( pd, handles, false, err );
2059 CPPUNIT_ASSERT( !MSQ_CHKERR(err) );
2060 CPPUNIT_ASSERT( !handles.empty() );
2061 // move the free vertex out of the plane
2062 pd.move_vertex( Vector3D( 0, 0, 0.2 ), 0, err );
2063 // find an evaluation that depends on the free vertex
2064 for (i = 0; i < handles.size(); ++i) {
2065 rval = qm->evaluate_with_gradient( pd, handles[i], qmval, indices, grad, err );
2066 CPPUNIT_ASSERT(!MSQ_CHKERR(err));
2067 CPPUNIT_ASSERT(rval);
2068 if (!grad.empty())
2069 break;
2070 }
2071 CPPUNIT_ASSERT_EQUAL( (size_t)1, indices.size() ); // one free vertex in patch
2072 CPPUNIT_ASSERT_EQUAL( (size_t)0, indices.front() ); // one free vertex in patch
2073 CPPUNIT_ASSERT_EQUAL( (size_t)1, grad.size() );
2074
2075 Vector3D v = qm->get_negate_flag() * Vector3D( 0, 0, 1 );
2076 CPPUNIT_ASSERT( v % grad[0] > DBL_EPSILON );
2077 }
2078
2079 if (type_is_supported(QUADRILATERAL)) {
2080 // get a patch
2081 get_ideal_quads( pd );
2082 qm->get_evaluations( pd, handles, false, err );
2083 CPPUNIT_ASSERT( !MSQ_CHKERR(err) );
2084 CPPUNIT_ASSERT( !handles.empty() );
2085 // move the free vertex out of the plane
2086 pd.move_vertex( Vector3D( 0, 0, 0.2 ), 0, err );
2087 // find an evaluation that depends on the free vertex
2088 for (i = 0; i < handles.size(); ++i) {
2089 rval = qm->evaluate_with_gradient( pd, handles[i], qmval, indices, grad, err );
2090 CPPUNIT_ASSERT(!MSQ_CHKERR(err));
2091 CPPUNIT_ASSERT(rval);
2092 if (!grad.empty())
2093 break;
2094 }
2095 CPPUNIT_ASSERT_EQUAL( (size_t)1, indices.size() ); // one free vertex in patch
2096 CPPUNIT_ASSERT_EQUAL( (size_t)0, indices.front() ); // one free vertex in patch
2097 CPPUNIT_ASSERT_EQUAL( (size_t)1, grad.size() );
2098
2099 Vector3D v = qm->get_negate_flag() * Vector3D( 0, 0, 1 );
2100 CPPUNIT_ASSERT( v % grad[0] > DBL_EPSILON );
2101 }
2102
2103 }
2104
test_gradient_reflects_quality(QualityMetric * qm)2105 void QualityMetricTester::test_gradient_reflects_quality( QualityMetric* qm )
2106 {
2107 MsqPrintError err( std::cout );
2108 PatchData pd;
2109 if (mSettings)
2110 pd.attach_settings(mSettings);
2111 std::vector<size_t> handles, indices;
2112 std::vector<Vector3D> grad;
2113 double qmval;
2114 bool rval;
2115 bool decreasing = false;
2116 if (qm->get_negate_flag() != 1) {
2117 CPPUNIT_ASSERT_EQUAL( qm->get_negate_flag(), -1 );
2118 decreasing = true;
2119 }
2120
2121 CPPUNIT_ASSERT( !types.empty() );
2122 for (size_t i = 0; i < types.size(); ++i) {
2123 // get a patch
2124 get_ideal_element( types[i], false, pd );
2125 qm->get_evaluations( pd, handles, false, err );
2126 CPPUNIT_ASSERT( !MSQ_CHKERR(err) );
2127 CPPUNIT_ASSERT( !handles.empty() );
2128 // make element non-ideal
2129 const size_t* verts = pd.element_by_index(0).get_vertex_index_array();
2130 const MsqVertex* coords = pd.get_vertex_array(err);
2131 size_t n = pd.element_by_index(0).vertex_count();
2132 Vector3D sum(0,0,0);
2133 for (unsigned j = 0; j < n; ++j)
2134 sum += coords[verts[j]];
2135 const Vector3D init_pos = coords[verts[0]];
2136 const Vector3D new_pos = 0.5 * coords[verts[0]] + sum/(2*n);
2137 pd.set_vertex_coordinates( new_pos, verts[0], err );
2138 // evaluate for non-ideal element
2139 rval = qm->evaluate_with_gradient( pd, handles[0], qmval, indices, grad, err );
2140 CPPUNIT_ASSERT(!MSQ_CHKERR(err));
2141 CPPUNIT_ASSERT( rval );
2142 // find gradient for vertex 0
2143 std::vector<size_t>::iterator it = std::find( indices.begin(), indices.end(), verts[0] );
2144 CPPUNIT_ASSERT( it != indices.end() );
2145 Vector3D g = grad[it - indices.begin()];
2146 // check gradient direction
2147 Vector3D v = init_pos - coords[verts[0]];
2148 double dotprod = v % g;
2149 if (decreasing)
2150 CPPUNIT_ASSERT( dotprod > 1e-6 );
2151 else
2152 CPPUNIT_ASSERT( dotprod < -1e-6 );
2153 }
2154 }
2155
test_measures_vertex_quality(QualityMetric * qm)2156 void QualityMetricTester::test_measures_vertex_quality( QualityMetric* qm )
2157 {
2158 MsqPrintError err( std::cout );
2159 PatchData pd;
2160 if (mSettings)
2161 pd.attach_settings(mSettings);
2162 double qmval1, qmval2;
2163 bool rval;
2164 bool increasing = true;
2165 if (qm->get_negate_flag() != 1) {
2166 CPPUNIT_ASSERT_EQUAL( qm->get_negate_flag(), -1 );
2167 increasing = true;
2168 }
2169
2170 unsigned count = 0;
2171 for (size_t i = 0; i < 3; ++i) {
2172 // get a patch
2173 switch (i) {
2174 case 0:
2175 if (!type_is_supported(TRIANGLE)) continue;
2176 get_ideal_tris( pd, false ); break;
2177 case 1:
2178 if (!type_is_supported(QUADRILATERAL)) continue;
2179 get_ideal_quads( pd ); break;
2180 case 2:
2181 if (!type_is_supported(HEXAHEDRON)) continue;
2182 get_ideal_hexes( pd ); break;
2183 }
2184 ++count;
2185
2186 // evaluate for ideal element
2187 rval = qm->evaluate( pd, 0, qmval1, err );
2188 CPPUNIT_ASSERT(!MSQ_CHKERR(err));
2189 CPPUNIT_ASSERT( rval );
2190 // move center vertex
2191 pd.move_vertex( Vector3D(0.3,0.3,0.3), 0, err );
2192 // evaluate for non-ideal element
2193 rval = qm->evaluate( pd, 0, qmval2, err );
2194 CPPUNIT_ASSERT(!MSQ_CHKERR(err));
2195 CPPUNIT_ASSERT( rval );
2196 // check values
2197 if (increasing) {
2198 CPPUNIT_ASSERT( qmval2 > qmval1 );
2199 }
2200 else {
2201 CPPUNIT_ASSERT( qmval2 < qmval1 );
2202 }
2203 }
2204 // make sure we tested something
2205 CPPUNIT_ASSERT( count > 0 );
2206 }
2207
test_vertex_gradient_reflects_quality(QualityMetric * qm)2208 void QualityMetricTester::test_vertex_gradient_reflects_quality( QualityMetric* qm )
2209 {
2210 MsqPrintError err( std::cout );
2211 PatchData pd;
2212 if (mSettings)
2213 pd.attach_settings(mSettings);
2214 std::vector<size_t> indices;
2215 std::vector<Vector3D> grad;
2216 double qmval;
2217 bool rval;
2218 bool increasing = true;
2219 if (qm->get_negate_flag() != 1) {
2220 CPPUNIT_ASSERT_EQUAL( qm->get_negate_flag(), -1 );
2221 increasing = false;
2222 }
2223
2224 unsigned count = 0;
2225 for (size_t i = 0; i < 3; ++i) {
2226 // get a patch
2227 switch (i) {
2228 case 0:
2229 if (!type_is_supported(TRIANGLE)) continue;
2230 get_ideal_tris( pd, false ); break;
2231 case 1:
2232 if (!type_is_supported(QUADRILATERAL)) continue;
2233 get_ideal_quads( pd ); break;
2234 case 2:
2235 if (!type_is_supported(HEXAHEDRON)) continue;
2236 get_ideal_hexes( pd ); break;
2237 }
2238
2239 // move center vertex
2240 pd.move_vertex( Vector3D(0.3,0.3,0.3), 0, err );
2241
2242 // evaluate for ideal non-element
2243 rval = qm->evaluate_with_gradient( pd, 0, qmval, indices, grad, err );
2244 CPPUNIT_ASSERT(!MSQ_CHKERR(err));
2245 CPPUNIT_ASSERT( rval );
2246 // find gradient for vertex 0
2247 if (indices.empty())
2248 continue;
2249 CPPUNIT_ASSERT_EQUAL( (size_t)1, indices.size() ); // only 1 free vertex
2250 CPPUNIT_ASSERT_EQUAL( (size_t)0, indices[0] ); // and that vertex is vertex 0
2251 ++count;
2252 Vector3D g = grad[0];
2253 // check gradient direction
2254 if (increasing)
2255 CPPUNIT_ASSERT( g[2] > 1e-6 );
2256 else
2257 CPPUNIT_ASSERT( g[2] < -1e-6 );
2258 }
2259 // make sure we tested something
2260 CPPUNIT_ASSERT( count > 0 );
2261 }
2262
test_ideal_element_zero_gradient(QualityMetric * qm,bool unit_area)2263 void QualityMetricTester::test_ideal_element_zero_gradient( QualityMetric* qm, bool unit_area )
2264 {
2265 MsqPrintError err( std::cout );
2266 PatchData pd;
2267 if (mSettings)
2268 pd.attach_settings(mSettings);
2269 std::vector<size_t> handles, indices;
2270 std::vector<Vector3D> grad;
2271 double qmval;
2272 bool rval;
2273
2274 CPPUNIT_ASSERT( !types.empty() );
2275 for (size_t i = 0; i < types.size(); ++i) {
2276 // get a patch
2277 get_ideal_element( types[i], unit_area, pd );
2278 qm->get_evaluations( pd, handles, false, err );
2279 CPPUNIT_ASSERT( !MSQ_CHKERR(err) );
2280 CPPUNIT_ASSERT( !handles.empty() );
2281 for (size_t j = 0; j < handles.size(); ++j) {
2282 // evaluate
2283 rval = qm->evaluate_with_gradient( pd, handles[j], qmval, indices, grad, err );
2284 CPPUNIT_ASSERT(!MSQ_CHKERR(err));
2285 CPPUNIT_ASSERT( rval );
2286 CPPUNIT_ASSERT_EQUAL( indices.size(), grad.size() );
2287 // check that all gradients are zero
2288 for (size_t k = 0; k < grad.size(); ++k) {
2289 CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D(0,0,0), grad[k], 1e-4 );
2290 } // for(k < grad.size())
2291 } // for(j < handles.size())
2292 } // for(i < types.size())
2293 }
2294
test_ideal_element_zero_vertex_gradient(QualityMetric * qm,bool unit_area)2295 void QualityMetricTester::test_ideal_element_zero_vertex_gradient(
2296 QualityMetric* qm, bool unit_area )
2297 {
2298 MsqPrintError err( std::cout );
2299 PatchData pd;
2300 if (mSettings)
2301 pd.attach_settings(mSettings);
2302 std::vector<size_t> handles, indices;
2303 std::vector<Vector3D> grad;
2304 double qmval;
2305 bool rval;
2306
2307 unsigned count = 0;
2308 for (size_t i = 0; i < 3; ++i) {
2309 // get a patch
2310 switch (i) {
2311 case 0:
2312 if (!type_is_supported(TRIANGLE)) continue;
2313 get_ideal_tris( pd, unit_area ); break;
2314 case 1:
2315 if (!type_is_supported(QUADRILATERAL)) continue;
2316 get_ideal_quads( pd ); break;
2317 case 2:
2318 if (!type_is_supported(HEXAHEDRON)) continue;
2319 get_ideal_hexes( pd ); break;
2320 }
2321
2322 qm->get_evaluations( pd, handles, false, err );
2323 CPPUNIT_ASSERT( !MSQ_CHKERR(err) );
2324 CPPUNIT_ASSERT( !handles.empty() );
2325 for (size_t j = 0; j < handles.size(); ++j) {
2326 // evaluate
2327 rval = qm->evaluate_with_gradient( pd, handles[j], qmval, indices, grad, err );
2328 CPPUNIT_ASSERT(!MSQ_CHKERR(err));
2329 CPPUNIT_ASSERT( rval );
2330 CPPUNIT_ASSERT_EQUAL( indices.size(), grad.size() );
2331 if (grad.empty())
2332 continue;
2333 ++count;
2334 CPPUNIT_ASSERT_EQUAL( (size_t)1, grad.size() ); // only 1 free vertex in patch
2335 // check that all gradients are zero
2336 CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D(0,0,0), grad[0], 1e-4 );
2337 } // for(j < handles.size())
2338 } // for(i < 3)
2339 // make sure we tested something
2340 CPPUNIT_ASSERT( count > 0 );
2341 }
2342
allocate_matrix(unsigned n)2343 static double** allocate_matrix( unsigned n )
2344 {
2345 unsigned num_ptr = n;
2346 if (num_ptr%2) ++num_ptr;
2347
2348 void* storage = malloc( n*n*sizeof(double) + num_ptr*sizeof(double*) );
2349 double** ptrs = (double**)storage;
2350 double* data = (double*)(ptrs+num_ptr);
2351 for (unsigned i = 0; i < n; ++i)
2352 ptrs[i] = data + i*n;
2353 return ptrs;
2354 }
2355
value(const Matrix3D * blocks,unsigned n,unsigned i,unsigned j)2356 static double value( const Matrix3D* blocks, unsigned n, unsigned i, unsigned j )
2357 {
2358 if (i > j)
2359 std::swap(i,j);
2360
2361 unsigned mi = i / 3;
2362 unsigned mj = j / 3;
2363 unsigned idx = n*mi - mi*(mi+1)/2 + mj;
2364 return blocks[idx][i%3][j%3];
2365 }
2366 /*
2367 static bool positive_definite( const Matrix3D* blocks, unsigned n )
2368 {
2369 // Do Cholesky-Banachiewicz decompositon of the matrix.
2370 // If this results in any imaginary values for diagonal terms,
2371 // then the matrix is not positive-definite.
2372 const int N = 3*n;
2373 double** mat = allocate_matrix(N);
2374 int i, j, k;
2375 double s;
2376 for (i = 0; i < N; ++i)
2377 {
2378 for (j = 0; j < i; ++j)
2379 {
2380 s = value(blocks,n,i,j);
2381 for (k = 0; k < j - 1; ++k)
2382 s -= mat[i][k]*mat[j][k];
2383 mat[i][j] = s / mat[j][j];
2384 }
2385 s = value(blocks,n,i,i);
2386 for (k = 0; k < i - 1; ++k)
2387 s -= mat[i][k]*mat[i][k];
2388 if (s < 0.0)
2389 {
2390 free( mat );
2391 return false;
2392 }
2393 mat[i][i] = sqrt(s);
2394 }
2395 return true;
2396 }
2397 */
2398
2399 /** Test that Hessian is positive-definite for ideal elements */
test_ideal_element_positive_definite_Hessian(QualityMetric * qm,bool unit_area)2400 void QualityMetricTester::test_ideal_element_positive_definite_Hessian(
2401 QualityMetric* qm, bool unit_area )
2402 {
2403 MsqPrintError err( std::cout );
2404 PatchData pd;
2405 if (mSettings)
2406 pd.attach_settings(mSettings);
2407 std::vector<size_t> handles, indices;
2408 std::vector<Vector3D> grad;
2409 std::vector<Matrix3D> Hess;
2410 double qmval;
2411 bool rval;
2412
2413 CPPUNIT_ASSERT( !types.empty() );
2414 for (size_t i = 0; i < types.size(); ++i) {
2415 // get a patch
2416 get_ideal_element( types[i], unit_area, pd );
2417 qm->get_evaluations( pd, handles, false, err );
2418 CPPUNIT_ASSERT( !MSQ_CHKERR(err) );
2419 CPPUNIT_ASSERT( !handles.empty() );
2420 for (size_t j = 0; j < handles.size(); ++j) {
2421 // evaluate
2422 rval = qm->evaluate_with_Hessian( pd, handles[j], qmval, indices, grad, Hess, err );
2423 CPPUNIT_ASSERT(!MSQ_CHKERR(err));
2424 CPPUNIT_ASSERT( rval );
2425 const size_t n = indices.size();
2426 CPPUNIT_ASSERT_EQUAL( n, grad.size() );
2427 CPPUNIT_ASSERT_EQUAL( n*(n+1)/2, Hess.size() );
2428 // if necessary, adjust Hessian for negate_flag();
2429 if (qm->get_negate_flag() == -1)
2430 for (size_t k = 0; k < Hess.size(); ++k)
2431 Hess[k] *= -1.0;
2432 // check that all diagonal blocks are positive definite
2433 size_t h_idx = 0;
2434 for (size_t k = 0; k < n; h_idx+=(n-k), ++k) {
2435 CPPUNIT_ASSERT( Hess[h_idx].positive_definite() );
2436 }
2437 // check that the entire Hessian is positive definite
2438 //CPPUNIT_ASSERT( positive_definite( arrptr(Hess), n ) );
2439 } // for(j < handles.size())
2440 } // for(i < types.size())
2441 }
2442
2443 /** Test that diagonal bocks of Hessian are symmetric */
test_symmetric_Hessian_diagonal_blocks(QualityMetric * qm)2444 void QualityMetricTester::test_symmetric_Hessian_diagonal_blocks( QualityMetric* qm )
2445 {
2446 MsqPrintError err( std::cout );
2447 PatchData pd;
2448 if (mSettings)
2449 pd.attach_settings(mSettings);
2450 std::vector<size_t> handles, indices;
2451 std::vector<Vector3D> grad;
2452 std::vector<Matrix3D> Hess;
2453 double qmval;
2454 bool rval;
2455
2456 CPPUNIT_ASSERT( !types.empty() );
2457 for (size_t i = 0; i < types.size(); ++i) {
2458 // get a patch
2459 get_ideal_element( types[i], false, pd );
2460 qm->get_evaluations( pd, handles, false, err );
2461 CPPUNIT_ASSERT( !MSQ_CHKERR(err) );
2462 CPPUNIT_ASSERT( !handles.empty() );
2463 for (size_t j = 0; j < handles.size(); ++j) {
2464 // evaluate
2465 rval = qm->evaluate_with_Hessian( pd, handles[j], qmval, indices, grad, Hess, err );
2466 CPPUNIT_ASSERT(!MSQ_CHKERR(err));
2467 CPPUNIT_ASSERT( rval );
2468 const size_t n = indices.size();
2469 CPPUNIT_ASSERT_EQUAL( n, grad.size() );
2470 CPPUNIT_ASSERT_EQUAL( n*(n+1)/2, Hess.size() );
2471 // check that diagonal Hessian blocks are symmetric
2472 size_t h_idx = 0;
2473 for (size_t k = 0; k < n; h_idx+=(n-k), ++k) {
2474 CPPUNIT_ASSERT_DOUBLES_EQUAL( Hess[h_idx][0][1], Hess[h_idx][1][0], 5e-3 );
2475 CPPUNIT_ASSERT_DOUBLES_EQUAL( Hess[h_idx][0][2], Hess[h_idx][2][0], 5e-3 );
2476 CPPUNIT_ASSERT_DOUBLES_EQUAL( Hess[h_idx][1][2], Hess[h_idx][2][1], 5e-3 );
2477 } // for(k < Hess.size())
2478 } // for(j < handles.size())
2479 } // for(i < types.size())
2480 }
2481
test_gradient_with_fixed_vertex(QualityMetric * qm,const Settings * settings)2482 void QualityMetricTester::test_gradient_with_fixed_vertex( QualityMetric* qm,
2483 const Settings* settings )
2484 {
2485 CPPUNIT_ASSERT(!types.empty());
2486 for (unsigned i = 0; i < types.size(); ++i)
2487 test_gradient_with_fixed_vertex( types[i], qm, settings );
2488 }
2489
test_hessian_with_fixed_vertex(QualityMetric * qm,const Settings * settings)2490 void QualityMetricTester::test_hessian_with_fixed_vertex( QualityMetric* qm,
2491 const Settings* settings )
2492 {
2493 CPPUNIT_ASSERT(!types.empty());
2494 for (unsigned i = 0; i < types.size(); ++i)
2495 test_hessian_with_fixed_vertex( types[i], qm, settings );
2496 }
2497
test_diagonal_with_fixed_vertex(QualityMetric * qm,const Settings * settings)2498 void QualityMetricTester::test_diagonal_with_fixed_vertex( QualityMetric* qm,
2499 const Settings* settings )
2500 {
2501 CPPUNIT_ASSERT(!types.empty());
2502 for (unsigned i = 0; i < types.size(); ++i)
2503 test_diagonal_with_fixed_vertex( types[i], qm, settings );
2504 }
2505
test_gradient_with_fixed_vertex(EntityTopology type,QualityMetric * qm,const Settings * settings)2506 void QualityMetricTester::test_gradient_with_fixed_vertex( EntityTopology type,
2507 QualityMetric* qm,
2508 const Settings* settings )
2509 {
2510 MsqPrintError err( std::cout );
2511 PatchData pd;
2512 if (settings)
2513 pd.attach_settings(settings);
2514 else if (mSettings)
2515 pd.attach_settings(mSettings);
2516
2517 std::vector<size_t> handles, indices, indices_fixed, conn;
2518 std::vector<Vector3D> grad, grad_fixed;
2519 double val, val_fixed;
2520 bool rval;
2521 const int n = TopologyInfo::corners( type );
2522
2523 get_nonideal_element( type, pd );
2524 pd.element_by_index(0).get_vertex_indices( conn );
2525
2526 qm->get_evaluations( pd, handles, false, err );
2527 ASSERT_NO_ERROR(err);
2528 // For sample-based metrics, it is difficult to set up such that
2529 // both surface and volume elements have only one sample point.
2530 // Make things easier by skipping element types with no sample points.
2531 if (handles.empty())
2532 return;
2533 CPPUNIT_ASSERT_EQUAL_MESSAGE( "test only valid for element-based metrics",
2534 (size_t)1, handles.size() );
2535 size_t h = handles[0];
2536
2537 rval = qm->evaluate_with_gradient( pd, h, val, indices, grad, err );
2538 ASSERT_NO_ERROR(err);
2539 CPPUNIT_ASSERT(rval);
2540 CPPUNIT_ASSERT_EQUAL( (size_t)n, indices.size() );
2541 CPPUNIT_ASSERT_EQUAL( (size_t)n, grad.size() );
2542
2543 for (int i = 0; i < n; ++i) {
2544 get_nonideal_element( type, pd, i );
2545
2546 qm->get_evaluations( pd, handles, false, err );
2547 ASSERT_NO_ERROR(err);
2548 CPPUNIT_ASSERT_EQUAL_MESSAGE( "test only valid for element-based metrics",
2549 (size_t)1, handles.size() );
2550 h = handles[0];
2551
2552 rval = qm->evaluate_with_gradient( pd, h, val_fixed, indices_fixed, grad_fixed, err );
2553 ASSERT_NO_ERROR(err);
2554 CPPUNIT_ASSERT(rval);
2555 CPPUNIT_ASSERT_EQUAL( (size_t)1, indices_fixed.size() );
2556 CPPUNIT_ASSERT_EQUAL( (size_t)1, grad_fixed.size() );
2557
2558 CPPUNIT_ASSERT_DOUBLES_EQUAL( val, val_fixed, 1e-5 );
2559 CPPUNIT_ASSERT_VECTORS_EQUAL( grad[conn[i]], grad_fixed.front(), 1e-5 );
2560 }
2561 }
2562
2563
test_hessian_with_fixed_vertex(EntityTopology type,QualityMetric * qm,const Settings * settings)2564 void QualityMetricTester::test_hessian_with_fixed_vertex( EntityTopology type,
2565 QualityMetric* qm,
2566 const Settings* settings )
2567 {
2568 MsqPrintError err( std::cout );
2569 PatchData pd;
2570 if (settings)
2571 pd.attach_settings(settings);
2572 else if (mSettings)
2573 pd.attach_settings(mSettings);
2574 std::vector<size_t> handles, indices, indices_fixed, conn;
2575 std::vector<Vector3D> grad, grad_fixed;
2576 std::vector<Matrix3D> hess, hess_fixed;
2577 double val, val_fixed;
2578 bool rval;
2579 const int n = TopologyInfo::corners( type );
2580
2581 get_nonideal_element( type, pd );
2582 pd.element_by_index(0).get_vertex_indices( conn );
2583
2584 qm->get_evaluations( pd, handles, false, err );
2585 ASSERT_NO_ERROR(err);
2586 // For sample-based metrics, it is difficult to set up such that
2587 // both surface and volume elements have only one sample point.
2588 // Make things easier by skipping element types with no sample points.
2589 if (handles.empty())
2590 return;
2591 CPPUNIT_ASSERT_EQUAL_MESSAGE( "test only valid for element-based metrics",
2592 (size_t)1, handles.size() );
2593 size_t h = handles[0];
2594
2595 rval = qm->evaluate_with_Hessian( pd, h, val, indices, grad, hess, err );
2596 ASSERT_NO_ERROR(err);
2597 CPPUNIT_ASSERT(rval);
2598 CPPUNIT_ASSERT_EQUAL( (size_t)n, indices.size() );
2599 CPPUNIT_ASSERT_EQUAL( (size_t)n, grad.size() );
2600 CPPUNIT_ASSERT_EQUAL( (size_t)(n*(n+1)/2), hess.size() );
2601
2602 for (int i = 0; i < n; ++i) {
2603 get_nonideal_element( type, pd, i );
2604
2605 qm->get_evaluations( pd, handles, false, err );
2606 ASSERT_NO_ERROR(err);
2607 CPPUNIT_ASSERT_EQUAL_MESSAGE( "test only valid for element-based metrics",
2608 (size_t)1, handles.size() );
2609 h = handles[0];
2610
2611 rval = qm->evaluate_with_Hessian( pd, h, val_fixed, indices_fixed, grad_fixed, hess_fixed, err );
2612 ASSERT_NO_ERROR(err);
2613 CPPUNIT_ASSERT(rval);
2614 CPPUNIT_ASSERT_EQUAL( (size_t)1, indices_fixed.size() );
2615 CPPUNIT_ASSERT_EQUAL( (size_t)1, grad_fixed.size() );
2616 CPPUNIT_ASSERT_EQUAL( (size_t)1, hess_fixed.size() );
2617
2618 const size_t j = conn[i];
2619 CPPUNIT_ASSERT_DOUBLES_EQUAL( val, val_fixed, 1e-5 );
2620 CPPUNIT_ASSERT_VECTORS_EQUAL( grad[j], grad_fixed.front(), 1e-5 );
2621 CPPUNIT_ASSERT_MATRICES_EQUAL( hess[n*j - j*(j-1)/2], hess_fixed.front(), 1e-4 );
2622 }
2623 }
2624
test_diagonal_with_fixed_vertex(EntityTopology type,QualityMetric * qm,const Settings * settings)2625 void QualityMetricTester::test_diagonal_with_fixed_vertex( EntityTopology type,
2626 QualityMetric* qm,
2627 const Settings* settings )
2628 {
2629 MsqPrintError err( std::cout );
2630 PatchData pd;
2631 if (settings)
2632 pd.attach_settings(settings);
2633 else if (mSettings)
2634 pd.attach_settings(mSettings);
2635 std::vector<size_t> handles, indices, indices_fixed, conn;
2636 std::vector<Vector3D> grad, grad_fixed;
2637 std::vector<SymMatrix3D> hess, hess_fixed;
2638 double val, val_fixed;
2639 bool rval;
2640 const int n = TopologyInfo::corners( type );
2641
2642 get_nonideal_element( type, pd );
2643 pd.element_by_index(0).get_vertex_indices( conn );
2644
2645 qm->get_evaluations( pd, handles, false, err );
2646 ASSERT_NO_ERROR(err);
2647 // For sample-based metrics, it is difficult to set up such that
2648 // both surface and volume elements have only one sample point.
2649 // Make things easier by skipping element types with no sample points.
2650 if (handles.empty())
2651 return;
2652 CPPUNIT_ASSERT_EQUAL_MESSAGE( "test only valid for element-based metrics",
2653 (size_t)1, handles.size() );
2654 size_t h = handles[0];
2655
2656 rval = qm->evaluate_with_Hessian_diagonal( pd, h, val, indices, grad, hess, err );
2657 ASSERT_NO_ERROR(err);
2658 CPPUNIT_ASSERT(rval);
2659 CPPUNIT_ASSERT_EQUAL( (size_t)n, indices.size() );
2660 CPPUNIT_ASSERT_EQUAL( (size_t)n, grad.size() );
2661 CPPUNIT_ASSERT_EQUAL( (size_t)n, hess.size() );
2662
2663 for (int i = 0; i < n; ++i) {
2664 get_nonideal_element( type, pd, i );
2665
2666 qm->get_evaluations( pd, handles, false, err );
2667 ASSERT_NO_ERROR(err);
2668 CPPUNIT_ASSERT_EQUAL_MESSAGE( "test only valid for element-based metrics",
2669 (size_t)1, handles.size() );
2670 h = handles[0];
2671
2672 rval = qm->evaluate_with_Hessian_diagonal( pd, h, val_fixed, indices_fixed, grad_fixed, hess_fixed, err );
2673 ASSERT_NO_ERROR(err);
2674 CPPUNIT_ASSERT(rval);
2675 CPPUNIT_ASSERT_EQUAL( (size_t)1, indices_fixed.size() );
2676 CPPUNIT_ASSERT_EQUAL( (size_t)1, grad_fixed.size() );
2677 CPPUNIT_ASSERT_EQUAL( (size_t)1, hess_fixed.size() );
2678
2679 const size_t j = conn[i];
2680 CPPUNIT_ASSERT_DOUBLES_EQUAL( val, val_fixed, 1e-5 );
2681 CPPUNIT_ASSERT_VECTORS_EQUAL( grad[j], grad_fixed.front(), 1e-5 );
2682 CPPUNIT_ASSERT_MATRICES_EQUAL( hess[j], hess_fixed.front(), 1e-4 );
2683 }
2684 }
2685