MOAB: Mesh Oriented datABase  (version 5.4.1)
main.cpp
Go to the documentation of this file.
00001 /* *****************************************************************
00002     MESQUITE -- The Mesh Quality Improvement Toolkit
00003 
00004     Copyright 2010 Sandia National Laboratories.  Developed at the
00005     University of Wisconsin--Madison under SNL contract number
00006     624796.  The U.S. Government and the University of Wisconsin
00007     retain certain rights to 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     (2010) [email protected]
00024 
00025   ***************************************************************** */
00026 
00027 /** \file main.cpp
00028  *  \brief
00029  *  \author Jason Kraftcheck
00030  */
00031 
00032 #include "TestUtil.hpp"
00033 #include "MeshImpl.hpp"
00034 #include "PlanarDomain.hpp"
00035 #include "TriLagrangeShape.hpp"
00036 #include "QuadLagrangeShape.hpp"
00037 #include "IdealShapeTarget.hpp"
00038 #include "TShapeNB1.hpp"
00039 #include "TQualityMetric.hpp"
00040 #include "SteepestDescent.hpp"
00041 #include "SlaveBoundaryVertices.hpp"
00042 #include "PMeanPTemplate.hpp"
00043 #include "InstructionQueue.hpp"
00044 #include "MsqVertex.hpp"
00045 #include "MsqError.hpp"
00046 #include "IdealShapeTarget.hpp"
00047 #include "QualityAssessor.hpp"
00048 #include "PatchData.hpp"
00049 
00050 #include <iostream>
00051 #include <algorithm>
00052 #include <cassert>
00053 
00054 using namespace MBMesquite;
00055 
00056 std::string DEFAULT_INPUT_FILE = std::string( STRINGIFY( SRCDIR ) ) + "/input.vtk";
00057 
00058 void usage( const char* argv0 )
00059 {
00060     std::cout << "Usage: " << argv0 << " [<input_file>] [-o <output_file_base>]" << std::endl;
00061     exit( 1 );
00062 }
00063 
00064 int check_slaved_coords( Mesh& mesh,
00065                          Settings::HigherOrderSlaveMode mode,
00066                          std::string& name,
00067                          bool have_slaved_flag,
00068                          MsqError& err );
00069 
00070 // compare vertex coodinates between two topologically equivalent meshes
00071 int compare_node_coords( Mesh& mesh1, Mesh& mesh2, MsqError& err );
00072 
00073 // verify that no element corner node is marked as slaved
00074 int check_no_slaved_corners( Mesh& mesh, MsqError& err );
00075 
00076 // compare slaved flags in mesh to slaved state in global patch data
00077 int check_global_patch_slaved( Mesh& mesh, MsqError& err );
00078 
00079 // tag vertices marked as slaved in the PatchData
00080 void tag_patch_slaved( Mesh& mesh, Settings::HigherOrderSlaveMode mode, MsqError& err );
00081 
00082 int main( int argc, char* argv[] )
00083 {
00084     const char* input_file       = DEFAULT_INPUT_FILE.c_str();
00085     const char* output_file_base = 0;
00086     bool expect_output_base      = false;
00087     for( int i = 1; i < argc; ++i )
00088     {
00089         if( expect_output_base )
00090         {
00091             output_file_base   = argv[i];
00092             expect_output_base = false;
00093         }
00094         else if( !strcmp( argv[i], "-o" ) )
00095             expect_output_base = true;
00096         else if( DEFAULT_INPUT_FILE.compare( input_file ) )
00097             usage( argv[0] );
00098         else
00099             input_file = argv[i];
00100     }
00101     if( expect_output_base ) usage( argv[0] );
00102 
00103     MsqPrintError err( std::cerr );
00104     SlaveBoundaryVertices slaver( 1 );
00105     TShapeNB1 tmetric;
00106     IdealShapeTarget target;
00107     TQualityMetric metric( &target, &tmetric );
00108     PMeanPTemplate of( 1.0, &metric );
00109     SteepestDescent improver( &of );
00110     TerminationCriterion inner;
00111     inner.add_absolute_vertex_movement( 1e-3 );
00112     improver.set_inner_termination_criterion( &inner );
00113     QualityAssessor assess( &metric );
00114     InstructionQueue q;
00115     q.set_master_quality_improver( &improver, err );
00116     q.add_quality_assessor( &assess, err );
00117 
00118     TriLagrangeShape trishape;
00119     QuadLagrangeShape quadshape;
00120     q.set_mapping_function( &trishape );
00121     q.set_mapping_function( &quadshape );
00122 
00123     const int NUM_MODES = 4;
00124 
00125     Settings::HigherOrderSlaveMode modes[NUM_MODES] = { Settings::SLAVE_NONE, Settings::SLAVE_ALL,
00126                                                         Settings::SLAVE_CALCULATED, Settings::SLAVE_FLAG };
00127 
00128     std::string names[NUM_MODES] = { "NONE", "ALL", "CALCULATED", "FLAG" };
00129 
00130     MeshImpl meshes[NUM_MODES];
00131     std::vector< MeshDomainAssoc > meshes_and_domains;
00132 
00133     bool have_slaved_flag = true;
00134     std::vector< bool > flag( 1 );
00135     for( int i = 0; i < NUM_MODES; ++i )
00136     {
00137         std::cout << std::endl
00138                   << "-----------------------------------------------" << std::endl
00139                   << "     Mode: " << names[i] << std::endl
00140                   << "-----------------------------------------------" << std::endl;
00141 
00142         meshes[i].read_vtk( input_file, err );
00143         if( err ) return 1;
00144 
00145         if( modes[i] == Settings::SLAVE_CALCULATED )
00146         {
00147             q.add_vertex_slaver( &slaver, err );
00148         }
00149         else if( modes[i] == Settings::SLAVE_FLAG )
00150         {
00151             std::vector< Mesh::VertexHandle > verts;
00152             meshes[i].get_all_vertices( verts, err );
00153             if( err ) return 1;
00154             meshes[i].vertices_get_slaved_flag( arrptr( verts ), flag, 1, err );
00155             if( err )
00156             {
00157                 have_slaved_flag = false;
00158                 std::cout << "Skipped because input file does not contain slaved attribute" << std::endl;
00159                 err.clear();
00160                 continue;
00161             }
00162         }
00163 
00164         if( have_slaved_flag && modes[i] == Settings::SLAVE_FLAG )
00165         {
00166             if( !check_no_slaved_corners( meshes[i], err ) || err ) return 1;
00167             if( !check_global_patch_slaved( meshes[i], err ) || err ) return 1;
00168         }
00169 
00170         PlanarDomain plane;
00171         plane.fit_vertices( &meshes[i], err );
00172         if( err ) return 1;
00173 
00174         q.set_slaved_ho_node_mode( modes[i] );
00175         meshes_and_domains.push_back( MeshDomainAssoc( &meshes[i], &plane ) );
00176         q.run_instructions( &meshes_and_domains[i], err );
00177         if( err ) return 1;
00178 
00179         if( modes[i] == Settings::SLAVE_CALCULATED )
00180         {
00181             q.remove_vertex_slaver( &slaver, err );
00182         }
00183 
00184         if( output_file_base )
00185         {
00186             tag_patch_slaved( meshes[i], modes[i], err );
00187             std::string name( output_file_base );
00188             name += "-";
00189             name += names[i];
00190             name += ".vtk";
00191             meshes[i].write_vtk( name.c_str(), err );
00192             if( err ) return 1;
00193         }
00194     }
00195 
00196     int exit_code = 0;
00197     if( !DEFAULT_INPUT_FILE.compare( input_file ) )
00198     {
00199         for( int i = 0; i < NUM_MODES; ++i )
00200         {
00201             std::cout << std::endl
00202                       << "-----------------------------------------------" << std::endl
00203                       << "     Mode: " << names[i] << std::endl
00204                       << "-----------------------------------------------" << std::endl;
00205 
00206             exit_code += check_slaved_coords( meshes[i], modes[i], names[i], have_slaved_flag, err );
00207             if( err ) return 1;
00208         }
00209 
00210         // flags should correspond to same slaved nodes as calculated,
00211         // so resulting meshes should be identical.
00212         if( have_slaved_flag )
00213         {
00214             int flag_idx = std::find( modes, modes + NUM_MODES, Settings::SLAVE_FLAG ) - modes;
00215             int calc_idx = std::find( modes, modes + NUM_MODES, Settings::SLAVE_CALCULATED ) - modes;
00216             exit_code += compare_node_coords( meshes[flag_idx], meshes[calc_idx], err );
00217             if( err ) return 1;
00218         }
00219     }
00220 
00221     return exit_code;
00222 }
00223 
00224 Vector3D get_slaved_coords( Mesh& mesh, Mesh::VertexHandle vertex, MsqError& err )
00225 {
00226     std::vector< Mesh::ElementHandle > elem;
00227     std::vector< size_t > off;
00228     mesh.vertices_get_attached_elements( &vertex, 1, elem, off, err );
00229     if( MSQ_CHKERR( err ) ) return Vector3D( 0.0 );
00230 
00231     std::vector< Mesh::VertexHandle > verts;
00232     mesh.elements_get_attached_vertices( arrptr( elem ), 1, verts, off, err );
00233     if( MSQ_CHKERR( err ) ) return Vector3D( 0.0 );
00234 
00235     EntityTopology type;
00236     mesh.elements_get_topologies( arrptr( elem ), &type, 1, err );
00237     if( MSQ_CHKERR( err ) ) return Vector3D( 0.0 );
00238 
00239     size_t idx = std::find( verts.begin(), verts.end(), vertex ) - verts.begin();
00240     unsigned side_dim, side_num;
00241     TopologyInfo::side_from_higher_order( type, verts.size(), idx, side_dim, side_num, err );
00242     if( MSQ_CHKERR( err ) ) return Vector3D( 0.0 );
00243 
00244     // just return the mean of the corner vertices defining the side.
00245     // this should always be correct for mid-edge nodes but is a bit
00246     // dubious for mid-face or mid-region nodes.  But our default input
00247     // file (the only file for which we should be in this function) will
00248     // have only mid-edge nodes.
00249     unsigned n;
00250     const unsigned* side_vtx = TopologyInfo::side_vertices( type, side_dim, side_num, n, err );
00251     if( MSQ_CHKERR( err ) ) return Vector3D( 0.0 );
00252 
00253     Vector3D sum( 0.0 );
00254     for( unsigned i = 0; i < n; ++i )
00255     {
00256         MsqVertex coords;
00257         Mesh::VertexHandle vtx = verts[side_vtx[i]];
00258         mesh.vertices_get_coordinates( &vtx, &coords, 1, err );
00259         if( MSQ_CHKERR( err ) ) return Vector3D( 0.0 );
00260         sum += coords;
00261     }
00262     return sum / n;
00263 }
00264 
00265 int check_slaved_coords( Mesh& mesh,
00266                          Settings::HigherOrderSlaveMode mode,
00267                          std::string& name,
00268                          bool have_slaved_flag,
00269                          MsqError& err )
00270 {
00271     const double EPSILON = 1e-4;
00272 
00273     // We can distinguish between nodes near the boundary and those
00274     // further inside based on the slaved attribute flag in the
00275     // the default input file.  The default input file should always
00276     // have it defined
00277     if( !have_slaved_flag )
00278     {
00279         std::cerr << "slaved flag not specified in input file. Cannot validate results." << std::endl;
00280         return 1;
00281     }
00282 
00283     // Get list of all higher-order nodes
00284     std::vector< Mesh::ElementHandle > elems;
00285     std::vector< Mesh::VertexHandle > verts;
00286     std::vector< size_t > offsets;
00287     mesh.get_all_elements( elems, err );
00288     MSQ_ERRZERO( err );
00289     std::vector< EntityTopology > types( elems.size() );
00290     mesh.elements_get_topologies( arrptr( elems ), arrptr( types ), elems.size(), err );
00291     MSQ_ERRZERO( err );
00292     mesh.elements_get_attached_vertices( arrptr( elems ), elems.size(), verts, offsets, err );
00293     MSQ_ERRZERO( err );
00294     std::vector< Mesh::VertexHandle >::iterator r, e, w = verts.begin();
00295     for( size_t i = 0; i < elems.size(); ++i )
00296     {
00297         r = verts.begin() + offsets[i] + TopologyInfo::corners( types[i] );
00298         e = verts.begin() + offsets[i + 1];
00299         w = std::copy( r, e, w );
00300     }
00301     std::sort( verts.begin(), w );
00302     verts.erase( std::unique( verts.begin(), w ), verts.end() );
00303 
00304     // Get lists of slaved and non-slaved free vertices, where 'slaved'
00305     // are those vertices far from the boundary that would be slaved if
00306     // mode == SLAVE_FLAG, not those that were actually slaved during
00307     // the optimization
00308     std::vector< bool > fixed, slaved;
00309     mesh.vertices_get_fixed_flag( arrptr( verts ), fixed, verts.size(), err );
00310     if( MSQ_CHKERR( err ) )
00311     {
00312         return 1;
00313     }
00314     mesh.vertices_get_slaved_flag( arrptr( verts ), slaved, verts.size(), err );
00315     if( MSQ_CHKERR( err ) )
00316     {
00317         return 1;
00318     }
00319     std::vector< Mesh::VertexHandle > free, slave;
00320     for( size_t i = 0; i < verts.size(); ++i )
00321     {
00322         if( !fixed[i] && slaved[i] )
00323             slave.push_back( verts[i] );
00324         else if( !fixed[i] && !slaved[i] )
00325             free.push_back( verts[i] );
00326     }
00327 
00328     // get all coordinates
00329     std::vector< MsqVertex > free_coords( free.size() ), slave_coords( slave.size() );
00330     mesh.vertices_get_coordinates( arrptr( free ), arrptr( free_coords ), free.size(), err );
00331     MSQ_ERRZERO( err );
00332     mesh.vertices_get_coordinates( arrptr( slave ), arrptr( slave_coords ), slave.size(), err );
00333     MSQ_ERRZERO( err );
00334 
00335     int error_count = 0;
00336 
00337     // Expect vertices near boundary to be at slave vertex positions
00338     // if mode was SLAVE_ALL
00339     if( mode == Settings::SLAVE_ALL )
00340     {
00341         for( size_t i = 0; i < free.size(); ++i )
00342         {
00343             Vector3D exp = get_slaved_coords( mesh, free[i], err );
00344             MSQ_ERRZERO( err );
00345             exp -= free_coords[i];
00346             if( exp.length() > EPSILON )
00347             {
00348                 std::cerr << "Slaved vertex " << (size_t)free[i] << " not at expected slaved location" << std::endl;
00349                 ++error_count;
00350             }
00351         }
00352     }
00353     // Otherwise, given the default input mesh, at least some of the
00354     // vertices should be somewhere other than the slaved position
00355     else
00356     {
00357         int not_at_slaved_count = 0;
00358 
00359         for( size_t i = 0; i < free.size(); ++i )
00360         {
00361             Vector3D exp = get_slaved_coords( mesh, free[i], err );
00362             MSQ_ERRZERO( err );
00363             exp -= free_coords[i];
00364             if( exp.length() > EPSILON )
00365             {
00366                 ++not_at_slaved_count;
00367             }
00368         }
00369 
00370         if( 0 == not_at_slaved_count )
00371         {
00372             std::cerr << "All non-slaved vertices at slaved vertex locations" << std::endl;
00373             error_count += free.size();
00374         }
00375     }
00376 
00377     // expect all interior vertices to be at roughly the slaved location
00378     if( mode != Settings::SLAVE_NONE )
00379     {
00380         for( size_t i = 0; i < slave.size(); ++i )
00381         {
00382             Vector3D exp = get_slaved_coords( mesh, slave[i], err );
00383             MSQ_ERRZERO( err );
00384             exp -= slave_coords[i];
00385             if( exp.length() > EPSILON )
00386             {
00387                 std::cerr << "Interior vertex " << (size_t)slave[i] << " not at expected slaved location" << std::endl;
00388                 ++error_count;
00389             }
00390         }
00391     }
00392 
00393     return error_count;
00394 }
00395 
00396 int compare_node_coords( Mesh& mesh1, Mesh& mesh2, MsqError& err )
00397 {
00398     const double EPSILON = 1e-4;
00399 
00400     std::vector< Mesh::VertexHandle > verts1, verts2;
00401     mesh1.get_all_vertices( verts1, err );
00402     MSQ_ERRZERO( err );
00403     mesh2.get_all_vertices( verts2, err );
00404     MSQ_ERRZERO( err );
00405     std::vector< MsqVertex > coords1( verts1.size() ), coords2( verts2.size() );
00406     mesh1.vertices_get_coordinates( arrptr( verts1 ), arrptr( coords1 ), verts1.size(), err );
00407     MSQ_ERRZERO( err );
00408     mesh2.vertices_get_coordinates( arrptr( verts2 ), arrptr( coords2 ), verts2.size(), err );
00409     MSQ_ERRZERO( err );
00410 
00411     int error_count = 0;
00412     assert( verts1.size() == verts2.size() );
00413     for( size_t i = 0; i < verts1.size(); ++i )
00414     {
00415         assert( verts1[i] == verts2[i] );
00416         Vector3D diff = coords1[i] - coords2[i];
00417         if( diff.length() > EPSILON )
00418         {
00419             std::cerr << "Vertex coordinates differ between calculated and flagged "
00420                          "meshes for vertex "
00421                       << (size_t)verts1[i] << std::endl;
00422             ++error_count;
00423         }
00424     }
00425 
00426     return error_count;
00427 }
00428 
00429 int check_no_slaved_corners( Mesh& mesh, MsqError& err )
00430 {
00431     std::vector< Mesh::ElementHandle > elems;
00432     std::vector< Mesh::VertexHandle > verts;
00433     std::vector< EntityTopology > types;
00434     std::vector< size_t > offsets;
00435     mesh.get_all_elements( elems, err );
00436     MSQ_ERRZERO( err );
00437     types.resize( elems.size() );
00438     mesh.elements_get_topologies( arrptr( elems ), arrptr( types ), elems.size(), err );
00439     MSQ_ERRZERO( err );
00440     mesh.elements_get_attached_vertices( arrptr( elems ), elems.size(), verts, offsets, err );
00441     MSQ_ERRZERO( err );
00442 
00443     std::vector< bool > slaved;
00444     mesh.vertices_get_slaved_flag( arrptr( verts ), slaved, verts.size(), err );
00445     MSQ_ERRZERO( err );
00446 
00447     int error_count = 0;
00448     for( size_t i = 0; i < elems.size(); ++i )
00449     {
00450         unsigned n = TopologyInfo::corners( types[i] );
00451         for( unsigned j = 0; j < n; ++j )
00452         {
00453             if( slaved[offsets[i] + j] )
00454             {
00455                 std::cerr << "Element " << (size_t)elems[i] << " corner " << j << " is slaved" << std::endl;
00456                 ++error_count;
00457             }
00458         }
00459     }
00460 
00461     return error_count == 0;
00462 }
00463 
00464 int check_global_patch_slaved( Mesh& mesh, MsqError& err )
00465 {
00466     Settings s;
00467     s.set_slaved_ho_node_mode( Settings::SLAVE_FLAG );
00468     MeshDomainAssoc mesh_and_domain = MeshDomainAssoc( &mesh, 0 );
00469     Instruction::initialize_vertex_byte( &mesh_and_domain, &s, err );
00470     MSQ_ERRZERO( err );
00471 
00472     PatchData pd;
00473     pd.attach_settings( &s );
00474     pd.set_mesh( &mesh );
00475     pd.fill_global_patch( err );
00476     MSQ_ERRZERO( err );
00477 
00478     std::vector< bool > fixed, slaved;
00479     mesh.vertices_get_fixed_flag( pd.get_vertex_handles_array(), fixed, pd.num_nodes(), err );
00480     MSQ_ERRZERO( err );
00481     mesh.vertices_get_slaved_flag( pd.get_vertex_handles_array(), slaved, pd.num_nodes(), err );
00482     MSQ_ERRZERO( err );
00483 
00484     const size_t first_free   = 0;
00485     const size_t first_slaved = pd.num_free_vertices();
00486     const size_t first_fixed  = pd.num_free_vertices() + pd.num_slave_vertices();
00487     int error_count           = 0;
00488     for( size_t i = first_free; i < first_slaved; ++i )
00489     {
00490         if( fixed[i] )
00491         {
00492             std::cerr << "Vertex " << (size_t)pd.get_vertex_handles_array()[i]
00493                       << " is fixed in mesh but free in PatchData" << std::endl;
00494             ++error_count;
00495         }
00496         if( slaved[i] )
00497         {
00498             std::cerr << "Vertex " << (size_t)pd.get_vertex_handles_array()[i]
00499                       << " is slaved in mesh but free in PatchData" << std::endl;
00500             ++error_count;
00501         }
00502     }
00503     for( size_t i = first_slaved; i < first_fixed; ++i )
00504     {
00505         if( fixed[i] )
00506         {
00507             std::cerr << "Vertex " << (size_t)pd.get_vertex_handles_array()[i]
00508                       << " is fixed in mesh but slaved in PatchData" << std::endl;
00509             ++error_count;
00510         }
00511         else if( !slaved[i] )
00512         {
00513             std::cerr << "Vertex " << (size_t)pd.get_vertex_handles_array()[i]
00514                       << " is free in Mesh but slaved in PatchData" << std::endl;
00515             ++error_count;
00516         }
00517     }
00518     for( size_t i = first_fixed; i < pd.num_nodes(); ++i )
00519     {
00520         if( !fixed[i] )
00521         {
00522             std::cerr << "Vertex " << (size_t)pd.get_vertex_handles_array()[i]
00523                       << " is not fixed in mesh but is in PatchData" << std::endl;
00524             ++error_count;
00525         }
00526     }
00527     return 0 == error_count;
00528 }
00529 
00530 void tag_patch_slaved( Mesh& mesh, Settings::HigherOrderSlaveMode mode, MsqError& err )
00531 {
00532     int zero      = 0;
00533     TagHandle tag = mesh.tag_create( "pd_slaved", Mesh::INT, 1, &zero, err );MSQ_ERRRTN( err );
00534 
00535     Settings s;
00536     s.set_slaved_ho_node_mode( mode );
00537     PatchData pd;
00538     pd.attach_settings( &s );
00539     pd.set_mesh( &mesh );
00540     pd.fill_global_patch( err );MSQ_ERRRTN( err );
00541 
00542     const Mesh::VertexHandle* verts = pd.get_vertex_handles_array() + pd.num_free_vertices();
00543     std::vector< int > ones( pd.num_slave_vertices(), 1 );
00544     mesh.tag_set_vertex_data( tag, pd.num_slave_vertices(), verts, arrptr( ones ), err );MSQ_ERRRTN( err );
00545 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines