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