MOAB: Mesh Oriented datABase
(version 5.3.1)
|
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 <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, 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 }