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 diachin2@llnl.gov, djmelan@sandia.gov, mbrewer@sandia.gov,
24 pknupp@sandia.gov, tleurent@mcs.anl.gov, tmunson@mcs.anl.gov
25
26 ***************************************************************** */
27 // -*- Mode : c++; tab-width: 2; c-tab-always-indent: t; indent-tabs-mode: nil; c-basic-offset: 2 -*-
28 //
29 // SUMMARY:
30 // USAGE:
31 //
32 // AUTHOR: Michael Brewer
33 // ORG: Sandia National Labs
34 // E-MAIL: mbrewer@sandia.gov
35 //
36 // ORIG-DATE: Jan. 29, 2003
37 // LAST-MOD: 25-Feb-04 at 10:49:32 by Thomas Leurent
38 //
39 // DESCRIPTION:
40 // ============
41 /*! \file PlanarGeometryTest.cpp
42
43 Regression testing using the planar geometry capabilities in
44 SimplifiedGeometryEngine.
45 */
46 // DESCRIP-END.
47 //
48
49 #include "meshfiles.h"
50
51 #include "PatchDataInstances.hpp"
52 #include "cppunit/extensions/HelperMacros.h"
53 #include <math.h>
54
55 #include "Mesquite.hpp"
56 #include "MsqError.hpp"
57 #include "InstructionQueue.hpp"
58 #include "PatchData.hpp"
59 //#include "StoppingCriterion.hpp"
60 #include "QualityAssessor.hpp"
61
62 #include "ConditionNumberQualityMetric.hpp"
63 #include "LPtoPTemplate.hpp"
64 #include "EdgeLengthQualityMetric.hpp"
65 #include "LaplacianSmoother.hpp"
66 #include "LInfTemplate.hpp"
67 #include "SteepestDescent.hpp"
68 #include "ConjugateGradient.hpp"
69 #include "AspectRatioGammaQualityMetric.hpp"
70 #include "UntangleBetaQualityMetric.hpp"
71 #include "MultiplyQualityMetric.hpp"
72 #include "PlanarDomain.hpp"
73
74 #include "MsqFPE.hpp"
75
76 #include "UnitUtil.hpp"
77
78 #include "MeshImpl.hpp"
79 #include <iostream>
80 using std::cout;
81 using std::endl;
82
83 using namespace MBMesquite;
84
85 class PlanarGeometryTest : public CppUnit::TestFixture
86 {
87 private:
88 CPPUNIT_TEST_SUITE(PlanarGeometryTest);
89 //run steepest descent on the tangled_tri.vtk mesh
90 CPPUNIT_TEST (test_plane_tri_tangled);
91 //run cg on tangled_quad.vtk mesh
92 CPPUNIT_TEST (test_plane_quad_tangled);
93 //run cg with asm metric on tri mesh in y=-5 plane
94 CPPUNIT_TEST (test_plane_tri_xz);
95 // test fit plane
96 CPPUNIT_TEST (test_fit_plane);
97 CPPUNIT_TEST_SUITE_END();
98
99 private:
100 double qualTol;//double used for double comparisons
101 int pF;//PRINT_FLAG
102 public:
setUp()103 void setUp()
104 {
105 //pF=1;//PRINT_FLAG IS ON
106 pF=0;//PRINT_FLAG IS OFF
107 //tolerance double
108 qualTol=MSQ_MIN;
109 }
110
tearDown()111 void tearDown()
112 {
113 }
114
115 public:
PlanarGeometryTest()116 PlanarGeometryTest()
117 {}
118
test_plane_tri_tangled()119 void test_plane_tri_tangled()
120 {
121 MBMesquite::MsqPrintError err(cout);
122 MBMesquite::MeshImpl mesh;
123
124 // This test doesn't use InstructionQueue, so
125 // we need to set up trapping of floating-point
126 // exceptions ourself.
127 MsqFPE fpe_trap( true );
128
129 mesh.read_vtk(MESH_FILES_DIR "2D/vtk/tris/tangled/tangled_tri.vtk", err);
130 CPPUNIT_ASSERT(!err);
131
132 //create geometry: plane z=5, normal (0,0,1)
133 Vector3D pnt(0,0,5);
134 Vector3D s_norm(0,0,1);
135 MBMesquite::PlanarDomain msq_geom(s_norm, pnt);
136
137 // creates an intruction queue
138 InstructionQueue queue1, queue2;
139
140 // creates a mean ratio quality metric ...
141 ConditionNumberQualityMetric shape;
142 UntangleBetaQualityMetric untan(.1);
143
144 // ... and builds an objective function with it (untangle)
145 LInfTemplate untan_func(&untan);
146 LPtoPTemplate shape_func(&shape,2,err);
147 //Make sure no errors
148 CPPUNIT_ASSERT(!err);
149 // creates the steepest descent optimization procedures
150 SteepestDescent pass1( &untan_func );
151 SteepestDescent pass2( &shape_func );
152 pass1.use_element_on_vertex_patch();
153 pass2.use_global_patch();
154 //Make sure no errors
155 CPPUNIT_ASSERT(!err);
156 QualityAssessor stop_qa=QualityAssessor( &untan );
157 QualityAssessor qa=QualityAssessor( &shape );
158 if(pF==0){
159 stop_qa.disable_printing_results();
160 qa.disable_printing_results();
161 }
162 //**********Set stopping criterion untangle ver small ********
163 //StoppingCriterion sc_qa(&stop_qa,-100,MSQ_MIN);
164 //pass1->set_stopping_criterion(&sc_qa);
165 TerminationCriterion sc_of, sc_inner;
166 sc_of.add_iteration_limit( 10 );
167 pass1.set_outer_termination_criterion(&sc_of);
168 sc_inner.add_iteration_limit( 6 );
169 //sc_inner.add_absolute_gradient_L2_norm( 0.01 );
170 pass1.set_inner_termination_criterion(&sc_inner);
171
172 //**********Set stopping criterion 5 iterates ****************
173 //StoppingCriterion sc5(StoppingCriterion::NUMBER_OF_PASSES,5);
174 //pass2->set_stopping_criterion(&sc5);
175 TerminationCriterion sc5;
176 sc5.add_iteration_limit( 5 );
177 pass2.set_inner_termination_criterion(&sc5);
178 //TerminationCriterion sc_inner;
179 //sc_inner.add_iteration_limit( 5 );
180 //pass2->set_inner_termination_criterion(&sc_inner);
181 //pass2->set_maximum_iteration(5);
182
183 queue1.set_master_quality_improver(&pass1, err); CPPUNIT_ASSERT(!err);
184 queue2.set_master_quality_improver(&pass2, err); CPPUNIT_ASSERT(!err);
185 //********************UNTANGLE*******************************
186 //Make sure no errors
187 CPPUNIT_ASSERT(!err);
188 // launches optimization on mesh_set1
189 MeshDomainAssoc mesh_and_domain = MeshDomainAssoc(&mesh, &msq_geom);
190 double orig_qa_val=stop_qa.loop_over_mesh(&mesh_and_domain, 0, err);
191 //Make sure no errors
192 CPPUNIT_ASSERT(!err);
193 queue1.run_instructions(&mesh_and_domain, err); CPPUNIT_ASSERT(!err);
194 //Make sure no errors
195 CPPUNIT_ASSERT(!err);
196 double fin_qa_val=stop_qa.loop_over_mesh(&mesh_and_domain, 0, err);
197 //Make sure no errors
198 CPPUNIT_ASSERT(!err);
199 //make sure 'quality' improved
200 CPPUNIT_ASSERT( (fin_qa_val-orig_qa_val) <= 0.0 );
201 //make sure sc_qa really was satisfied
202 CPPUNIT_ASSERT( fin_qa_val <= MSQ_MIN );
203
204 //********************SMOOTH*******************************
205 //Make sure no errors
206 CPPUNIT_ASSERT(!err);
207 // launches optimization on mesh_set1
208 orig_qa_val=qa.loop_over_mesh(&mesh_and_domain, 0, err);
209 //Make sure no errors
210 CPPUNIT_ASSERT(!err);
211 queue2.run_instructions(&mesh_and_domain, err); CPPUNIT_ASSERT(!err);
212 //Make sure no errors
213 CPPUNIT_ASSERT(!err);
214 fin_qa_val=qa.loop_over_mesh(&mesh_and_domain, 0, err);
215 //Make sure no errors
216 CPPUNIT_ASSERT(!err);
217 //make sure 'quality' improved
218 CPPUNIT_ASSERT( (fin_qa_val-orig_qa_val) <= 0.0 );
219 print_timing_diagnostics(cout);
220 }
221
test_plane_quad_tangled()222 void test_plane_quad_tangled()
223 {
224 MBMesquite::MeshImpl mesh;
225 MsqPrintError err(cout);
226 mesh.read_vtk(MESH_FILES_DIR "2D/vtk/quads/tangled/tangled_quad.vtk", err);
227 CPPUNIT_ASSERT(!err);
228
229 //create geometry: plane z=5, normal (0,0,1)
230 Vector3D pnt(0,0,5);
231 Vector3D s_norm(0,0,1);
232 MBMesquite::PlanarDomain msq_geom(s_norm, pnt);
233
234 // creates an intruction queue
235 InstructionQueue queue1, queue2;
236
237 //creates a mean ratio quality metric ...
238 ConditionNumberQualityMetric shape;
239 UntangleBetaQualityMetric untan(.1);
240
241 // ... and builds an objective function with it (untangle)
242 LInfTemplate untan_func(&untan);
243 LPtoPTemplate shape_func(&shape,2,err);
244 //Make sure no errors
245 CPPUNIT_ASSERT(!err);
246 // creates the cg optimization procedures
247 ConjugateGradient pass1( &untan_func, err );
248 ConjugateGradient pass2( &shape_func, err );
249 pass1.use_element_on_vertex_patch();
250 pass2.use_global_patch();
251 //Make sure no errors
252 CPPUNIT_ASSERT(!err);
253 QualityAssessor stop_qa=QualityAssessor( &untan );
254 QualityAssessor qa=QualityAssessor( &shape );
255 //turn off printing if print flag not set.
256 if(pF==0){
257 stop_qa.disable_printing_results();
258 qa.disable_printing_results();
259 }
260 //**********Set stopping criterion untangle ver small ********
261 //StoppingCriterion sc_qa(&stop_qa,-100,MSQ_MIN);
262 //pass1->set_stopping_criterion(&sc_qa);
263 TerminationCriterion sc_of;
264 sc_of.add_iteration_limit( 10 );
265 pass1.set_outer_termination_criterion(&sc_of);
266
267 //**********Set stopping criterion 5 iterates ****************
268 //StoppingCriterion sc5(StoppingCriterion::NUMBER_OF_PASSES,5);
269 //pass2->set_stopping_criterion(&sc5);
270 TerminationCriterion sc5;
271 sc5.add_iteration_limit( 5 );
272 pass2.set_inner_termination_criterion(&sc5);
273 //pass2->set_maximum_iteration(5);
274 //TerminationCriterion sc_inner;
275 //sc_inner.add_iteration_limit( 5 );
276 //pass2->set_inner_termination_criterion(&sc_inner);
277 queue1.set_master_quality_improver(&pass1, err); CPPUNIT_ASSERT(!err);
278 queue2.set_master_quality_improver(&pass2, err); CPPUNIT_ASSERT(!err);
279 //********************UNTANGLE*******************************
280 //Make sure no errors
281 CPPUNIT_ASSERT(!err);
282 // launches optimization on mesh_set1
283 MeshDomainAssoc mesh_and_domain = MeshDomainAssoc(&mesh, &msq_geom);
284 double orig_qa_val=stop_qa.loop_over_mesh(&mesh_and_domain, 0, err);
285 //Make sure no errors
286 CPPUNIT_ASSERT(!err);
287 queue1.run_instructions(&mesh_and_domain, err); CPPUNIT_ASSERT(!err);
288 //Make sure no errors
289 CPPUNIT_ASSERT(!err);
290 double fin_qa_val=stop_qa.loop_over_mesh(&mesh_and_domain, 0, err);
291 //Make sure no errors
292 CPPUNIT_ASSERT(!err);
293 //make sure 'quality' improved
294 CPPUNIT_ASSERT( (fin_qa_val-orig_qa_val) <= 0.0 );
295 //make sure sc_qa really was satisfied
296 CPPUNIT_ASSERT( fin_qa_val <= MSQ_MIN );
297
298 //********************SMOOTH*******************************
299 //Make sure no errors
300 CPPUNIT_ASSERT(!err);
301 // launches optimization on mesh_set1
302 orig_qa_val=qa.loop_over_mesh(&mesh_and_domain, 0, err);
303 //Make sure no errors
304 CPPUNIT_ASSERT(!err);
305 queue2.run_instructions(&mesh_and_domain, err); CPPUNIT_ASSERT(!err);
306 //Make sure no errors
307 CPPUNIT_ASSERT(!err);
308 fin_qa_val=qa.loop_over_mesh(&mesh_and_domain, 0, err);
309 //Make sure no errors
310 CPPUNIT_ASSERT(!err);
311 //make sure 'quality' improved
312 CPPUNIT_ASSERT( (fin_qa_val-orig_qa_val) <= 0.0 );
313 print_timing_diagnostics(cout);
314 }
315
test_plane_tri_xz()316 void test_plane_tri_xz()
317 {
318 MsqPrintError err(cout);
319 MBMesquite::MeshImpl mesh;
320 mesh.read_vtk(MESH_FILES_DIR "2D/vtk/tris/untangled/tri_5_xz.vtk", err);
321 CPPUNIT_ASSERT(!err);
322
323 //create geometry: plane y=5, normal (0,1,0)
324 Vector3D pnt(0,-5,0);
325 Vector3D s_norm(0,-1,0);
326 MBMesquite::PlanarDomain msq_geom(s_norm, pnt);
327
328 // creates an intruction queue
329 InstructionQueue queue1;
330
331 //creates a asm quality metric ...
332 ConditionNumberQualityMetric smooth;
333
334 // ... and builds an objective function with it (untangle)
335 LPtoPTemplate smooth_func(&smooth,1,err);
336 //Make sure no errors
337 CPPUNIT_ASSERT(!err);
338 // creates the cg optimization procedures
339 ConjugateGradient pass1( &smooth_func, err );
340 //pass1->set_patch_type(PatchData::ELEMENTS_ON_VERTEX_PATCH, err,1 ,1);
341 pass1.use_global_patch();
342 pass1.set_debugging_level(1);
343 //Make sure no errors
344 CPPUNIT_ASSERT(!err);
345 QualityAssessor qa=QualityAssessor( &smooth );
346
347 //**********Set stopping criterion 5 iterates ****************
348 TerminationCriterion sc5;
349 sc5.add_iteration_limit( 5 );
350 pass1.set_inner_termination_criterion(&sc5);
351 //StoppingCriterion sc5(StoppingCriterion::NUMBER_OF_PASSES,5);
352 //pass1->set_stopping_criterion(&sc5);
353 TerminationCriterion sc_inner;
354 sc_inner.add_iteration_limit( 5 );
355 pass1.set_inner_termination_criterion(&sc_inner);
356 //pass1->set_maximum_iteration(5);
357
358 queue1.set_master_quality_improver(&pass1, err); CPPUNIT_ASSERT(!err);
359 //********************UNTANGLE*******************************
360 //Make sure no errors
361 CPPUNIT_ASSERT(!err);
362 // launches optimization on mesh_set1
363 MeshDomainAssoc mesh_and_domain = MeshDomainAssoc(&mesh, &msq_geom);
364 double orig_qa_val=qa.loop_over_mesh(&mesh_and_domain, 0, err);
365 //Make sure no errors
366 CPPUNIT_ASSERT(!err);
367 queue1.run_instructions(&mesh_and_domain, err); CPPUNIT_ASSERT(!err);
368 //Make sure no errors
369 CPPUNIT_ASSERT(!err);
370 double fin_qa_val=qa.loop_over_mesh(&mesh_and_domain, 0, err);
371 //Make sure no errors
372 CPPUNIT_ASSERT(!err);
373 //make sure 'quality' improved
374 CPPUNIT_ASSERT( (fin_qa_val-orig_qa_val) <= 0.0 );
375 print_timing_diagnostics(cout);
376 }
377
378 void test_fit_plane();
379 };
380
381
382 CPPUNIT_TEST_SUITE_NAMED_REGISTRATION(PlanarGeometryTest, "PlanarGeometryTest");
383 CPPUNIT_TEST_SUITE_NAMED_REGISTRATION(PlanarGeometryTest, "Regression");
384
test_fit_plane()385 void PlanarGeometryTest::test_fit_plane()
386 {
387 MsqPrintError err(std::cerr);
388 PlanarDomain plane(PlanarDomain::XY,-1);
389 const double epsilon = 1e-8;
390
391 MeshImpl mesh1;
392 mesh1.read_vtk(MESH_FILES_DIR "2D/vtk/tris/untangled/bad_circle_tri.vtk", err);
393 ASSERT_NO_ERROR(err);
394 plane.fit_vertices( &mesh1, err, epsilon );
395 ASSERT_NO_ERROR(err);
396 CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D(0,0,1), plane.get_normal(), epsilon );
397 CPPUNIT_ASSERT_DOUBLES_EQUAL( 5, plane.get_coeff(), epsilon );
398
399 MeshImpl mesh2;
400 mesh2.read_vtk(MESH_FILES_DIR "2D/vtk/tris/untangled/equil_tri.vtk", err);
401 ASSERT_NO_ERROR(err);
402 plane.fit_vertices( &mesh2, err, epsilon );
403 ASSERT_NO_ERROR(err);
404 CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D(0,0,1), plane.get_normal(), epsilon );
405 CPPUNIT_ASSERT_DOUBLES_EQUAL( 0, plane.get_coeff(), epsilon );
406
407 MeshImpl mesh3;
408 mesh3.read_vtk(MESH_FILES_DIR "2D/vtk/quads/untangled/quads_4by2.vtk", err);
409 ASSERT_NO_ERROR(err);
410 plane.fit_vertices( &mesh3, err, epsilon );
411 ASSERT_NO_ERROR(err);
412 CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D(0,0,1), plane.get_normal(), epsilon );
413 CPPUNIT_ASSERT_DOUBLES_EQUAL( -2, plane.get_coeff(), epsilon );
414
415 MeshImpl mesh4;
416 mesh4.read_vtk(MESH_FILES_DIR "2D/vtk/tris/untangled/tri_5_xz.vtk", err);
417 ASSERT_NO_ERROR(err);
418 plane.fit_vertices( &mesh4, err, epsilon );
419 ASSERT_NO_ERROR(err);
420 CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D(0,-1,0), plane.get_normal(), epsilon );
421 CPPUNIT_ASSERT_DOUBLES_EQUAL( -5, plane.get_coeff(), epsilon );
422 }
423