1 /* *****************************************************************
2 MESQUITE -- The Mesh Quality Improvement Toolkit
3
4 Copyright 2004 Sandia Corporation and Argonne National
5 Laboratory. Under the terms of Contract DE-AC04-94AL85000
6 with Sandia Corporation, the U.S. Government retains certain
7 rights in this software.
8
9 This library is free software; you can redistribute it and/or
10 modify it under the terms of the GNU Lesser General Public
11 License as published by the Free Software Foundation; either
12 version 2.1 of the License, or (at your option) any later version.
13
14 This library is distributed in the hope that it will be useful,
15 but WITHOUT ANY WARRANTY; without even the implied warranty of
16 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 Lesser General Public License for more details.
18
19 You should have received a copy of the GNU Lesser General Public License
20 (lgpl.txt) along with this library; if not, write to the Free Software
21 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
22
23 kraftche@cae.wisc.edu
24
25 ***************************************************************** */
26 #include <iostream>
27 #include <sstream>
28 using std::cout;
29 using std::cerr;
30 using std::endl;
31
32 #include <stdlib.h>
33
34 #include "Mesquite.hpp"
35 #include "MeshImpl.hpp"
36 #include "MsqError.hpp"
37 #include "InstructionQueue.hpp"
38 #include "TerminationCriterion.hpp"
39 #include "QualityAssessor.hpp"
40 #include "PlanarDomain.hpp"
41
42 #include "IdealWeightInverseMeanRatio.hpp"
43 #include "ElementPMeanP.hpp"
44 #include "TInverseMeanRatio.hpp"
45 #include "TQualityMetric.hpp"
46 #include "TOffset.hpp"
47 #include "TOffset.hpp"
48
49 #include "IdealShapeTarget.hpp"
50 #include "PMeanPTemplate.hpp"
51 #include "SteepestDescent.hpp"
52 #include "ConjugateGradient.hpp"
53 #include "QuasiNewton.hpp"
54
55 #include "MeshWriter.hpp"
56 #include "MsqTimer.hpp"
57 #include "TestUtil.hpp"
58
59 using namespace MBMesquite;
60
61 std::string DEFAULT_INPUT = TestDir + "/2D/vtk/tris/untangled/equil_tri2.vtk";
62
63 /* Print usage or help information: exits if err == true */
usage(const char * argv0=0,bool err=true)64 void usage( const char* argv0 = 0, bool err = true )
65 {
66 std::ostream& s = err ? cerr : cout;
67 const char* defname = "main";
68 if (!argv0)
69 argv0 = defname;
70
71 s << "Usage: " << defname << " [-n|-c|-q] [-e] [-t] [-a] [-N] [{-v|-g|-p} <output_file>] [<input_file>]" << endl
72 << " " << defname << " -h" << std::endl;
73
74 if (err)
75 exit(1);
76
77 s << " -n : Use SteepestDescent solver (default)" << endl
78 << " -c : Use ConjugateGradient solver" << endl
79 << " -q : Use QuasiNewton solver" << endl
80 << " -e : Test IdealWeightInverseMeanRatio metric" << endl
81 << " -t : Test InverseMeanRatio3D target metric" << endl
82 << " -a : Test ElementPMeanP(InverseMeanRatio3D)" << endl
83 << " -N : Test InverseMeanRatio3D with finite difference derivatives" << endl
84 << " -C : Compare InverseMeanRatio3D and IdealWeightInverseMeanRatio" << endl
85 << " -v : Write output mesh to specified VTK file" << endl
86 << " -g : Write output mesh to specified GnuPlot file" << endl
87 << " -p : Write output mesh to specified EPS file" << endl
88 << " -P : Write solver plot data to specified file." << endl
89 << "Default is all metrics if non specified." << endl
90 << "Default input file: " << DEFAULT_INPUT << endl;
91 }
92
93 const char* vtk_file = 0; /* vtk output file name */
94 const char* eps_file = 0; /* eps output file name */
95 const char* gpt_file = 0; /* GNUPlot output file name */
96 const char* plot_file = 0; /* Time-dependent plot of solver data */
97
98 enum Solver { STEEP_DESCENT, CONJ_GRAD, QUASI_NEWT };
99
100 /* Run an optimization: returns average quality of final mesh */
101 double run( QualityMetric* metric,
102 Solver solver,
103 const char* input_file,
104 double& seconds_out,
105 int& iterations_out );
106
107 /* Force finite difference approximation of derivatives for
108 * an existing quality metric */
109 class NumericQM : public QualityMetric {
110 private:
111 QualityMetric* realMetric;
112 public:
NumericQM(QualityMetric * real_metric)113 NumericQM( QualityMetric* real_metric ) : realMetric(real_metric) {}
114 virtual MetricType get_metric_type() const;
115 virtual std::string get_name() const;
116 virtual int get_negate_flag() const;
117 virtual void get_evaluations( PatchData& pd,
118 std::vector<size_t>& handles,
119 bool free_vertices_only,
120 MsqError& err );
121 virtual bool evaluate( PatchData& pd,
122 size_t handle,
123 double& value,
124 MsqError& err );
125 virtual bool evaluate_with_indices( PatchData& pd,
126 size_t handle,
127 double& value,
128 std::vector<size_t>& indices,
129 MsqError& err );
130 };
get_metric_type() const131 QualityMetric::MetricType NumericQM::get_metric_type() const
132 { return realMetric->get_metric_type(); }
get_name() const133 std::string NumericQM::get_name() const{
134 std::string r = realMetric->get_name();
135 r += " (FD)";
136 return r;
137 }
138
139 /* At each evaluation of this metric, compare the values resulting
140 * from the evaluation of two other metrics: flag an error if they
141 * differ or return the common result if the same */
142 class CompareMetric : public QualityMetric
143 {
144 private:
145 QualityMetric *metric1, *metric2;
146 bool maskPlane;
147 int maskAxis;
148 std::vector<size_t> m2Handles;
149 std::vector<Vector3D> m2Grad;
150 std::vector<SymMatrix3D> m2Diag;
151 std::vector<Matrix3D> m2Hess;
152 static const double epsilon;
153 public:
CompareMetric(QualityMetric * qm1,QualityMetric * qm2,bool mask_qm2_coord=false)154 CompareMetric( QualityMetric* qm1,
155 QualityMetric* qm2,
156 bool mask_qm2_coord = false )
157 : metric1(qm1), metric2(qm2), maskPlane(mask_qm2_coord), maskAxis(-1)
158 {}
159
160 MetricType get_metric_type() const;
161
162 std::string get_name() const;
163
164 int get_negate_flag() const;
165
166 void get_evaluations( PatchData& pd,
167 std::vector<size_t>& handles,
168 bool free_vertices_only,
169 MsqError& err );
170
171 bool evaluate( PatchData& pd,
172 size_t handle,
173 double& value,
174 MsqError& err );
175
176 bool evaluate_with_indices( PatchData& pd,
177 size_t handle,
178 double& value,
179 std::vector<size_t>& indices,
180 MsqError& err );
181
182 bool evaluate_with_gradient( PatchData& pd,
183 size_t handle,
184 double& value,
185 std::vector<size_t>& indices,
186 std::vector<Vector3D>& gradient,
187 MsqError& err );
188
189 bool evaluate_with_Hessian_diagonal( PatchData& pd,
190 size_t handle,
191 double& value,
192 std::vector<size_t>& indices,
193 std::vector<Vector3D>& gradient,
194 std::vector<SymMatrix3D>& Hessian_diagonal,
195 MsqError& err );
196
197 bool evaluate_with_Hessian( PatchData& pd,
198 size_t handle,
199 double& value,
200 std::vector<size_t>& indices,
201 std::vector<Vector3D>& gradient,
202 std::vector<Matrix3D>& Hessian,
203 MsqError& err );
204
205 void get_mask_axis( PatchData& pd );
206 bool equal( Vector3D grad1, const Vector3D& grad2 ) const;
207 bool equal( SymMatrix3D hess1, const SymMatrix3D& hess2 ) const;
208 bool equal( Matrix3D hess1, const Matrix3D& hess2 ) const;
209 };
210 const double CompareMetric::epsilon = 5e-2;
211
212 /* Parse command line options and call 'run' */
main(int argc,char * argv[])213 int main( int argc, char* argv[] )
214 {
215 Solver solver = STEEP_DESCENT;
216 bool do_non_target_metric = false;
217 bool do_new_target_metric = false;
218 bool do_new_target_average = false;
219 bool do_new_target_numeric = false;
220 bool do_compare_metric = false;
221 const char* input_file = DEFAULT_INPUT.c_str();
222 bool no_more_flags = false;
223
224 std::list<const char**> exp_list;
225 for (int i = 1; i < argc; ++i) {
226 if (argv[i][0] == '-' && !no_more_flags) {
227 for (int k = 1; argv[i][k]; ++k) {
228 switch (argv[i][k]) {
229 case 'n': solver = STEEP_DESCENT; break;
230 case 'c': solver = CONJ_GRAD; break;
231 case 'q': solver = QUASI_NEWT; break;
232 case 'e': do_non_target_metric = true; break;
233 case 't': do_new_target_metric = true; break;
234 case 'a': do_new_target_average = true; break;
235 case 'N': do_new_target_numeric = true; break;
236 case 'C': do_compare_metric = true; break;
237 case 'v': exp_list.push_back(&vtk_file); break;
238 case 'g': exp_list.push_back(&gpt_file); break;
239 case 'p': exp_list.push_back(&eps_file); break;
240 case 'P': exp_list.push_back(&plot_file); break;
241 case 'h': usage(argv[0],false); return 0;
242 case '-': no_more_flags = true; break;
243 default:
244 cerr << "Invalid flag: '" << argv[i][k] << "'" << endl;
245 usage(argv[0]);
246 }
247 }
248 }
249 else if (!exp_list.empty()) {
250 const char** ptr = exp_list.front();
251 exp_list.pop_front();
252 *ptr = argv[i];
253 }
254 else {
255 if (!DEFAULT_INPUT.compare(input_file)) {
256 cerr << "Unexpected argument: \"" << argv[i] << '"' << endl;
257 usage(argv[0]);
258 }
259 input_file = argv[i];
260 }
261 }
262
263 int count = 0;
264 if (do_non_target_metric)
265 ++count;
266 if (do_new_target_metric)
267 ++count;
268 if (do_new_target_average)
269 ++count;
270 if (do_new_target_numeric)
271 ++count;
272 if (do_compare_metric)
273 ++count;
274
275 if (!count) {
276 do_compare_metric = true;
277 count = 1;
278 }
279
280 if ((vtk_file || gpt_file || eps_file) && count != 1) {
281 cerr << "Error: Cannot write output file if running multiple tests" << endl;
282 return 2;
283 }
284
285 IdealWeightInverseMeanRatio non_target_metric;
286 IdealShapeTarget new_target;
287 TInverseMeanRatio tmp;
288 TOffset tmp_off( 1.0, &tmp ); // target metrics are zero-based
289 TQualityMetric new_target_metric( &new_target, &tmp_off );
290 ElementPMeanP new_target_average( 1.0, &new_target_metric );
291 NumericQM new_target_numeric( &new_target_metric );
292 CompareMetric comp_metric( &non_target_metric, &new_target_average, true );
293
294 std::ostringstream os;
295 double secs,qual;
296 if (do_non_target_metric) {
297 qual = run( &non_target_metric, solver, input_file, secs, count );
298 os << "IdealWeightInverseMeanRatio: " << qual << " after " << count << " iterations in " << secs << " seconds" << endl;
299 }
300 if (do_new_target_metric) {
301 qual = run( &new_target_metric, solver, input_file, secs, count );
302 os << "TQualityMetric : " << qual << " after " << count << " iterations in " << secs << " seconds" << endl;
303 }
304 if (do_new_target_average) {
305 qual = run( &new_target_average, solver, input_file, secs, count );
306 os << "ElementPMeanP : " << qual << " after " << count << " iterations in " << secs << " seconds" << endl;
307 }
308 if (do_new_target_numeric) {
309 qual = run( &new_target_numeric, solver, input_file, secs, count );
310 os << "TQualityMetric (FD) : " << qual << " after " << count << " iterations in " << secs << " seconds" << endl;
311 }
312 if (do_compare_metric) {
313 qual = run( &comp_metric, solver, input_file, secs, count );
314 os << "Metric comparison : " << qual << " after " << count << " iterations in " << secs << " seconds" << endl;
315 }
316
317 cout << endl << os.str() << endl;
318 return 0;
319 }
320
321
322
run(QualityMetric * metric,Solver solver_type,const char * input_file,double & seconds_out,int & iterations_out)323 double run( QualityMetric* metric,
324 Solver solver_type,
325 const char* input_file,
326 double& seconds_out,
327 int& iterations_out )
328 {
329 MsqPrintError err(cerr);
330 IdealWeightInverseMeanRatio qa_metric;
331 TerminationCriterion inner, outer;
332 outer.add_iteration_limit( 1 );
333 inner.add_absolute_vertex_movement( 1e-4 );
334 inner.add_iteration_limit( 100 );
335 PMeanPTemplate of( 1.0, metric );
336 QualityAssessor qa( &qa_metric );
337 qa.add_quality_assessment( metric );
338 InstructionQueue q;
339 SteepestDescent steep(&of);
340 QuasiNewton quasi(&of);
341 ConjugateGradient conj(&of);
342 VertexMover* solver = 0;
343 switch (solver_type) {
344 case STEEP_DESCENT: solver = &steep; break;
345 case QUASI_NEWT:solver = &quasi;break;
346 case CONJ_GRAD: solver = &conj; break;
347 }
348 q.set_master_quality_improver( solver, err );
349 q.add_quality_assessor( &qa, err );
350 solver->set_inner_termination_criterion(&inner);
351 solver->set_outer_termination_criterion(&outer);
352
353
354 if (plot_file)
355 inner.write_iterations( plot_file, err );
356
357 MeshImpl mesh;
358 mesh.read_vtk( input_file, err );
359 if (err) {
360 cerr << "Failed to read input file: \"" << input_file << '"' << endl;
361 exit(1);
362 }
363
364 std::vector<Mesh::VertexHandle> handles;
365 mesh.get_all_vertices( handles, err );
366 if (handles.empty()) {
367 cerr << "no veritces in mesh" << endl;
368 exit(1);
369 }
370 std::vector<MsqVertex> coords(handles.size());
371 mesh.vertices_get_coordinates( arrptr(handles), arrptr(coords), handles.size(), err );
372 Vector3D min(HUGE_VAL), max(-HUGE_VAL);
373 for (size_t i = 0; i < coords.size(); ++i) {
374 for (int j = 0; j < 3; ++j) {
375 if (coords[i][j] < min[j])
376 min[j] = coords[i][j];
377 if (coords[i][j] > max[j])
378 max[j] = coords[i][j];
379 }
380 }
381
382 Vector3D size = max - min;
383 PlanarDomain* domain = 0;
384 if (size[0] < 1e-4)
385 domain = new PlanarDomain( PlanarDomain::YZ, min[0] );
386 else if (size[1] < 1e-4)
387 domain = new PlanarDomain( PlanarDomain::XZ, min[1] );
388 else if (size[2] < 1e-4)
389 domain = new PlanarDomain( PlanarDomain::XY, min[2] );
390
391 Timer timer;
392 MeshDomainAssoc mesh_and_domain = MeshDomainAssoc(&mesh, domain);
393 q.run_instructions( &mesh_and_domain, err );
394 seconds_out = timer.since_birth();
395 if (err) {
396 cerr << "Optimization failed." << endl << err << endl;
397 abort();
398 }
399
400 if (vtk_file) {
401 MeshWriter::write_vtk( &mesh, vtk_file, err );
402 if (err)
403 cerr << vtk_file << ": failed to write file." << endl;
404 }
405 if (gpt_file) {
406 MeshWriter::write_gnuplot( &mesh, gpt_file, err );
407 if (err)
408 cerr << gpt_file << ": failed to write file." << endl;
409 }
410 if (eps_file) {
411 PlanarDomain xy(PlanarDomain::XY);
412 MeshWriter::Projection proj( domain ? domain : &xy );
413 MeshWriter::write_eps( &mesh, eps_file, proj, err );
414 if (err)
415 cerr << eps_file << ": failed to write file." << endl;
416 }
417 delete domain;
418
419 iterations_out = inner.get_iteration_count();
420
421 const QualityAssessor::Assessor* a = qa.get_results( &qa_metric );
422 return a->get_average();
423 }
424
get_negate_flag() const425 int NumericQM::get_negate_flag() const
426 { return realMetric->get_negate_flag(); }
427
get_evaluations(PatchData & pd,std::vector<size_t> & handles,bool free_vertices_only,MsqError & err)428 void NumericQM::get_evaluations( PatchData& pd,
429 std::vector<size_t>& handles,
430 bool free_vertices_only,
431 MsqError& err )
432 { realMetric->get_evaluations( pd, handles, free_vertices_only, err ); }
433
evaluate(PatchData & pd,size_t handle,double & value,MsqError & err)434 bool NumericQM::evaluate( PatchData& pd,
435 size_t handle,
436 double& value,
437 MsqError& err )
438 { return realMetric->evaluate( pd, handle, value, err ); }
439
evaluate_with_indices(PatchData & pd,size_t handle,double & value,std::vector<size_t> & indices,MsqError & err)440 bool NumericQM::evaluate_with_indices( PatchData& pd,
441 size_t handle,
442 double& value,
443 std::vector<size_t>& indices,
444 MsqError& err )
445 { return realMetric->evaluate_with_indices( pd, handle, value, indices, err ); }
446
447
get_metric_type() const448 QualityMetric::MetricType CompareMetric::get_metric_type() const
449 {
450 MetricType t1 = metric1->get_metric_type();
451 assert(metric2->get_metric_type() == t1);
452 return t1;
453 }
454
get_name() const455 std::string CompareMetric::get_name() const
456 {
457 std::string n = metric1->get_name();
458 n += " =? ";
459 n += metric2->get_name();
460 return n;
461 }
462
get_negate_flag() const463 int CompareMetric::get_negate_flag() const
464 {
465 assert(metric1->get_negate_flag() == metric2->get_negate_flag());
466 return metric1->get_negate_flag();
467 }
468
get_evaluations(PatchData & pd,std::vector<size_t> & handles,bool free_vertices_only,MsqError & err)469 void CompareMetric::get_evaluations( PatchData& pd,
470 std::vector<size_t>& handles,
471 bool free_vertices_only,
472 MsqError& err )
473 {
474 if (maskPlane)
475 get_mask_axis(pd);
476
477 m2Handles.clear();
478 metric1->get_evaluations( pd, handles, free_vertices_only, err ); MSQ_ERRRTN(err);
479 metric2->get_evaluations( pd, m2Handles, free_vertices_only, err ); MSQ_ERRRTN(err);
480 bool same = (handles.size() == m2Handles.size());
481 std::sort( m2Handles.begin(), m2Handles.end() );
482 for (std::vector<size_t>::iterator i = handles.begin(); i != handles.end(); ++i)
483 if (!std::binary_search( m2Handles.begin(), m2Handles.end(), *i ))
484 same = false;
485 if (!same) {
486 MSQ_SETERR(err)("Metrics have incompatible lists of evaluation handles.\n",
487 MsqError::INVALID_STATE);
488 }
489 }
490
491
evaluate(PatchData & pd,size_t handle,double & value,MsqError & err)492 bool CompareMetric::evaluate( PatchData& pd,
493 size_t handle,
494 double& value,
495 MsqError& err )
496 {
497 double m2val;
498 bool r1, r2;
499 r1 = metric1->evaluate( pd, handle, value, err ); MSQ_ERRZERO(err);
500 r2 = metric2->evaluate( pd, handle, m2val, err ); MSQ_ERRZERO(err);
501 if (r1 != r2 || (r1 && fabs(value - m2val) > epsilon)) {
502 MSQ_SETERR(err)(MsqError::INVALID_STATE,
503 "Metrics returned different values for handle %lu "
504 "in evaluate:\n"
505 "\t%s %f vs. %s %f\n", (unsigned long)handle,
506 r1?"true":"false",value,r2?"true":"false",m2val);
507 }
508
509 return r1 && !err;
510 }
511
512
evaluate_with_indices(PatchData & pd,size_t handle,double & value,std::vector<size_t> & indices,MsqError & err)513 bool CompareMetric::evaluate_with_indices( PatchData& pd,
514 size_t handle,
515 double& value,
516 std::vector<size_t>& indices,
517 MsqError& err )
518 {
519 double m2val;
520 bool r1, r2;
521 m2Handles.clear();
522 r1 = metric1->evaluate_with_indices( pd, handle, value, indices, err ); MSQ_ERRZERO(err);
523 r2 = metric2->evaluate_with_indices( pd, handle, m2val, m2Handles, err ); MSQ_ERRZERO(err);
524 if (r1 != r2 || (r1 && fabs(value - m2val) > epsilon)) {
525 MSQ_SETERR(err)(MsqError::INVALID_STATE,
526 "Metrics returned different values for handle %lu "
527 "in evaluate_with_indices:\n"
528 "\t%s %f vs. %s %f\n", (unsigned long)handle,
529 r1?"true":"false",value,r2?"true":"false",m2val);
530 }
531 else {
532 bool same = (indices.size() == m2Handles.size());
533 std::sort( m2Handles.begin(), m2Handles.end() );
534 for (std::vector<size_t>::iterator i = indices.begin(); i != indices.end(); ++i)
535 if (!std::binary_search( m2Handles.begin(), m2Handles.end(), *i ))
536 same = false;
537 if (!same) {
538 MSQ_SETERR(err)(MsqError::INVALID_STATE,
539 "Metrics returned incompatible lists of vertex indices"
540 " for handle %lu in evaluate_with_indices\n.",
541 (unsigned long)handle );
542 }
543 }
544
545 return r1 && !err;
546 }
547
548
evaluate_with_gradient(PatchData & pd,size_t handle,double & value,std::vector<size_t> & indices,std::vector<Vector3D> & gradient,MsqError & err)549 bool CompareMetric::evaluate_with_gradient( PatchData& pd,
550 size_t handle,
551 double& value,
552 std::vector<size_t>& indices,
553 std::vector<Vector3D>& gradient,
554 MsqError& err )
555 {
556 double m2val;
557 bool r1, r2;
558 m2Handles.clear();
559 m2Grad.clear();
560 r1 = metric1->evaluate_with_gradient( pd, handle, value, indices, gradient, err ); MSQ_ERRZERO(err);
561 r2 = metric2->evaluate_with_gradient( pd, handle, m2val, m2Handles, m2Grad, err ); MSQ_ERRZERO(err);
562 if (r1 != r2 || (r1 && fabs(value - m2val) > epsilon)) {
563 MSQ_SETERR(err)(MsqError::INVALID_STATE,
564 "Metrics returned different values for handle %lu in "
565 "evaluate_with_gradient:\n"
566 "\t%s %f vs. %s %f\n", (unsigned long)handle,
567 r1?"true":"false",value,r2?"true":"false",m2val);
568 }
569 else {
570 std::vector<size_t>::const_iterator i, j;
571 std::vector<Vector3D>::const_iterator r, s;
572 int grad_diff = 0;
573 bool same = (indices.size() == m2Handles.size());
574 std::sort( m2Handles.begin(), m2Handles.end() );
575 for (i = indices.begin(); i != indices.end(); ++i) {
576 j = std::lower_bound( m2Handles.begin(), m2Handles.end(), *i );
577 if (j == m2Handles.end() || *j != *i) {
578 same = false;
579 continue;
580 }
581
582 r = gradient.begin() + (i - indices.begin());
583 s = m2Grad.begin() + (j - m2Handles.begin());
584 if (!equal(*r,*s))
585 ++grad_diff;
586 }
587
588 if (!same) {
589 MSQ_SETERR(err)(MsqError::INVALID_STATE,
590 "Metrics returned incompatible lists of vertex indices"
591 " for handle %lu in evaluate_with_gradient\n.",
592 (unsigned long)handle );
593 }
594 else if (grad_diff) {
595 MSQ_SETERR(err)(MsqError::INVALID_STATE,
596 "Metrics returned different gradient vectors for "
597 " %d of %u vertices for handle %lu in "
598 "evaluate_with_gradient\n.",
599 grad_diff, (unsigned)gradient.size(),
600 (unsigned long)handle );
601 }
602 }
603
604 return r1 && !err;
605 }
606
evaluate_with_Hessian_diagonal(PatchData & pd,size_t handle,double & value,std::vector<size_t> & indices,std::vector<Vector3D> & gradient,std::vector<SymMatrix3D> & diagonal,MsqError & err)607 bool CompareMetric::evaluate_with_Hessian_diagonal( PatchData& pd,
608 size_t handle,
609 double& value,
610 std::vector<size_t>& indices,
611 std::vector<Vector3D>& gradient,
612 std::vector<SymMatrix3D>& diagonal,
613 MsqError& err )
614 {
615 double m2val;
616 bool r1, r2;
617 m2Handles.clear();
618 m2Grad.clear();
619 m2Diag.clear();
620 r1 = metric1->evaluate_with_Hessian_diagonal( pd, handle, value, indices, gradient, diagonal, err ); MSQ_ERRZERO(err);
621 r2 = metric2->evaluate_with_Hessian_diagonal( pd, handle, m2val, m2Handles, m2Grad, m2Diag, err ); MSQ_ERRZERO(err);
622 if (r1 != r2 || (r1 && fabs(value - m2val) > epsilon)) {
623 MSQ_SETERR(err)(MsqError::INVALID_STATE,
624 "Metrics returned different values for handle %lu in "
625 "evaluate_with_Hessian_diagonal:\n"
626 "\t%s %f vs. %s %f\n", (unsigned long)handle,
627 r1?"true":"false",value,r2?"true":"false",m2val);
628 }
629 else {
630 std::vector<size_t>::const_iterator i, j;
631 std::vector<Vector3D>::const_iterator r, s;
632 std::vector<SymMatrix3D>::const_iterator u, v;
633 int grad_diff = 0, hess_diff = 0;
634 bool same = (indices.size() == m2Handles.size());
635 std::sort( m2Handles.begin(), m2Handles.end() );
636 for (i = indices.begin(); i != indices.end(); ++i) {
637 j = std::lower_bound( m2Handles.begin(), m2Handles.end(), *i );
638 if (j == m2Handles.end() || *j != *i) {
639 same = false;
640 continue;
641 }
642
643 r = gradient.begin() + (i - indices.begin());
644 s = m2Grad.begin() + (j - m2Handles.begin());
645 if (!equal(*r,*s))
646 ++grad_diff;
647
648 u = diagonal.begin() + (i - indices.begin());
649 v = m2Diag.begin() + (j - m2Handles.begin());
650 if (!equal(*u,*v))
651 ++hess_diff;
652 }
653
654 if (!same) {
655 MSQ_SETERR(err)(MsqError::INVALID_STATE,
656 "Metrics returned incompatible lists of vertex indices"
657 " for handle %lu in evaluate_with_Hessian_diagonal\n.",
658 (unsigned long)handle );
659 }
660 else if (grad_diff) {
661 MSQ_SETERR(err)(MsqError::INVALID_STATE,
662 "Metrics returned different gradient vectors for "
663 " %d of %u vertices for handle %lu in "
664 "evaluate_with_Hessian_diagonal\n.",
665 grad_diff, (unsigned)gradient.size(),
666 (unsigned long)handle );
667 }
668 else if (hess_diff) {
669 MSQ_SETERR(err)(MsqError::INVALID_STATE,
670 "Metrics returned different Hessian blocks for "
671 " %d of %u vertices for handle %lu in "
672 "evaluate_with_Hessian_diagonal\n.",
673 hess_diff, (unsigned)diagonal.size(),
674 (unsigned long)handle );
675 }
676 }
677
678 return r1 && !err;
679 }
680
evaluate_with_Hessian(PatchData & pd,size_t handle,double & value,std::vector<size_t> & indices,std::vector<Vector3D> & gradient,std::vector<Matrix3D> & Hessian,MsqError & err)681 bool CompareMetric::evaluate_with_Hessian( PatchData& pd,
682 size_t handle,
683 double& value,
684 std::vector<size_t>& indices,
685 std::vector<Vector3D>& gradient,
686 std::vector<Matrix3D>& Hessian,
687 MsqError& err )
688 {
689 double m2val;
690 bool r1, r2;
691 m2Handles.clear();
692 m2Grad.clear();
693 m2Hess.clear();
694 r1 = metric1->evaluate_with_Hessian( pd, handle, value, indices, gradient, Hessian, err ); MSQ_ERRZERO(err);
695 r2 = metric2->evaluate_with_Hessian( pd, handle, m2val, m2Handles, m2Grad, m2Hess, err ); MSQ_ERRZERO(err);
696 if (r1 != r2 || fabs(value - m2val) > epsilon) {
697 MSQ_SETERR(err)(MsqError::INVALID_STATE,
698 "Metrics returned different values for handle %lu in "
699 "evaluate_with_Hessian:\n"
700 "\t%s %f vs. %s %f\n", (unsigned long)handle,
701 r1?"true":"false",value,r2?"true":"false",m2val);
702 }
703 else {
704 std::vector<size_t>::const_iterator i, j;
705 std::vector<Vector3D>::const_iterator r, s;
706 int grad_diff = 0, hess_diff = 0;
707 bool same = (indices.size() == m2Handles.size());
708 std::sort( m2Handles.begin(), m2Handles.end() );
709 for (i = indices.begin(); i != indices.end(); ++i) {
710 j = std::lower_bound( m2Handles.begin(), m2Handles.end(), *i );
711 if (j == m2Handles.end() || *j != *i) {
712 same = false;
713 continue;
714 }
715
716 r = gradient.begin() + (i - indices.begin());
717 s = m2Grad.begin() + (j - m2Handles.begin());
718 if (!equal(*r,*s)) {
719 ++grad_diff;
720 // call again for so debugger can step into it after failure is found
721 std::vector<size_t> i2;
722 std::vector<Vector3D> g2;
723 std::vector<Matrix3D> h2;
724 metric2->evaluate_with_Hessian(pd, handle, m2val, i2, g2, h2, err );
725 }
726 }
727
728 if (!same) {
729 MSQ_SETERR(err)(MsqError::INVALID_STATE,
730 "Metrics returned incompatible lists of vertex indices"
731 " for handle %lu in evaluate_with_Hessian\n.",
732 (unsigned long)handle );
733 }
734 else if (grad_diff) {
735 MSQ_SETERR(err)(MsqError::INVALID_STATE,
736 "Metrics returned different gradient vectors for "
737 " %d of %u vertices for handle %lu in "
738 "evaluate_with_Hessian\n.",
739 grad_diff, (unsigned)gradient.size(),
740 (unsigned long)handle );
741 }
742 else {
743 size_t row, col, row2, col2, idx, idx2;
744 for (row = idx = 0; row < indices.size(); ++row) {
745 row2 = std::lower_bound( m2Handles.begin(), m2Handles.end(), indices[row] ) - m2Handles.begin();
746 for (col = row; col < indices.size(); ++col, ++idx) {
747 col2 = std::lower_bound( m2Handles.begin(), m2Handles.end(), indices[col] ) - m2Handles.begin();
748 if (row2 <= col2) {
749 idx2 = indices.size()*row2 - row2*(row2+1)/2 + col2;
750 if (!equal(Hessian[idx], m2Hess[idx2]))
751 ++hess_diff;
752 }
753 else {
754 idx2 = indices.size()*col2 - col2*(col2+1)/2 + row2;
755 if (!equal(Hessian[idx], transpose(m2Hess[idx2])))
756 ++hess_diff;
757 }
758 }
759 }
760
761 if (hess_diff) {
762 MSQ_SETERR(err)(MsqError::INVALID_STATE,
763 "Metrics returned different Hessian blocks for "
764 " %d of %u vertices for handle %lu in "
765 "evaluate_with_Hessian\n.",
766 hess_diff, (unsigned)Hessian.size(),
767 (unsigned long)handle );
768 }
769 }
770 }
771
772 return r1 && !err;
773 }
774
get_mask_axis(PatchData & pd)775 void CompareMetric::get_mask_axis( PatchData& pd )
776 {
777 maskAxis = -1;
778 PlanarDomain* dom = reinterpret_cast<PlanarDomain*>(pd.get_domain());
779 int bits = 0;
780 if (dom) {
781 Vector3D n = dom->get_normal();
782 for (int i = 0; i < 3; ++i)
783 if (fabs(n[i]) < epsilon)
784 bits |= (1 << i);
785 switch (bits) {
786 case 3: maskAxis = 2; break;
787 case 5: maskAxis = 1; break;
788 case 6: maskAxis = 0; break;
789 }
790 }
791 }
792
equal(Vector3D grad1,const Vector3D & grad2) const793 bool CompareMetric::equal( Vector3D grad1, const Vector3D& grad2 ) const
794 {
795 if (maskAxis >= 0)
796 grad1[maskAxis] = grad2[maskAxis];
797 return (grad1 - grad2).length_squared() <= epsilon*epsilon;
798 }
799
equal(SymMatrix3D hess1,const SymMatrix3D & hess2) const800 bool CompareMetric::equal( SymMatrix3D hess1, const SymMatrix3D& hess2 ) const
801 {
802 if (maskAxis >= 0) {
803 for (unsigned i = 0; i < 3; ++i) {
804 hess1(maskAxis,i) = hess2(maskAxis,i);
805 hess1(i,maskAxis) = hess2(i,maskAxis);
806 }
807 }
808 return Frobenius_2(hess1-hess2) <= epsilon*epsilon;
809 }
810
equal(Matrix3D hess1,const Matrix3D & hess2) const811 bool CompareMetric::equal( Matrix3D hess1, const Matrix3D& hess2 ) const
812 {
813 if (maskAxis >= 0) {
814 for (unsigned i = 0; i < 3; ++i) {
815 hess1[maskAxis][i] = hess2[maskAxis][i];
816 hess1[i][maskAxis] = hess2[i][maskAxis];
817 }
818 }
819 return Frobenius_2(hess1-hess2) <= epsilon*epsilon;
820 }
821