MOAB: Mesh Oriented datABase
(version 5.4.1)
|
00001 #include "MeshImpl.hpp" 00002 #include "MsqError.hpp" 00003 #include "MsqVertex.hpp" 00004 00005 #include <iostream> 00006 #include <cstdlib> 00007 #include <cassert> 00008 #include <map> 00009 #include <vector> 00010 #include <algorithm> 00011 00012 const double PERTURB_FRACT = 0.4; 00013 00014 using namespace MBMesquite; 00015 00016 void usage( const char* argv0 ) 00017 { 00018 std::cerr << "Usage: " << argv0 << " <depth> <input_file> <output_file>" << std::endl; 00019 exit( 1 ); 00020 } 00021 00022 // Routine to create initial mesh for test. 00023 // o Marks vertices at a greater topological depth than the specified 00024 // value as slaved. 00025 // o Perturbs higher-order vertices on skin towards element center 00026 // o Marks skin vertices as fixed 00027 int main( int argc, char* argv[] ) 00028 { 00029 if( argc != 4 ) usage( argv[0] ); 00030 00031 char* endptr = 0; 00032 const long n = strtol( argv[1], &endptr, 0 ); 00033 if( *endptr || n < 0 ) usage( argv[0] ); 00034 00035 // read input mesh 00036 MeshImpl mesh; 00037 MsqPrintError err( std::cerr ); 00038 mesh.read_vtk( argv[2], err ); 00039 if( err ) return 1; 00040 00041 // get skin vertices 00042 mesh.mark_skin_fixed( err, true ); 00043 if( err ) return 1; 00044 std::vector< Mesh::VertexHandle > verts; 00045 mesh.get_all_vertices( verts, err ); 00046 if( err ) return 1; 00047 std::vector< bool > fixed; 00048 mesh.vertices_get_fixed_flag( arrptr( verts ), fixed, verts.size(), err ); 00049 if( err ) return 1; 00050 std::vector< Mesh::VertexHandle > skin; 00051 for( size_t i = 0; i < verts.size(); ++i ) 00052 if( fixed[i] ) skin.push_back( verts[i] ); 00053 00054 // create map for vertex depth, and initialize to 0 for skin vertices 00055 std::map< Mesh::VertexHandle, int > depth; 00056 std::map< Mesh::VertexHandle, int >::iterator d_iter; 00057 for( size_t i = 0; i < skin.size(); ++i ) 00058 depth[skin[i]] = 0; 00059 00060 // get all elements 00061 std::vector< Mesh::ElementHandle > curr, next; 00062 std::vector< Mesh::ElementHandle > conn; 00063 std::vector< size_t > off; 00064 mesh.get_all_elements( next, err ); 00065 00066 // build sorted list of higher-order vertices 00067 std::vector< Mesh::VertexHandle > higher_order; 00068 for( size_t i = 0; i < next.size(); ++i ) 00069 { 00070 Mesh::ElementHandle elem = next[i]; 00071 conn.clear(); 00072 mesh.elements_get_attached_vertices( &elem, 1, conn, off, err ); 00073 if( err ) return 1; 00074 EntityTopology type; 00075 mesh.elements_get_topologies( &elem, &type, 1, err ); 00076 std::copy( conn.begin() + TopologyInfo::corners( type ), conn.end(), std::back_inserter( higher_order ) ); 00077 } 00078 std::sort( higher_order.begin(), higher_order.end() ); 00079 higher_order.erase( std::unique( higher_order.begin(), higher_order.end() ), higher_order.end() ); 00080 00081 // build depth map for all vertices 00082 while( !next.empty() ) 00083 { 00084 curr.swap( next ); 00085 next.clear(); 00086 while( !curr.empty() ) 00087 { 00088 Mesh::ElementHandle elem = curr.back(); 00089 curr.pop_back(); 00090 00091 conn.clear(); 00092 mesh.elements_get_attached_vertices( &elem, 1, conn, off, err ); 00093 if( err ) return 1; 00094 00095 int min = std::numeric_limits< int >::max(); 00096 for( size_t i = 0; i < conn.size(); ++i ) 00097 { 00098 d_iter = depth.find( conn[i] ); 00099 if( d_iter != depth.end() && d_iter->second < min ) min = d_iter->second; 00100 } 00101 00102 if( min == std::numeric_limits< int >::max() ) 00103 { 00104 next.push_back( elem ); 00105 continue; 00106 } 00107 00108 for( size_t i = 0; i < conn.size(); ++i ) 00109 { 00110 d_iter = depth.find( conn[i] ); 00111 00112 if( d_iter == depth.end() || d_iter->second > min + 1 ) depth[conn[i]] = min + 1; 00113 } 00114 } 00115 } 00116 00117 // write depth map to tag for debugging purposes 00118 std::vector< int > depth_vals( verts.size() ); 00119 for( size_t i = 0; i < verts.size(); ++i ) 00120 depth_vals[i] = depth[verts[i]]; 00121 TagHandle tag = mesh.tag_create( "depth", Mesh::INT, 1, 0, err ); 00122 if( err ) return 1; 00123 mesh.tag_set_vertex_data( tag, verts.size(), arrptr( verts ), arrptr( depth_vals ), err ); 00124 if( err ) return 1; 00125 00126 // set tag specifying slaved vertices 00127 for( size_t i = 0; i < verts.size(); ++i ) 00128 if( std::binary_search( higher_order.begin(), higher_order.end(), verts[i] ) ) 00129 depth_vals[i] = depth[verts[i]] > n; 00130 else 00131 depth_vals[i] = 0; 00132 tag = mesh.tag_create( "slaved", Mesh::INT, 1, 0, err ); 00133 if( err ) return 1; 00134 mesh.tag_set_vertex_data( tag, verts.size(), arrptr( verts ), arrptr( depth_vals ), err ); 00135 if( err ) return 1; 00136 00137 // perturb mid-edge nodes along boundary 00138 std::vector< MsqVertex > coords; 00139 for( size_t i = 0; i < skin.size(); ++i ) 00140 { 00141 if( !std::binary_search( higher_order.begin(), higher_order.end(), skin[i] ) ) continue; 00142 00143 curr.clear(); 00144 mesh.vertices_get_attached_elements( &skin[i], 1, curr, off, err ); 00145 if( err ) return 1; 00146 assert( curr.size() == 1 ); 00147 conn.clear(); 00148 mesh.elements_get_attached_vertices( arrptr( curr ), 1, conn, off, err ); 00149 if( err ) return 1; 00150 00151 // estimate element center 00152 coords.resize( conn.size() ); 00153 mesh.vertices_get_coordinates( arrptr( conn ), arrptr( coords ), conn.size(), err ); 00154 if( err ) return 1; 00155 00156 Vector3D mean( 0.0 ); 00157 for( size_t j = 0; j < coords.size(); ++j ) 00158 mean += coords[j]; 00159 mean /= coords.size(); 00160 00161 size_t idx = std::find( conn.begin(), conn.end(), skin[i] ) - conn.begin(); 00162 assert( idx < conn.size() ); 00163 Vector3D init = coords[idx]; 00164 Vector3D pos = ( 1 - PERTURB_FRACT ) * init + PERTURB_FRACT * mean; 00165 mesh.vertex_set_coordinates( skin[i], pos, err ); 00166 if( err ) return 1; 00167 } 00168 00169 mesh.write_vtk( argv[3], err ); 00170 if( err ) return 1; 00171 00172 return 0; 00173 }