MOAB: Mesh Oriented datABase  (version 5.2.1)
PlanarGeometryTest.cpp
Go to the documentation of this file.
00001 /* *****************************************************************
00002     MESQUITE -- The Mesh Quality Improvement Toolkit
00003 
00004     Copyright 2004 Sandia Corporation and Argonne National
00005     Laboratory.  Under the terms of Contract DE-AC04-94AL85000
00006     with Sandia Corporation, the U.S. Government retains certain
00007     rights in this software.
00008 
00009     This library is free software; you can redistribute it and/or
00010     modify it under the terms of the GNU Lesser General Public
00011     License as published by the Free Software Foundation; either
00012     version 2.1 of the License, or (at your option) any later version.
00013 
00014     This library is distributed in the hope that it will be useful,
00015     but WITHOUT ANY WARRANTY; without even the implied warranty of
00016     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00017     Lesser General Public License for more details.
00018 
00019     You should have received a copy of the GNU Lesser General Public License
00020     (lgpl.txt) along with this library; if not, write to the Free Software
00021     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00022 
00023     diachin2@llnl.gov, djmelan@sandia.gov, mbrewer@sandia.gov,
00024     pknupp@sandia.gov, tleurent@mcs.anl.gov, tmunson@mcs.anl.gov
00025 
00026   ***************************************************************** */
00027 // -*- Mode : c++; tab-width: 2; c-tab-always-indent: t; indent-tabs-mode: nil; c-basic-offset: 2
00028 // -*-
00029 //
00030 //   SUMMARY:
00031 //     USAGE:
00032 //
00033 //    AUTHOR: Michael Brewer
00034 //       ORG: Sandia National Labs
00035 //    E-MAIL: mbrewer@sandia.gov
00036 //
00037 // ORIG-DATE: Jan. 29, 2003
00038 //  LAST-MOD: 25-Feb-04 at 10:49:32 by Thomas Leurent
00039 //
00040 // DESCRIPTION:
00041 // ============
00042 /*! \file PlanarGeometryTest.cpp
00043 
00044 Regression testing using the planar geometry capabilities in
00045 SimplifiedGeometryEngine.
00046  */
00047 // DESCRIP-END.
00048 //
00049 
00050 #include "meshfiles.h"
00051 
00052 #include "PatchDataInstances.hpp"
00053 #include "cppunit/extensions/HelperMacros.h"
00054 #include <math.h>
00055 
00056 #include "Mesquite.hpp"
00057 #include "MsqError.hpp"
00058 #include "InstructionQueue.hpp"
00059 #include "PatchData.hpp"
00060 //#include "StoppingCriterion.hpp"
00061 #include "QualityAssessor.hpp"
00062 
00063 #include "ConditionNumberQualityMetric.hpp"
00064 #include "LPtoPTemplate.hpp"
00065 #include "EdgeLengthQualityMetric.hpp"
00066 #include "LaplacianSmoother.hpp"
00067 #include "LInfTemplate.hpp"
00068 #include "SteepestDescent.hpp"
00069 #include "ConjugateGradient.hpp"
00070 #include "AspectRatioGammaQualityMetric.hpp"
00071 #include "UntangleBetaQualityMetric.hpp"
00072 #include "MultiplyQualityMetric.hpp"
00073 #include "PlanarDomain.hpp"
00074 
00075 #include "MsqFPE.hpp"
00076 
00077 #include "UnitUtil.hpp"
00078 
00079 #include "MeshImpl.hpp"
00080 #include <iostream>
00081 using std::cout;
00082 using std::endl;
00083 
00084 using namespace MBMesquite;
00085 
00086 class PlanarGeometryTest : public CppUnit::TestFixture
00087 {
00088   private:
00089     CPPUNIT_TEST_SUITE( PlanarGeometryTest );
00090     // run steepest descent on the tangled_tri.vtk mesh
00091     CPPUNIT_TEST( test_plane_tri_tangled );
00092     // run cg on tangled_quad.vtk mesh
00093     CPPUNIT_TEST( test_plane_quad_tangled );
00094     // run cg with asm metric on tri mesh in y=-5 plane
00095     CPPUNIT_TEST( test_plane_tri_xz );
00096     // test fit plane
00097     CPPUNIT_TEST( test_fit_plane );
00098     CPPUNIT_TEST_SUITE_END();
00099 
00100   private:
00101     double qualTol;  // double used for double comparisons
00102     int pF;          // PRINT_FLAG
00103   public:
00104     void setUp()
00105     {
00106         // pF=1;//PRINT_FLAG IS ON
00107         pF = 0;  // PRINT_FLAG IS OFF
00108                  // tolerance double
00109         qualTol = MSQ_MIN;
00110     }
00111 
00112     void tearDown() {}
00113 
00114   public:
00115     PlanarGeometryTest() {}
00116 
00117     void test_plane_tri_tangled()
00118     {
00119         MBMesquite::MsqPrintError err( cout );
00120         MBMesquite::MeshImpl mesh;
00121 
00122         // This test doesn't use InstructionQueue, so
00123         // we need to set up trapping of floating-point
00124         // exceptions ourself.
00125         MsqFPE fpe_trap( true );
00126 
00127         mesh.read_vtk( MESH_FILES_DIR "2D/vtk/tris/tangled/tangled_tri.vtk", err );
00128         CPPUNIT_ASSERT( !err );
00129 
00130         // create geometry: plane z=5, normal (0,0,1)
00131         Vector3D pnt( 0, 0, 5 );
00132         Vector3D s_norm( 0, 0, 1 );
00133         MBMesquite::PlanarDomain msq_geom( s_norm, pnt );
00134 
00135         // creates an intruction queue
00136         InstructionQueue queue1, queue2;
00137 
00138         // creates a mean ratio quality metric ...
00139         ConditionNumberQualityMetric shape;
00140         UntangleBetaQualityMetric untan( .1 );
00141 
00142         // ... and builds an objective function with it (untangle)
00143         LInfTemplate untan_func( &untan );
00144         LPtoPTemplate shape_func( &shape, 2, err );
00145         // Make sure no errors
00146         CPPUNIT_ASSERT( !err );
00147         // creates the steepest descent optimization procedures
00148         SteepestDescent pass1( &untan_func );
00149         SteepestDescent pass2( &shape_func );
00150         pass1.use_element_on_vertex_patch();
00151         pass2.use_global_patch();
00152         // Make sure no errors
00153         CPPUNIT_ASSERT( !err );
00154         QualityAssessor stop_qa = QualityAssessor( &untan );
00155         QualityAssessor qa      = QualityAssessor( &shape );
00156         if( pF == 0 )
00157         {
00158             stop_qa.disable_printing_results();
00159             qa.disable_printing_results();
00160         }
00161         //**********Set stopping criterion  untangle ver small ********
00162         // StoppingCriterion sc_qa(&stop_qa,-100,MSQ_MIN);
00163         // pass1->set_stopping_criterion(&sc_qa);
00164         TerminationCriterion sc_of, sc_inner;
00165         sc_of.add_iteration_limit( 10 );
00166         pass1.set_outer_termination_criterion( &sc_of );
00167         sc_inner.add_iteration_limit( 6 );
00168         // sc_inner.add_absolute_gradient_L2_norm( 0.01 );
00169         pass1.set_inner_termination_criterion( &sc_inner );
00170 
00171         //**********Set stopping criterion  5 iterates ****************
00172         // StoppingCriterion sc5(StoppingCriterion::NUMBER_OF_PASSES,5);
00173         // pass2->set_stopping_criterion(&sc5);
00174         TerminationCriterion sc5;
00175         sc5.add_iteration_limit( 5 );
00176         pass2.set_inner_termination_criterion( &sc5 );
00177         // TerminationCriterion sc_inner;
00178         // sc_inner.add_iteration_limit( 5 );
00179         // pass2->set_inner_termination_criterion(&sc_inner);
00180         // pass2->set_maximum_iteration(5);
00181 
00182         queue1.set_master_quality_improver( &pass1, err );
00183         CPPUNIT_ASSERT( !err );
00184         queue2.set_master_quality_improver( &pass2, err );
00185         CPPUNIT_ASSERT( !err );
00186         //********************UNTANGLE*******************************
00187         // Make sure no errors
00188         CPPUNIT_ASSERT( !err );
00189         // launches optimization on mesh_set1
00190         MeshDomainAssoc mesh_and_domain = MeshDomainAssoc( &mesh, &msq_geom );
00191         double orig_qa_val              = stop_qa.loop_over_mesh( &mesh_and_domain, 0, err );
00192         // Make sure no errors
00193         CPPUNIT_ASSERT( !err );
00194         queue1.run_instructions( &mesh_and_domain, err );
00195         CPPUNIT_ASSERT( !err );
00196         // Make sure no errors
00197         CPPUNIT_ASSERT( !err );
00198         double fin_qa_val = stop_qa.loop_over_mesh( &mesh_and_domain, 0, err );
00199         // Make sure no errors
00200         CPPUNIT_ASSERT( !err );
00201         // make sure 'quality' improved
00202         CPPUNIT_ASSERT( ( fin_qa_val - orig_qa_val ) <= 0.0 );
00203         // make sure sc_qa really was satisfied
00204         CPPUNIT_ASSERT( fin_qa_val <= MSQ_MIN );
00205 
00206         //********************SMOOTH*******************************
00207         // Make sure no errors
00208         CPPUNIT_ASSERT( !err );
00209         // launches optimization on mesh_set1
00210         orig_qa_val = qa.loop_over_mesh( &mesh_and_domain, 0, err );
00211         // Make sure no errors
00212         CPPUNIT_ASSERT( !err );
00213         queue2.run_instructions( &mesh_and_domain, err );
00214         CPPUNIT_ASSERT( !err );
00215         // Make sure no errors
00216         CPPUNIT_ASSERT( !err );
00217         fin_qa_val = qa.loop_over_mesh( &mesh_and_domain, 0, err );
00218         // Make sure no errors
00219         CPPUNIT_ASSERT( !err );
00220         // make sure 'quality' improved
00221         CPPUNIT_ASSERT( ( fin_qa_val - orig_qa_val ) <= 0.0 );
00222         print_timing_diagnostics( cout );
00223     }
00224 
00225     void test_plane_quad_tangled()
00226     {
00227         MBMesquite::MeshImpl mesh;
00228         MsqPrintError err( cout );
00229         mesh.read_vtk( MESH_FILES_DIR "2D/vtk/quads/tangled/tangled_quad.vtk", err );
00230         CPPUNIT_ASSERT( !err );
00231 
00232         // create geometry: plane z=5, normal (0,0,1)
00233         Vector3D pnt( 0, 0, 5 );
00234         Vector3D s_norm( 0, 0, 1 );
00235         MBMesquite::PlanarDomain msq_geom( s_norm, pnt );
00236 
00237         // creates an intruction queue
00238         InstructionQueue queue1, queue2;
00239 
00240         // creates a mean ratio quality metric ...
00241         ConditionNumberQualityMetric shape;
00242         UntangleBetaQualityMetric untan( .1 );
00243 
00244         // ... and builds an objective function with it (untangle)
00245         LInfTemplate untan_func( &untan );
00246         LPtoPTemplate shape_func( &shape, 2, err );
00247         // Make sure no errors
00248         CPPUNIT_ASSERT( !err );
00249         // creates the cg optimization procedures
00250         ConjugateGradient pass1( &untan_func, err );
00251         ConjugateGradient pass2( &shape_func, err );
00252         pass1.use_element_on_vertex_patch();
00253         pass2.use_global_patch();
00254         // Make sure no errors
00255         CPPUNIT_ASSERT( !err );
00256         QualityAssessor stop_qa = QualityAssessor( &untan );
00257         QualityAssessor qa      = QualityAssessor( &shape );
00258         // turn off printing if print flag not set.
00259         if( pF == 0 )
00260         {
00261             stop_qa.disable_printing_results();
00262             qa.disable_printing_results();
00263         }
00264         //**********Set stopping criterion  untangle ver small ********
00265         // StoppingCriterion sc_qa(&stop_qa,-100,MSQ_MIN);
00266         // pass1->set_stopping_criterion(&sc_qa);
00267         TerminationCriterion sc_of;
00268         sc_of.add_iteration_limit( 10 );
00269         pass1.set_outer_termination_criterion( &sc_of );
00270 
00271         //**********Set stopping criterion  5 iterates ****************
00272         // StoppingCriterion sc5(StoppingCriterion::NUMBER_OF_PASSES,5);
00273         // pass2->set_stopping_criterion(&sc5);
00274         TerminationCriterion sc5;
00275         sc5.add_iteration_limit( 5 );
00276         pass2.set_inner_termination_criterion( &sc5 );
00277         // pass2->set_maximum_iteration(5);
00278         // TerminationCriterion sc_inner;
00279         // sc_inner.add_iteration_limit( 5 );
00280         // pass2->set_inner_termination_criterion(&sc_inner);
00281         queue1.set_master_quality_improver( &pass1, err );
00282         CPPUNIT_ASSERT( !err );
00283         queue2.set_master_quality_improver( &pass2, err );
00284         CPPUNIT_ASSERT( !err );
00285         //********************UNTANGLE*******************************
00286         // Make sure no errors
00287         CPPUNIT_ASSERT( !err );
00288         // launches optimization on mesh_set1
00289         MeshDomainAssoc mesh_and_domain = MeshDomainAssoc( &mesh, &msq_geom );
00290         double orig_qa_val              = stop_qa.loop_over_mesh( &mesh_and_domain, 0, err );
00291         // Make sure no errors
00292         CPPUNIT_ASSERT( !err );
00293         queue1.run_instructions( &mesh_and_domain, err );
00294         CPPUNIT_ASSERT( !err );
00295         // Make sure no errors
00296         CPPUNIT_ASSERT( !err );
00297         double fin_qa_val = stop_qa.loop_over_mesh( &mesh_and_domain, 0, err );
00298         // Make sure no errors
00299         CPPUNIT_ASSERT( !err );
00300         // make sure 'quality' improved
00301         CPPUNIT_ASSERT( ( fin_qa_val - orig_qa_val ) <= 0.0 );
00302         // make sure sc_qa really was satisfied
00303         CPPUNIT_ASSERT( fin_qa_val <= MSQ_MIN );
00304 
00305         //********************SMOOTH*******************************
00306         // Make sure no errors
00307         CPPUNIT_ASSERT( !err );
00308         // launches optimization on mesh_set1
00309         orig_qa_val = qa.loop_over_mesh( &mesh_and_domain, 0, err );
00310         // Make sure no errors
00311         CPPUNIT_ASSERT( !err );
00312         queue2.run_instructions( &mesh_and_domain, err );
00313         CPPUNIT_ASSERT( !err );
00314         // Make sure no errors
00315         CPPUNIT_ASSERT( !err );
00316         fin_qa_val = qa.loop_over_mesh( &mesh_and_domain, 0, err );
00317         // Make sure no errors
00318         CPPUNIT_ASSERT( !err );
00319         // make sure 'quality' improved
00320         CPPUNIT_ASSERT( ( fin_qa_val - orig_qa_val ) <= 0.0 );
00321         print_timing_diagnostics( cout );
00322     }
00323 
00324     void test_plane_tri_xz()
00325     {
00326         MsqPrintError err( cout );
00327         MBMesquite::MeshImpl mesh;
00328         mesh.read_vtk( MESH_FILES_DIR "2D/vtk/tris/untangled/tri_5_xz.vtk", err );
00329         CPPUNIT_ASSERT( !err );
00330 
00331         // create geometry: plane y=5, normal (0,1,0)
00332         Vector3D pnt( 0, -5, 0 );
00333         Vector3D s_norm( 0, -1, 0 );
00334         MBMesquite::PlanarDomain msq_geom( s_norm, pnt );
00335 
00336         // creates an intruction queue
00337         InstructionQueue queue1;
00338 
00339         // creates a asm quality metric ...
00340         ConditionNumberQualityMetric smooth;
00341 
00342         // ... and builds an objective function with it (untangle)
00343         LPtoPTemplate smooth_func( &smooth, 1, err );
00344         // Make sure no errors
00345         CPPUNIT_ASSERT( !err );
00346         // creates the cg optimization procedures
00347         ConjugateGradient pass1( &smooth_func, err );
00348         // pass1->set_patch_type(PatchData::ELEMENTS_ON_VERTEX_PATCH, err,1 ,1);
00349         pass1.use_global_patch();
00350         pass1.set_debugging_level( 1 );
00351         // Make sure no errors
00352         CPPUNIT_ASSERT( !err );
00353         QualityAssessor qa = QualityAssessor( &smooth );
00354 
00355         //**********Set stopping criterion  5 iterates ****************
00356         TerminationCriterion sc5;
00357         sc5.add_iteration_limit( 5 );
00358         pass1.set_inner_termination_criterion( &sc5 );
00359         // StoppingCriterion sc5(StoppingCriterion::NUMBER_OF_PASSES,5);
00360         // pass1->set_stopping_criterion(&sc5);
00361         TerminationCriterion sc_inner;
00362         sc_inner.add_iteration_limit( 5 );
00363         pass1.set_inner_termination_criterion( &sc_inner );
00364         // pass1->set_maximum_iteration(5);
00365 
00366         queue1.set_master_quality_improver( &pass1, err );
00367         CPPUNIT_ASSERT( !err );
00368         //********************UNTANGLE*******************************
00369         // Make sure no errors
00370         CPPUNIT_ASSERT( !err );
00371         // launches optimization on mesh_set1
00372         MeshDomainAssoc mesh_and_domain = MeshDomainAssoc( &mesh, &msq_geom );
00373         double orig_qa_val              = qa.loop_over_mesh( &mesh_and_domain, 0, err );
00374         // Make sure no errors
00375         CPPUNIT_ASSERT( !err );
00376         queue1.run_instructions( &mesh_and_domain, err );
00377         CPPUNIT_ASSERT( !err );
00378         // Make sure no errors
00379         CPPUNIT_ASSERT( !err );
00380         double fin_qa_val = qa.loop_over_mesh( &mesh_and_domain, 0, err );
00381         // Make sure no errors
00382         CPPUNIT_ASSERT( !err );
00383         // make sure 'quality' improved
00384         CPPUNIT_ASSERT( ( fin_qa_val - orig_qa_val ) <= 0.0 );
00385         print_timing_diagnostics( cout );
00386     }
00387 
00388     void test_fit_plane();
00389 };
00390 
00391 CPPUNIT_TEST_SUITE_NAMED_REGISTRATION( PlanarGeometryTest, "PlanarGeometryTest" );
00392 CPPUNIT_TEST_SUITE_NAMED_REGISTRATION( PlanarGeometryTest, "Regression" );
00393 
00394 void PlanarGeometryTest::test_fit_plane()
00395 {
00396     MsqPrintError err( std::cerr );
00397     PlanarDomain plane( PlanarDomain::XY, -1 );
00398     const double epsilon = 1e-8;
00399 
00400     MeshImpl mesh1;
00401     mesh1.read_vtk( MESH_FILES_DIR "2D/vtk/tris/untangled/bad_circle_tri.vtk", err );
00402     ASSERT_NO_ERROR( err );
00403     plane.fit_vertices( &mesh1, err, epsilon );
00404     ASSERT_NO_ERROR( err );
00405     CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( 0, 0, 1 ), plane.get_normal(), epsilon );
00406     CPPUNIT_ASSERT_DOUBLES_EQUAL( 5, plane.get_coeff(), epsilon );
00407 
00408     MeshImpl mesh2;
00409     mesh2.read_vtk( MESH_FILES_DIR "2D/vtk/tris/untangled/equil_tri.vtk", err );
00410     ASSERT_NO_ERROR( err );
00411     plane.fit_vertices( &mesh2, err, epsilon );
00412     ASSERT_NO_ERROR( err );
00413     CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( 0, 0, 1 ), plane.get_normal(), epsilon );
00414     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0, plane.get_coeff(), epsilon );
00415 
00416     MeshImpl mesh3;
00417     mesh3.read_vtk( MESH_FILES_DIR "2D/vtk/quads/untangled/quads_4by2.vtk", err );
00418     ASSERT_NO_ERROR( err );
00419     plane.fit_vertices( &mesh3, err, epsilon );
00420     ASSERT_NO_ERROR( err );
00421     CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( 0, 0, 1 ), plane.get_normal(), epsilon );
00422     CPPUNIT_ASSERT_DOUBLES_EQUAL( -2, plane.get_coeff(), epsilon );
00423 
00424     MeshImpl mesh4;
00425     mesh4.read_vtk( MESH_FILES_DIR "2D/vtk/tris/untangled/tri_5_xz.vtk", err );
00426     ASSERT_NO_ERROR( err );
00427     plane.fit_vertices( &mesh4, err, epsilon );
00428     ASSERT_NO_ERROR( err );
00429     CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( 0, -1, 0 ), plane.get_normal(), epsilon );
00430     CPPUNIT_ASSERT_DOUBLES_EQUAL( -5, plane.get_coeff(), epsilon );
00431 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines