Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
ExtrudePoly.cpp
Go to the documentation of this file.
1 /** @example ExtrudePoly.cpp
2  * Description: read a 2d mesh in plane, extrude to form layers of prisms (polyhedra) \n
3  *
4  * To run: ./ExtrudePoly [meshfile] [outfile] [nlayers] [thickness_layer] \n
5  */
6 
7 #include "moab/Core.hpp"
8 #include "moab/ReadUtilIface.hpp"
9 #include <iostream>
10 
11 using namespace moab;
12 using namespace std;
13 
14 #ifndef MESH_DIR
15 #define MESH_DIR "."
16 #endif
17 
18 // at most 20 edges per polygon
19 #define MAXEDGES 20
20 // Note: change the file name below to test a trivial "No such file or directory" error
21 string test_file_name = string( MESH_DIR ) + string( "/io/poly8-10.vtk" );
22 string output = string( "polyhedra.vtk" );
23 int layers = 1;
24 double layer_thick = 1.0;
25 
26 int main( int argc, char** argv )
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 }