MOAB: Mesh Oriented datABase  (version 5.4.1)
setup.cpp
Go to the documentation of this file.
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 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines