#include "moab/Core.hpp"
#include "moab/ReadUtilIface.hpp"
#include <iostream>

Include dependency graph for ExtrudePoly.cpp:

Macros
#define	MESH_DIR "."

#define	MAXEDGES 20

Functions
int	main (int argc, char **argv)

Variables
string	test_file_name = string( MESH_DIR ) + string( "/io/poly8-10.vtk" )

string	output = string( "polyhedra.vtk" )

int	layers = 1

double	layer_thick = 1.0

Macro Definition Documentation

◆ MAXEDGES

#define MAXEDGES 20

Examples: ExtrudePoly.cpp.

Definition at line 19 of file ExtrudePoly.cpp.

◆ MESH_DIR

#define MESH_DIR "."

Examples: ComputeTriDual.cpp, ExtrudePoly.cpp, HelloParMOAB.cpp, LaplacianSmoother.cpp, LloydRelaxation.cpp, LoadPartial.cpp, ReadWriteTest.cpp, ReduceExchangeTags.cpp, TestErrorHandling.cpp, TestErrorHandlingPar.cpp, TestExodusII.cpp, and VisTags.cpp.

Definition at line 15 of file ExtrudePoly.cpp.

Function Documentation

◆ main()

int main	(	int	argc,
		char **	argv
	)

Examples: ExtrudePoly.cpp.

Definition at line 26 of file ExtrudePoly.cpp.

   27 {
   28     // Get MOAB instance
   29     Interface* mb = new( std::nothrow ) Core;
   30     if( NULL == mb ) return 1;
   31  
   32     // Need option handling here for input filename
   33     if( argc > 1 )
   34     {
   35         // User has input a mesh file
   36         test_file_name = argv[1];
   37     }
   38     if( argc > 2 ) output = argv[2];
   39  
   40     if( argc > 3 ) layers = atoi( argv[3] );
   41  
   42     if( argc > 4 ) layer_thick = atof( argv[4] );
   43  
   44     std::cout << "Run: " << argv[0] << " " << test_file_name << " " << output << " " << layers << " " << layer_thick
   45               << "\n";
   46     // Load the mesh from vtk file
   47     ErrorCode rval = mb->load_mesh( test_file_name.c_str() );MB_CHK_ERR( rval );
   48  
   49     // Get verts entities, by type
   50     Range verts;
   51     rval = mb->get_entities_by_type( 0, MBVERTEX, verts );MB_CHK_ERR( rval );
   52  
   53     // Get faces, by dimension, so we stay generic to entity type
   54     Range faces;
   55     rval = mb->get_entities_by_dimension( 0, 2, faces );MB_CHK_ERR( rval );
   56     cout << "Number of vertices is " << verts.size() << endl;
   57     cout << "Number of faces is " << faces.size() << endl;
   58  
   59     Range edges;
   60     // create all edges
   61     rval = mb->get_adjacencies( faces, 1, true, edges, Interface::UNION );MB_CHK_ERR( rval );
   62  
   63     cout << "Number of edges is " << edges.size() << endl;
   64  
   65     std::vector< double > coords;
   66     int nvPerLayer = (int)verts.size();
   67     coords.resize( 3 * nvPerLayer );
   68  
   69     rval = mb->get_coords( verts, &coords[0] );MB_CHK_ERR( rval );
   70     // create first vertices
   71     Range* newVerts = new Range[layers + 1];
   72     newVerts[0]     = verts;  // just for convenience
   73     for( int ii = 0; ii < layers; ii++ )
   74     {
   75         for( int i = 0; i < nvPerLayer; i++ )
   76             coords[3 * i + 2] += layer_thick;
   77  
   78         rval = mb->create_vertices( &coords[0], nvPerLayer, newVerts[ii + 1] );MB_CHK_ERR( rval );
   79     }
   80     // for each edge, we will create layers quads
   81     int nquads = edges.size() * layers;
   82     ReadUtilIface* read_iface;
   83     rval = mb->query_interface( read_iface );MB_CHK_SET_ERR( rval, "Error in query_interface" );
   84  
   85     EntityHandle start_elem, *connect;
   86     // Create quads
   87     rval = read_iface->get_element_connect( nquads, 4, MBQUAD, 0, start_elem, connect );MB_CHK_SET_ERR( rval, "Error in get_element_connect" );
   88     int nedges = (int)edges.size();
   89  
   90     int indexConn = 0;
   91     for( int j = 0; j < nedges; j++ )
   92     {
   93         EntityHandle edge = edges[j];
   94  
   95         const EntityHandle* conn2 = NULL;
   96         int num_nodes;
   97         rval = mb->get_connectivity( edge, conn2, num_nodes );MB_CHK_ERR( rval );
   98         if( 2 != num_nodes ) MB_CHK_ERR( MB_FAILURE );
   99  
  100         int i0 = verts.index( conn2[0] );
  101         int i1 = verts.index( conn2[1] );
  102         for( int ii = 0; ii < layers; ii++ )
  103         {
  104             connect[indexConn++] = newVerts[ii][i0];
  105             connect[indexConn++] = newVerts[ii][i1];
  106             connect[indexConn++] = newVerts[ii + 1][i1];
  107             connect[indexConn++] = newVerts[ii + 1][i0];
  108         }
  109     }
  110  
  111     int nfaces                = (int)faces.size();
  112     EntityHandle* allPolygons = new EntityHandle[nfaces * ( layers + 1 )];
  113     for( int i = 0; i < nfaces; i++ )
  114         allPolygons[i] = faces[i];
  115  
  116     // vertices are parallel to the base vertices
  117     int indexVerts[MAXEDGES]                  = { 0 };  // polygons with at most MAXEDGES edges
  118     EntityHandle newConn[MAXEDGES]            = { 0 };
  119     EntityHandle polyhedronConn[MAXEDGES + 2] = { 0 };
  120  
  121     // edges will be used to determine the lateral faces of polyhedra (prisms)
  122     int indexEdges[MAXEDGES] = { 0 };  // index of edges in base polygon
  123     for( int j = 0; j < nfaces; j++ )
  124     {
  125         EntityHandle polyg = faces[j];
  126  
  127         const EntityHandle* connp = NULL;
  128         int num_nodes;
  129         rval = mb->get_connectivity( polyg, connp, num_nodes );MB_CHK_ERR( rval );
  130  
  131         for( int i = 0; i < num_nodes; i++ )
  132         {
  133             indexVerts[i]             = verts.index( connp[i] );
  134             int i1                    = ( i + 1 ) % num_nodes;
  135             EntityHandle edgeVerts[2] = { connp[i], connp[i1] };
  136             // get edge adjacent to these vertices
  137             Range adjEdges;
  138             rval = mb->get_adjacencies( edgeVerts, 2, 1, false, adjEdges );MB_CHK_ERR( rval );
  139             if( adjEdges.size() < 1 ) MB_CHK_SET_ERR( MB_FAILURE, " did not find edge " );
  140             indexEdges[i] = edges.index( adjEdges[0] );
  141             if( indexEdges[i] < 0 ) MB_CHK_SET_ERR( MB_FAILURE, "did not find edge in range" );
  142         }
  143  
  144         for( int ii = 0; ii < layers; ii++ )
  145         {
  146             // create a polygon on each layer
  147             for( int i = 0; i < num_nodes; i++ )
  148                 newConn[i] = newVerts[ii + 1][indexVerts[i]];  // vertices in layer ii+1
  149  
  150             rval = mb->create_element( MBPOLYGON, newConn, num_nodes, allPolygons[nfaces * ( ii + 1 ) + j] );MB_CHK_ERR( rval );
  151  
  152             // now create a polyhedra with top, bottom and lateral swept faces
  153             // first face is the bottom
  154             polyhedronConn[0] = allPolygons[nfaces * ii + j];
  155             // second face is the top
  156             polyhedronConn[1] = allPolygons[nfaces * ( ii + 1 ) + j];
  157             // next add lateral quads, in order of edges, using the start_elem
  158             // first layer of quads has EntityHandle from start_elem to start_elem+nedges-1
  159             // second layer of quads has EntityHandle from start_elem + nedges to
  160             // start_elem+2*nedges-1 ,etc
  161             for( int i = 0; i < num_nodes; i++ )
  162             {
  163                 polyhedronConn[2 + i] = start_elem + ii + layers * indexEdges[i];
  164             }
  165             // Create polyhedron
  166             EntityHandle polyhedron;
  167             rval = mb->create_element( MBPOLYHEDRON, polyhedronConn, 2 + num_nodes, polyhedron );MB_CHK_ERR( rval );
  168         }
  169     }
  170     rval = mb->write_file( output.c_str() );MB_CHK_ERR( rval );
  171  
  172     delete mb;
  173  
  174     return 0;
  175 }