MOAB: Mesh Oriented datABase  (version 5.3.0)
main.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: 3; c-tab-always-indent: t; indent-tabs-mode: nil; c-basic-offset: 3
00028 // -*-
00029 //
00030 //   SUMMARY:
00031 //     USAGE:
00032 //
00033 // ORIG-DATE: 19-Feb-02 at 10:57:52
00034 //  LAST-MOD: 23-Jul-03 at 18:09:51 by Thomas Leurent
00035 //
00036 //
00037 // DESCRIPTION:
00038 // ============
00039 /*! \file main.cpp
00040 
00041 describe main.cpp here
00042 
00043  */
00044 // DESCRIP-END.
00045 //
00046 
00047 #include <iostream>
00048 using std::cout;
00049 using std::endl;
00050 #include <cstdlib>
00051 
00052 #include "TestUtil.hpp"
00053 #include "Mesquite.hpp"
00054 #include "MeshImpl.hpp"
00055 #include "MsqError.hpp"
00056 #include "Vector3D.hpp"
00057 #include "InstructionQueue.hpp"
00058 #include "PatchData.hpp"
00059 #include "TerminationCriterion.hpp"
00060 #include "QualityAssessor.hpp"
00061 //#include "QuadLagrangeShape.hpp"
00062 
00063 // algorithms
00064 #include "IdealShapeTarget.hpp"
00065 #include "TQualityMetric.hpp"
00066 #include "TInverseMeanRatio.hpp"
00067 #include "ConditionNumberQualityMetric.hpp"
00068 #include "LPtoPTemplate.hpp"
00069 #include "LInfTemplate.hpp"
00070 #include "SteepestDescent.hpp"
00071 #include "ConjugateGradient.hpp"
00072 
00073 #include "PlanarDomain.hpp"
00074 using namespace MBMesquite;
00075 
00076 /* This is the input mesh topology
00077      (0)------(16)-----(1)------(17)-----(2)------(18)-----(3)
00078       |                 |                 |                 |
00079       |                 |                 |                 |
00080       |                 |                 |                 |
00081       |                 |                 |                 |
00082      (19)      0       (20)      1       (21)      2       (22)
00083       |                 |                 |                 |
00084       |                 |                 |                 |
00085       |                 |                 |                 |
00086       |                 |                 |                 |
00087      (4)------(23)-----(5)------(24)-----(6)------(25)-----(7)
00088       |                 |                 |                 |
00089       |                 |                 |                 |
00090       |                 |                 |                 |
00091       |                 |                 |                 |
00092      (26)      3       (27)      4       (28)      5       (29)
00093       |                 |                 |                 |
00094       |                 |                 |                 |
00095       |                 |                 |                 |
00096       |                 |                 |                 |
00097      (8)------(30)-----(9)------(31)-----(10)-----(32)-----(11)
00098       |                 |                 |                 |
00099       |                 |                 |                 |
00100       |                 |                 |                 |
00101       |                 |                 |                 |
00102      (33)      6       (34)      7       (35)      8       (36)
00103       |                 |                 |                 |
00104       |                 |                 |                 |
00105       |                 |                 |                 |
00106       |                 |                 |                 |
00107      (12)-----(37)-----(13)-----(38)-----(14)-----(39)-----(15)
00108 */
00109 
00110 std::string LINEAR_INPUT_FILE_NAME       = std::string( STRINGIFY( SRCDIR ) ) + "/linear_input.vtk";
00111 std::string QUADRATIC_INPUT_FILE_NAME    = std::string( STRINGIFY( SRCDIR ) ) + "/quadratic_input.vtk";
00112 std::string EXPECTED_LINAR_FILE_NAME     = std::string( STRINGIFY( SRCDIR ) ) + "/expected_linear_output.vtk";
00113 std::string EXPECTED_QUADRATIC_FILE_NAME = std::string( STRINGIFY( SRCDIR ) ) + "/expected_quadratic_output.vtk";
00114 std::string HOUR_INPUT_FILE_NAME         = std::string( STRINGIFY( SRCDIR ) ) + "/hour-quad8.vtk";
00115 std::string OUTPUT_FILE_NAME             = "smoothed_qudratic_mesh.vtk";
00116 const unsigned NUM_CORNER_VERTICES       = 16;
00117 const unsigned NUM_MID_NODES             = 24;
00118 const double SPATIAL_COMPARE_TOLERANCE   = 4e-6;
00119 
00120 void compare_nodes( size_t start_index, size_t end_index, Mesh* mesh1, Mesh* mesh2, MsqError& err )
00121 {
00122     size_t i, num_verts = end_index - start_index;
00123     std::vector< MsqVertex > verts1( num_verts ), verts2( num_verts );
00124     std::vector< Mesh::VertexHandle > handles1( num_verts ), handles2( num_verts );
00125 
00126     /* VertexIterator skips higher-order nodes.
00127        For now, just assume index == handle
00128 
00129       std::vector<Mesh::VertexHandle>::iterator handle_iter1, handle_iter2;
00130 
00131         // Skip start_index vertices
00132       VertexIterator* iter1 = mesh1->vertex_iterator( err ); MSQ_ERRRTN(err);
00133       VertexIterator* iter2 = mesh2->vertex_iterator( err ); MSQ_ERRRTN(err);
00134       for (i = 0; i < start_index; ++i)
00135       {
00136         if (iter1->is_at_end())
00137         {
00138           MSQ_SETERR(err)("start index out of range for first mesh set", MsqError::INVALID_ARG);
00139           return;
00140         }
00141         if (iter2->is_at_end())
00142         {
00143           MSQ_SETERR(err)("start index out of range for second mesh set", MsqError::INVALID_ARG);
00144           return;
00145         }
00146         iter1->operator++();
00147         iter2->operator++();
00148       }
00149 
00150         // Get handles for vertices
00151       handle_iter1 = handles1.begin();
00152       handle_iter2 = handles2.begin();
00153       for (i = start_index; i < end_index; ++i)
00154       {
00155         if (iter1->is_at_end())
00156         {
00157           MSQ_SETERR(err)("end index out of range for first mesh set", MsqError::INVALID_ARG);
00158           return;
00159         }
00160         *handle_iter1 = iter1->operator*();
00161         iter1->operator++();
00162         ++handle_iter1;
00163 
00164         if (iter2->is_at_end())
00165         {
00166           MSQ_SETERR(err)("end index out of range for second mesh set", MsqError::INVALID_ARG);
00167           return;
00168         }
00169         *handle_iter2 = iter2->operator*();
00170         iter2->operator++();
00171         ++handle_iter2;
00172       }
00173     */
00174     for( i = start_index; i < end_index; ++i )
00175         handles1[i - start_index] = handles2[i - start_index] = (void*)i;
00176 
00177     // Get coordinates from handles
00178     mesh1->vertices_get_coordinates( arrptr( handles1 ), arrptr( verts1 ), num_verts, err );MSQ_ERRRTN( err );
00179     mesh2->vertices_get_coordinates( arrptr( handles2 ), arrptr( verts2 ), num_verts, err );MSQ_ERRRTN( err );
00180 
00181     // Compare coordinates
00182     for( i = 0; i < num_verts; ++i )
00183     {
00184         const double diff = ( verts1[i] - verts2[i] ).length();
00185         if( diff > SPATIAL_COMPARE_TOLERANCE )
00186         {
00187             MSQ_SETERR( err )
00188             ( MsqError::INTERNAL_ERROR, "%u%s vertices differ. (%f,%f,%f) vs (%f,%f,%f)", (unsigned)( 1 + i ),
00189               i % 10 == 0   ? "st"
00190               : i % 10 == 1 ? "nd"
00191               : i % 10 == 2 ? "rd"
00192                             : "th",
00193               verts1[i][0], verts1[i][1], verts1[i][2], verts2[i][0], verts2[i][1], verts2[i][2] );
00194             return;
00195         }
00196     }
00197 }
00198 
00199 // code copied from testSuite/algorithm_test/main.cpp
00200 InstructionQueue* create_instruction_queue( MsqError& err )
00201 {
00202     // creates an intruction queue
00203     InstructionQueue* queue1 = new InstructionQueue;
00204 
00205     // creates a mean ratio quality metric ...
00206     // IdealWeightInverseMeanRatio* mean = new IdealWeightInverseMeanRatio(err); MSQ_ERRZERO(err);
00207     TargetCalculator* tc = new IdealShapeTarget;
00208     TQualityMetric* mean = new TQualityMetric( tc, 0, new TInverseMeanRatio );
00209 
00210     LPtoPTemplate* obj_func = new LPtoPTemplate( mean, 1, err );
00211     MSQ_ERRZERO( err );
00212 
00213     // creates the optimization procedures
00214     //   ConjugateGradient* pass1 = new ConjugateGradient( obj_func, err );
00215     SteepestDescent* pass1 = new SteepestDescent( obj_func );
00216 
00217     // perform optimization globally
00218     pass1->use_global_patch();
00219 
00220     // QualityAssessor* mean_qa = new QualityAssessor(mean);
00221 
00222     //**************Set termination criterion****************
00223 
00224     // perform 1 pass of the outer loop (this line isn't essential as it is
00225     // the default behavior).
00226     TerminationCriterion* tc_outer = new TerminationCriterion;
00227     tc_outer->add_iteration_limit( 1 );
00228     pass1->set_outer_termination_criterion( tc_outer );
00229 
00230     // perform the inner loop until a certain objective function value is
00231     // reached.  The exact value needs to be determined (about 18095).
00232     // As a safety, also stop if the time exceeds 10 minutes (600 seconds).
00233     TerminationCriterion* tc_inner = new TerminationCriterion;
00234     tc_inner->add_absolute_vertex_movement( 1e-6 );
00235     pass1->set_inner_termination_criterion( tc_inner );
00236 
00237     // adds 1 pass of pass1 to mesh_set1
00238     // queue1->add_quality_assessor(mean_qa,err); MSQ_ERRZERO(err);
00239     queue1->set_master_quality_improver( pass1, err );
00240     MSQ_ERRZERO( err );
00241     // queue1->add_quality_assessor(mean_qa,err); MSQ_ERRZERO(err);
00242 
00243     return queue1;
00244 }
00245 
00246 int do_test( bool slave )
00247 {
00248     MsqPrintError err( cout );
00249     //  QuadLagrangeShape quad9;
00250 
00251     // Create geometry
00252     Vector3D z( 0, 0, 1 ), o( 0, 0, 0 );
00253     PlanarDomain geom( z, o );
00254 
00255     // Read in linear input mesh
00256     cout << "Reading " << LINEAR_INPUT_FILE_NAME << endl;
00257     MeshImpl* linear_in = new MeshImpl;
00258     linear_in->read_vtk( LINEAR_INPUT_FILE_NAME.c_str(), err );
00259     if( MSQ_CHKERR( err ) ) return 1;
00260 
00261     // Read in expected linear results
00262     cout << "Reading " << EXPECTED_LINAR_FILE_NAME << endl;
00263     MeshImpl* linear_ex = new MeshImpl;
00264     linear_ex->read_vtk( EXPECTED_LINAR_FILE_NAME.c_str(), err );
00265     if( MSQ_CHKERR( err ) ) return 1;
00266 
00267     // Read in second copy of quadratic input mesh
00268     cout << "Reading " << QUADRATIC_INPUT_FILE_NAME << " again" << endl;
00269     MeshImpl* quadratic_in_2 = new MeshImpl;
00270     quadratic_in_2->read_vtk( QUADRATIC_INPUT_FILE_NAME.c_str(), err );
00271     if( MSQ_CHKERR( err ) ) return 1;
00272 
00273     // Read in expected quadratic results
00274     cout << "Reading " << EXPECTED_QUADRATIC_FILE_NAME << endl;
00275     MeshImpl* quadratic_ex = new MeshImpl;
00276     quadratic_ex->read_vtk( EXPECTED_QUADRATIC_FILE_NAME.c_str(), err );
00277     if( MSQ_CHKERR( err ) ) return 1;
00278 
00279     // Smooth linear mesh and check results
00280     cout << "Smoothing linear elements" << endl;
00281     InstructionQueue* q1 = create_instruction_queue( err );
00282     if( MSQ_CHKERR( err ) ) return 1;
00283     MeshDomainAssoc* mesh_and_domain = new MeshDomainAssoc( linear_in, &geom );
00284     q1->run_instructions( mesh_and_domain, err );
00285     if( MSQ_CHKERR( err ) ) return 1;
00286     delete mesh_and_domain;
00287     cout << "Checking results" << endl;
00288     compare_nodes( 0, NUM_CORNER_VERTICES, linear_in, linear_ex, err );
00289     if( MSQ_CHKERR( err ) )
00290     {
00291         MsqError tmperr;
00292         linear_in->write_vtk( "bad_mesh.vtk", tmperr );
00293         return 1;
00294     }
00295     delete q1;
00296     delete linear_in;
00297 
00298     // Smooth corner vertices and adjust mid-side nodes
00299     cout << "Smoothing quadratic elements" << endl;
00300     InstructionQueue* q3 = create_instruction_queue( err );
00301     if( MSQ_CHKERR( err ) ) return 1;
00302     if( !slave ) q3->set_slaved_ho_node_mode( Settings::SLAVE_NONE );
00303     //  q3->set_mapping_function( &quad9 );
00304     MeshDomainAssoc* mesh_and_domain2 = new MeshDomainAssoc( quadratic_in_2, &geom );
00305     q3->run_instructions( mesh_and_domain2, err );
00306     if( MSQ_CHKERR( err ) ) return 1;
00307     delete mesh_and_domain2;
00308     // Make sure corner vertices are the same as in the linear case
00309     cout << "Checking results" << endl;
00310     compare_nodes( 0, NUM_CORNER_VERTICES, quadratic_in_2, linear_ex, err );
00311     if( MSQ_CHKERR( err ) ) return 1;
00312     delete linear_ex;
00313 
00314     // Make sure mid-side vertices are updated correctly
00315     compare_nodes( NUM_CORNER_VERTICES, NUM_CORNER_VERTICES + NUM_MID_NODES, quadratic_in_2, quadratic_ex, err );
00316     if( MSQ_CHKERR( err ) )
00317     {
00318         MsqError tmperr;
00319         quadratic_in_2->write_vtk( "bad_mesh.vtk", tmperr );
00320         return 1;
00321     }
00322     delete q3;
00323     delete quadratic_in_2;
00324     delete quadratic_ex;
00325 
00326     if( MSQ_CHKERR( err ) ) return 1;
00327     return 0;
00328 }
00329 
00330 int do_smooth_ho()
00331 {
00332     MsqPrintError err( cout );
00333     //  QuadLagrangeShape quad9;
00334 
00335     // Create geometry
00336     PlanarDomain geom( PlanarDomain::XY );
00337 
00338     // Read in one copy of quadratic input mesh
00339     cout << "Reading " << HOUR_INPUT_FILE_NAME << endl;
00340     MeshImpl* quadratic_in = new MeshImpl;
00341     quadratic_in->read_vtk( HOUR_INPUT_FILE_NAME.c_str(), err );
00342     if( MSQ_CHKERR( err ) ) return 1;
00343 
00344     // Read in expected results
00345     // cout << "Reading " << HOUR_EXPECTED_FILE_NAME << endl;
00346     // MeshImpl* quadratic_ex = new MeshImpl;
00347     // quadratic_ex->read_vtk( results, err );
00348     // if (MSQ_CHKERR(err)) return 1;
00349 
00350     // Smooth linear mesh and check results
00351     cout << "Smoothing higher-order nodes" << endl;
00352     InstructionQueue* q1 = create_instruction_queue( err );
00353     if( MSQ_CHKERR( err ) ) return 1;
00354     q1->set_slaved_ho_node_mode( Settings::SLAVE_NONE );
00355     //  q1->set_mapping_function( &quad9 );
00356     MeshDomainAssoc mesh_and_domain = MeshDomainAssoc( quadratic_in, &geom );
00357     q1->run_instructions( &mesh_and_domain, err );
00358     if( MSQ_CHKERR( err ) ) return 1;
00359     cout << "Checking results" << endl;
00360     // compare_nodes( 0, NUM_CORNER_VERTICES + NUM_MID_NODES,
00361     //               quadratic_in, quadratic_ex, err );
00362     if( MSQ_CHKERR( err ) )
00363     {
00364         MsqError tmperr;
00365         quadratic_in->write_vtk( "bad_mesh.vtk", tmperr );
00366         return 1;
00367     }
00368     delete q1;
00369     quadratic_in->write_vtk( "smooth_ho.vtk", err );
00370     delete quadratic_in;
00371 
00372     if( MSQ_CHKERR( err ) ) return 1;
00373     return 0;
00374 }
00375 
00376 int main()
00377 {
00378     cout << "Running test with all higher-order nodes slaved." << endl;
00379     int result1 = do_test( true );
00380     cout << "Running test with no higher-order nodes slaved." << endl;
00381     int result2 = do_test( false );
00382     cout << "Running test with only ho-nodes free." << endl;
00383     int result3 = do_smooth_ho();
00384     return result1 + result2 + result3;
00385 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines