MOAB: Mesh Oriented datABase  (version 5.2.1)
ExtrudePoly.cpp
Go to the documentation of this file.
00001 /** @example ExtrudePoly.cpp
00002  * Description: read a 2d mesh in plane, extrude to form  layers of prisms (polyhedra) \n
00003  *
00004  * To run: ./ExtrudePoly [meshfile] [outfile] [nlayers] [thickness_layer] \n
00005  */
00006 
00007 #include "moab/Core.hpp"
00008 #include "moab/ReadUtilIface.hpp"
00009 #include <iostream>
00010 
00011 using namespace moab;
00012 using namespace std;
00013 
00014 #ifndef MESH_DIR
00015 #define MESH_DIR "."
00016 #endif
00017 
00018 // at most 20 edges per polygon
00019 #define MAXEDGES 20
00020 // Note: change the file name below to test a trivial "No such file or directory" error
00021 string test_file_name = string( MESH_DIR ) + string( "/io/poly8-10.vtk" );
00022 string output         = string( "polyhedra.vtk" );
00023 int layers            = 1;
00024 double layer_thick    = 1.0;
00025 
00026 int main( int argc, char** argv )
00027 {
00028     // Get MOAB instance
00029     Interface* mb = new( std::nothrow ) Core;
00030     if( NULL == mb ) return 1;
00031 
00032     // Need option handling here for input filename
00033     if( argc > 1 )
00034     {
00035         // User has input a mesh file
00036         test_file_name = argv[1];
00037     }
00038     if( argc > 2 ) output = argv[2];
00039 
00040     if( argc > 3 ) layers = atoi( argv[3] );
00041 
00042     if( argc > 4 ) layer_thick = atof( argv[4] );
00043 
00044     std::cout << "Run: " << argv[0] << " " << test_file_name << " " << output << " " << layers << " " << layer_thick
00045               << "\n";
00046     // Load the mesh from vtk file
00047     ErrorCode rval = mb->load_mesh( test_file_name.c_str() );MB_CHK_ERR( rval );
00048 
00049     // Get verts entities, by type
00050     Range verts;
00051     rval = mb->get_entities_by_type( 0, MBVERTEX, verts );MB_CHK_ERR( rval );
00052 
00053     // Get faces, by dimension, so we stay generic to entity type
00054     Range faces;
00055     rval = mb->get_entities_by_dimension( 0, 2, faces );MB_CHK_ERR( rval );
00056     cout << "Number of vertices is " << verts.size() << endl;
00057     cout << "Number of faces is " << faces.size() << endl;
00058 
00059     Range edges;
00060     // create all edges
00061     rval = mb->get_adjacencies( faces, 1, true, edges, Interface::UNION );MB_CHK_ERR( rval );
00062 
00063     cout << "Number of edges is " << edges.size() << endl;
00064 
00065     std::vector< double > coords;
00066     int nvPerLayer = (int)verts.size();
00067     coords.resize( 3 * nvPerLayer );
00068 
00069     rval = mb->get_coords( verts, &coords[0] );MB_CHK_ERR( rval );
00070     // create first vertices
00071     Range* newVerts = new Range[layers + 1];
00072     newVerts[0]     = verts;  // just for convenience
00073     for( int ii = 0; ii < layers; ii++ )
00074     {
00075         for( int i = 0; i < nvPerLayer; i++ )
00076             coords[3 * i + 2] += layer_thick;
00077 
00078         rval = mb->create_vertices( &coords[0], nvPerLayer, newVerts[ii + 1] );MB_CHK_ERR( rval );
00079     }
00080     // for each edge, we will create layers quads
00081     int nquads = edges.size() * layers;
00082     ReadUtilIface* read_iface;
00083     rval = mb->query_interface( read_iface );MB_CHK_SET_ERR( rval, "Error in query_interface" );
00084 
00085     EntityHandle start_elem, *connect;
00086     // Create quads
00087     rval = read_iface->get_element_connect( nquads, 4, MBQUAD, 0, start_elem, connect );MB_CHK_SET_ERR( rval, "Error in get_element_connect" );
00088     int nedges = (int)edges.size();
00089 
00090     int indexConn = 0;
00091     for( int j = 0; j < nedges; j++ )
00092     {
00093         EntityHandle edge = edges[j];
00094 
00095         const EntityHandle* conn2 = NULL;
00096         int num_nodes;
00097         rval = mb->get_connectivity( edge, conn2, num_nodes );MB_CHK_ERR( rval );
00098         if( 2 != num_nodes ) MB_CHK_ERR( MB_FAILURE );
00099 
00100         int i0 = verts.index( conn2[0] );
00101         int i1 = verts.index( conn2[1] );
00102         for( int ii = 0; ii < layers; ii++ )
00103         {
00104             connect[indexConn++] = newVerts[ii][i0];
00105             connect[indexConn++] = newVerts[ii][i1];
00106             connect[indexConn++] = newVerts[ii + 1][i1];
00107             connect[indexConn++] = newVerts[ii + 1][i0];
00108         }
00109     }
00110 
00111     int nfaces                = (int)faces.size();
00112     EntityHandle* allPolygons = new EntityHandle[nfaces * ( layers + 1 )];
00113     for( int i = 0; i < nfaces; i++ )
00114         allPolygons[i] = faces[i];
00115 
00116     // vertices are parallel to the base vertices
00117     int indexVerts[MAXEDGES]                  = { 0 };  // polygons with at most MAXEDGES edges
00118     EntityHandle newConn[MAXEDGES]            = { 0 };
00119     EntityHandle polyhedronConn[MAXEDGES + 2] = { 0 };
00120 
00121     // edges will be used to determine the lateral faces of polyhedra (prisms)
00122     int indexEdges[MAXEDGES] = { 0 };  // index of edges in base polygon
00123     for( int j = 0; j < nfaces; j++ )
00124     {
00125         EntityHandle polyg = faces[j];
00126 
00127         const EntityHandle* connp = NULL;
00128         int num_nodes;
00129         rval = mb->get_connectivity( polyg, connp, num_nodes );MB_CHK_ERR( rval );
00130 
00131         for( int i = 0; i < num_nodes; i++ )
00132         {
00133             indexVerts[i]             = verts.index( connp[i] );
00134             int i1                    = ( i + 1 ) % num_nodes;
00135             EntityHandle edgeVerts[2] = { connp[i], connp[i1] };
00136             // get edge adjacent to these vertices
00137             Range adjEdges;
00138             rval = mb->get_adjacencies( edgeVerts, 2, 1, false, adjEdges );MB_CHK_ERR( rval );
00139             if( adjEdges.size() < 1 ) MB_CHK_SET_ERR( MB_FAILURE, " did not find edge " );
00140             indexEdges[i] = edges.index( adjEdges[0] );
00141             if( indexEdges[i] < 0 ) MB_CHK_SET_ERR( MB_FAILURE, "did not find edge in range" );
00142         }
00143 
00144         for( int ii = 0; ii < layers; ii++ )
00145         {
00146             // create a polygon on each layer
00147             for( int i = 0; i < num_nodes; i++ )
00148                 newConn[i] = newVerts[ii + 1][indexVerts[i]];  // vertices in layer ii+1
00149 
00150             rval = mb->create_element( MBPOLYGON, newConn, num_nodes, allPolygons[nfaces * ( ii + 1 ) + j] );MB_CHK_ERR( rval );
00151 
00152             // now create a polyhedra with top, bottom and lateral swept faces
00153             // first face is the bottom
00154             polyhedronConn[0] = allPolygons[nfaces * ii + j];
00155             // second face is the top
00156             polyhedronConn[1] = allPolygons[nfaces * ( ii + 1 ) + j];
00157             // next add lateral quads, in order of edges, using the start_elem
00158             // first layer of quads has EntityHandle from start_elem to start_elem+nedges-1
00159             // second layer of quads has EntityHandle from start_elem + nedges to
00160             // start_elem+2*nedges-1 ,etc
00161             for( int i = 0; i < num_nodes; i++ )
00162             {
00163                 polyhedronConn[2 + i] = start_elem + ii + layers * indexEdges[i];
00164             }
00165             // Create polyhedron
00166             EntityHandle polyhedron;
00167             rval = mb->create_element( MBPOLYHEDRON, polyhedronConn, 2 + num_nodes, polyhedron );MB_CHK_ERR( rval );
00168         }
00169     }
00170     rval = mb->write_file( output.c_str() );MB_CHK_ERR( rval );
00171 
00172     delete mb;
00173 
00174     return 0;
00175 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines