MOAB: Mesh Oriented datABase  (version 5.2.1)
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" : i % 10 == 1 ? "nd" : i % 10 == 2 ? "rd" : "th", verts1[i][0], verts1[i][1],
00190               verts1[i][2], verts2[i][0], verts2[i][1], verts2[i][2] );
00191             return;
00192         }
00193     }
00194 }
00195 
00196 // code copied from testSuite/algorithm_test/main.cpp
00197 InstructionQueue* create_instruction_queue( MsqError& err )
00198 {
00199     // creates an intruction queue
00200     InstructionQueue* queue1 = new InstructionQueue;
00201 
00202     // creates a mean ratio quality metric ...
00203     // IdealWeightInverseMeanRatio* mean = new IdealWeightInverseMeanRatio(err); MSQ_ERRZERO(err);
00204     TargetCalculator* tc = new IdealShapeTarget;
00205     TQualityMetric* mean = new TQualityMetric( tc, 0, new TInverseMeanRatio );
00206 
00207     LPtoPTemplate* obj_func = new LPtoPTemplate( mean, 1, err );
00208     MSQ_ERRZERO( err );
00209 
00210     // creates the optimization procedures
00211     //   ConjugateGradient* pass1 = new ConjugateGradient( obj_func, err );
00212     SteepestDescent* pass1 = new SteepestDescent( obj_func );
00213 
00214     // perform optimization globally
00215     pass1->use_global_patch();
00216 
00217     // QualityAssessor* mean_qa = new QualityAssessor(mean);
00218 
00219     //**************Set termination criterion****************
00220 
00221     // perform 1 pass of the outer loop (this line isn't essential as it is
00222     // the default behavior).
00223     TerminationCriterion* tc_outer = new TerminationCriterion;
00224     tc_outer->add_iteration_limit( 1 );
00225     pass1->set_outer_termination_criterion( tc_outer );
00226 
00227     // perform the inner loop until a certain objective function value is
00228     // reached.  The exact value needs to be determined (about 18095).
00229     // As a safety, also stop if the time exceeds 10 minutes (600 seconds).
00230     TerminationCriterion* tc_inner = new TerminationCriterion;
00231     tc_inner->add_absolute_vertex_movement( 1e-6 );
00232     pass1->set_inner_termination_criterion( tc_inner );
00233 
00234     // adds 1 pass of pass1 to mesh_set1
00235     // queue1->add_quality_assessor(mean_qa,err); MSQ_ERRZERO(err);
00236     queue1->set_master_quality_improver( pass1, err );
00237     MSQ_ERRZERO( err );
00238     // queue1->add_quality_assessor(mean_qa,err); MSQ_ERRZERO(err);
00239 
00240     return queue1;
00241 }
00242 
00243 int do_test( bool slave )
00244 {
00245     MsqPrintError err( cout );
00246     //  QuadLagrangeShape quad9;
00247 
00248     // Create geometry
00249     Vector3D z( 0, 0, 1 ), o( 0, 0, 0 );
00250     PlanarDomain geom( z, o );
00251 
00252     // Read in linear input mesh
00253     cout << "Reading " << LINEAR_INPUT_FILE_NAME << endl;
00254     MeshImpl* linear_in = new MeshImpl;
00255     linear_in->read_vtk( LINEAR_INPUT_FILE_NAME.c_str(), err );
00256     if( MSQ_CHKERR( err ) ) return 1;
00257 
00258     // Read in expected linear results
00259     cout << "Reading " << EXPECTED_LINAR_FILE_NAME << endl;
00260     MeshImpl* linear_ex = new MeshImpl;
00261     linear_ex->read_vtk( EXPECTED_LINAR_FILE_NAME.c_str(), err );
00262     if( MSQ_CHKERR( err ) ) return 1;
00263 
00264     // Read in second copy of quadratic input mesh
00265     cout << "Reading " << QUADRATIC_INPUT_FILE_NAME << " again" << endl;
00266     MeshImpl* quadratic_in_2 = new MeshImpl;
00267     quadratic_in_2->read_vtk( QUADRATIC_INPUT_FILE_NAME.c_str(), err );
00268     if( MSQ_CHKERR( err ) ) return 1;
00269 
00270     // Read in expected quadratic results
00271     cout << "Reading " << EXPECTED_QUADRATIC_FILE_NAME << endl;
00272     MeshImpl* quadratic_ex = new MeshImpl;
00273     quadratic_ex->read_vtk( EXPECTED_QUADRATIC_FILE_NAME.c_str(), err );
00274     if( MSQ_CHKERR( err ) ) return 1;
00275 
00276     // Smooth linear mesh and check results
00277     cout << "Smoothing linear elements" << endl;
00278     InstructionQueue* q1 = create_instruction_queue( err );
00279     if( MSQ_CHKERR( err ) ) return 1;
00280     MeshDomainAssoc* mesh_and_domain = new MeshDomainAssoc( linear_in, &geom );
00281     q1->run_instructions( mesh_and_domain, err );
00282     if( MSQ_CHKERR( err ) ) return 1;
00283     delete mesh_and_domain;
00284     cout << "Checking results" << endl;
00285     compare_nodes( 0, NUM_CORNER_VERTICES, linear_in, linear_ex, err );
00286     if( MSQ_CHKERR( err ) )
00287     {
00288         MsqError tmperr;
00289         linear_in->write_vtk( "bad_mesh.vtk", tmperr );
00290         return 1;
00291     }
00292     delete q1;
00293     delete linear_in;
00294 
00295     // Smooth corner vertices and adjust mid-side nodes
00296     cout << "Smoothing quadratic elements" << endl;
00297     InstructionQueue* q3 = create_instruction_queue( err );
00298     if( MSQ_CHKERR( err ) ) return 1;
00299     if( !slave ) q3->set_slaved_ho_node_mode( Settings::SLAVE_NONE );
00300     //  q3->set_mapping_function( &quad9 );
00301     MeshDomainAssoc* mesh_and_domain2 = new MeshDomainAssoc( quadratic_in_2, &geom );
00302     q3->run_instructions( mesh_and_domain2, err );
00303     if( MSQ_CHKERR( err ) ) return 1;
00304     delete mesh_and_domain2;
00305     // Make sure corner vertices are the same as in the linear case
00306     cout << "Checking results" << endl;
00307     compare_nodes( 0, NUM_CORNER_VERTICES, quadratic_in_2, linear_ex, err );
00308     if( MSQ_CHKERR( err ) ) return 1;
00309     delete linear_ex;
00310 
00311     // Make sure mid-side vertices are updated correctly
00312     compare_nodes( NUM_CORNER_VERTICES, NUM_CORNER_VERTICES + NUM_MID_NODES, quadratic_in_2, quadratic_ex, err );
00313     if( MSQ_CHKERR( err ) )
00314     {
00315         MsqError tmperr;
00316         quadratic_in_2->write_vtk( "bad_mesh.vtk", tmperr );
00317         return 1;
00318     }
00319     delete q3;
00320     delete quadratic_in_2;
00321     delete quadratic_ex;
00322 
00323     if( MSQ_CHKERR( err ) ) return 1;
00324     return 0;
00325 }
00326 
00327 int do_smooth_ho()
00328 {
00329     MsqPrintError err( cout );
00330     //  QuadLagrangeShape quad9;
00331 
00332     // Create geometry
00333     PlanarDomain geom( PlanarDomain::XY );
00334 
00335     // Read in one copy of quadratic input mesh
00336     cout << "Reading " << HOUR_INPUT_FILE_NAME << endl;
00337     MeshImpl* quadratic_in = new MeshImpl;
00338     quadratic_in->read_vtk( HOUR_INPUT_FILE_NAME.c_str(), err );
00339     if( MSQ_CHKERR( err ) ) return 1;
00340 
00341     // Read in expected results
00342     // cout << "Reading " << HOUR_EXPECTED_FILE_NAME << endl;
00343     // MeshImpl* quadratic_ex = new MeshImpl;
00344     // quadratic_ex->read_vtk( results, err );
00345     // if (MSQ_CHKERR(err)) return 1;
00346 
00347     // Smooth linear mesh and check results
00348     cout << "Smoothing higher-order nodes" << endl;
00349     InstructionQueue* q1 = create_instruction_queue( err );
00350     if( MSQ_CHKERR( err ) ) return 1;
00351     q1->set_slaved_ho_node_mode( Settings::SLAVE_NONE );
00352     //  q1->set_mapping_function( &quad9 );
00353     MeshDomainAssoc mesh_and_domain = MeshDomainAssoc( quadratic_in, &geom );
00354     q1->run_instructions( &mesh_and_domain, err );
00355     if( MSQ_CHKERR( err ) ) return 1;
00356     cout << "Checking results" << endl;
00357     // compare_nodes( 0, NUM_CORNER_VERTICES + NUM_MID_NODES,
00358     //               quadratic_in, quadratic_ex, err );
00359     if( MSQ_CHKERR( err ) )
00360     {
00361         MsqError tmperr;
00362         quadratic_in->write_vtk( "bad_mesh.vtk", tmperr );
00363         return 1;
00364     }
00365     delete q1;
00366     quadratic_in->write_vtk( "smooth_ho.vtk", err );
00367     delete quadratic_in;
00368 
00369     if( MSQ_CHKERR( err ) ) return 1;
00370     return 0;
00371 }
00372 
00373 int main()
00374 {
00375     cout << "Running test with all higher-order nodes slaved." << endl;
00376     int result1 = do_test( true );
00377     cout << "Running test with no higher-order nodes slaved." << endl;
00378     int result2 = do_test( false );
00379     cout << "Running test with only ho-nodes free." << endl;
00380     int result3 = do_smooth_ho();
00381     return result1 + result2 + result3;
00382 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines