Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
ComputeTriDual.cpp
Go to the documentation of this file.
1 /** @example ComputeTriDual.cpp \n
2  * \brief Compute the polygonal (dual) mesh of the primal simplex grid in 2-D. \n
3  * <b>To run</b>: ComputeTriDual input_file output_file \n
4  * An example to demonstrate computing the dual polygonal grid of a triangulation
5  * on a spherical surface.
6  */
7 
8 #include <iostream>
9 #include <vector>
10 
11 // Include header for MOAB instance and tag conventions for
12 #include "moab/Interface.hpp"
13 #include "moab/Range.hpp"
14 #include "moab/ProgOptions.hpp"
15 #include "moab/CartVect.hpp"
16 #include "moab/Core.hpp"
17 #include "moab/ReadUtilIface.hpp"
18 
19 #ifdef MOAB_HAVE_MPI
20 #include "moab/ParallelComm.hpp"
21 #endif
22 
23 ///////////////////////////////////////////////////////////////////////////////
24 
26 {
27  moab::ErrorCode rval;
28 
29  moab::Range verts;
30  rval = mb->get_connectivity( cells, verts );MB_CHK_SET_ERR( rval, "Failed to get connectivity" );
31 
33  rval = mb->query_interface( iface );MB_CHK_SET_ERR( rval, "Can't get reader interface" );
34 
35  moab::EntityHandle startv;
36  std::vector< double* > ccenters;
37  rval = iface->get_node_coords( 3, cells.size(), 0, startv, ccenters );MB_CHK_SET_ERR( rval, "Can't get node coords" );
38 
39  moab::Tag gidTag = mb->globalId_tag();
40 
41  std::vector< int > gids( cells.size() );
42  rval = mb->get_coords( cells, ccenters[0], ccenters[1], ccenters[2] );MB_CHK_SET_ERR( rval, "Failed to get coordinates" );
43  rval = mb->tag_get_data( gidTag, cells, &gids[0] );MB_CHK_SET_ERR( rval, "Can't set global_id tag" );
44 
45  moab::Range dualverts( startv, startv + cells.size() - 1 );
46  rval = mb->add_entities( dual_set, dualverts );MB_CHK_SET_ERR( rval, "Can't add entities" );
47  rval = mb->tag_set_data( gidTag, dualverts, &gids[0] );MB_CHK_SET_ERR( rval, "Can't set global_id tag" );
48 
49 #define CC( ind ) moab::CartVect( ccenters[0][ind], ccenters[1][ind], ccenters[2][ind] )
50 #define CCXMY( ind, cv ) \
51  moab::CartVect( ccenters[0][ind] - ( cv )[0], ccenters[1][ind] - ( cv )[1], ccenters[2][ind] - ( cv )[2] )
52 
53  const moab::EntityHandle svtx = *( cells.begin() );
54 
55  // maintain block-wise polygon elements
56  moab::Range dualcells;
57  std::map< size_t, std::vector< moab::EntityHandle > > polygonConn;
58  // Generate new Face array
59  for( moab::Range::iterator vit = verts.begin(); vit != verts.end(); ++vit )
60  {
61  moab::EntityHandle vtx = *vit;
62 
63  std::vector< moab::EntityHandle > adjs;
64  // generate all required adjacencies
65  rval = mb->get_adjacencies( &vtx, 1, 2, true, adjs );MB_CHK_SET_ERR( rval, "Failed to get adjacencies" );
66  const size_t nEdges = adjs.size();
67 
68  double vxyz[3];
69  rval = mb->get_coords( &vtx, 1, vxyz );MB_CHK_SET_ERR( rval, "Failed to get coordinates" );
70 
71  // Reorient Faces
72  moab::CartVect nodeCentral( vxyz );
73  moab::CartVect noderef = CC( adjs[0] - svtx );
74  moab::CartVect node0 = noderef - nodeCentral;
75  double dNode0Mag = node0.length();
76 
77  moab::CartVect nodeCross = noderef * nodeCentral;
78 
79  // Determine the angles about the central Node of each Face Node
80  std::vector< double > dAngles;
81  dAngles.resize( nEdges );
82  dAngles[0] = 0.0;
83 
84  for( size_t j = 1; j < nEdges; j++ )
85  {
86  moab::CartVect nodeDiff = CCXMY( adjs[j] - svtx, nodeCentral );
87  double dNodeDiffMag = nodeDiff.length();
88 
89  double dSide = nodeCross % nodeDiff;
90  double dDotNorm = ( node0 % nodeDiff ) / ( dNode0Mag * dNodeDiffMag );
91 
92  // double dAngle;
93  if( dDotNorm > 1.0 )
94  {
95  dDotNorm = 1.0;
96  }
97 
98  dAngles[j] = acos( dDotNorm );
99 
100  if( dSide > 0.0 )
101  {
102  dAngles[j] = -dAngles[j] + 2.0 * M_PI;
103  }
104  }
105 
106  std::vector< moab::EntityHandle >& face = polygonConn[nEdges];
107  if( face.size() == 0 ) face.reserve( nEdges * verts.size() / 4 );
108  face.push_back( dualverts[adjs[0] - svtx] );
109 
110  // Orient each Face by putting Nodes in order of increasing angle
111  double dCurrentAngle = 0.0;
112  for( size_t j = 1; j < nEdges; j++ )
113  {
114  int ixNextNode = 1;
115  double dNextAngle = 2.0 * M_PI;
116 
117  for( size_t k = 1; k < nEdges; k++ )
118  {
119  if( ( dAngles[k] > dCurrentAngle ) && ( dAngles[k] < dNextAngle ) )
120  {
121  ixNextNode = k;
122  dNextAngle = dAngles[k];
123  }
124  }
125 
126  face.push_back( dualverts[adjs[ixNextNode] - svtx] );
127  dCurrentAngle = dNextAngle;
128  }
129  }
130 
131  std::map< size_t, std::vector< moab::EntityHandle > >::iterator eit;
132  for( eit = polygonConn.begin(); eit != polygonConn.end(); ++eit )
133  {
134  const size_t nVPerE = eit->first;
135  const size_t nElePerType = static_cast< int >( eit->second.size() / nVPerE );
136 
137  moab::EntityHandle starte; // Connectivity
138  moab::EntityHandle* conn;
139 
140  rval = iface->get_element_connect( nElePerType, eit->first, moab::MBPOLYGON, 0, starte, conn );MB_CHK_SET_ERR( rval, "Can't get element connectivity" );
141 
142  // copy the connectivity that we have accumulated
143  std::copy( eit->second.begin(), eit->second.end(), conn );
144  eit->second.clear();
145 
146  // add this polygon sequence to the aggregate data
147  moab::Range mbcells( starte, starte + nElePerType - 1 );
148  dualcells.merge( mbcells );
149  }
150 
151  // add the computed dual cells to mesh
152  rval = mb->add_entities( dual_set, dualcells );MB_CHK_SET_ERR( rval, "Can't add polygonal entities" );
153 
154  // Assign global IDs to all the dual cells - same as original vertices
155  assert( dualcells.size() == verts.size() );
156  gids.resize( verts.size() );
157  rval = mb->tag_get_data( gidTag, verts, &gids[0] );MB_CHK_SET_ERR( rval, "Can't set global_id tag" );
158  if( gids[0] == gids[1] && gids[0] < 0 )
159  {
160 #ifdef MOAB_HAVE_MPI
161  moab::ParallelComm pcomm( mb, MPI_COMM_WORLD );
162 
163  rval = pcomm.assign_global_ids( dual_set, 2, 1, false, true, true );MB_CHK_SET_ERR( rval, "Can't assign global_ids" );
164 #else
165  // No global ID assigned to input vertices.
166  // Can we use std::iota ??
167  for( size_t ix = 0; ix < gids.size(); ++ix )
168  gids[ix] = ix + 1;
169 
170  // set GID tag
171  rval = mb->tag_set_data( gidTag, dualcells, &gids[0] );MB_CHK_SET_ERR( rval, "Can't set global_id tag" );
172  std::cout << "GIDs: " << gids[0] << " " << gids[1] << " " << gids[2] << " " << gids[3] << " " << gids[4] << " "
173  << gids[5] << "\n";
174 #endif
175  }
176  gids.clear();
177 
178  // delete the original entities from the mesh
179  rval = mb->delete_entities( cells );MB_CHK_SET_ERR( rval, "Can't remove entities" );
180  rval = mb->delete_entities( verts );MB_CHK_SET_ERR( rval, "Can't remove entities" );
181 
182  return moab::MB_SUCCESS;
183 }
184 
185 ///////////////////////////////////////////////////////////////////////////////
186 
187 int main( int argc, char** argv )
188 {
189  moab::ErrorCode rval;
190  ProgOptions opts;
191 
192  // Get MOAB instance
193  moab::Interface* mb = new( std::nothrow ) moab::Core;
194  if( NULL == mb ) return 1;
195 
196  std::string inputFile = std::string( MESH_DIR ) + "/globalmpas_deltri.h5m";
197  opts.addOpt< std::string >( "inFile,i", "Specify the input file name string ", &inputFile );
198 #ifdef MOAB_HAVE_HDF5
199  std::string outputFile = "outDualMesh.h5m";
200 #else
201  std::string outputFile = "outDualMesh.vtk";
202 #endif
203  opts.addOpt< std::string >( "outFile,o", "Specify the input file name string ", &outputFile );
204 
205  opts.parseCommandLine( argc, argv );
206 
207  moab::EntityHandle triangle_set, dual_set;
208  rval = mb->create_meshset( moab::MESHSET_SET, triangle_set );MB_CHK_SET_ERR( rval, "Can't create new set" );
209  rval = mb->create_meshset( moab::MESHSET_SET, dual_set );MB_CHK_SET_ERR( rval, "Can't create new set" );
210 
211  // This file is in the mesh files directory
212  const char* readopts = "";
213  rval = mb->load_file( inputFile.c_str(), &triangle_set, readopts );MB_CHK_SET_ERR( rval, "Failed to read" );
214 
215  // get all cells of dimension 2;
216  moab::Range cells;
217  rval = mb->get_entities_by_dimension( triangle_set, 2, cells );MB_CHK_SET_ERR( rval, "Failed to get cells" );
218 
219  std::cout << "Original number of triangular cells : " << cells.size() << "\n";
220 
221  // call the routine to compute the dual grid
222  rval = compute_dual_mesh( mb, dual_set, cells );MB_CHK_SET_ERR( rval, "Failed to compute dual mesh" );
223 
224  // write the mesh to disk
225  rval = mb->write_file( outputFile.c_str() );MB_CHK_SET_ERR( rval, "Failed to write new file" );
226 
227  cells.clear();
228  rval = mb->get_entities_by_dimension( dual_set, 2, cells );MB_CHK_SET_ERR( rval, "Failed to get cells" );
229  std::cout << "Wrote the dual mesh: " << outputFile << " containing " << cells.size() << " polygonal cells\n";
230 
231  delete mb;
232 
233  return 0;
234 }