 |
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
|
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
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
25 moab::ErrorCode compute_dual_mesh( moab::Interface* mb, moab::EntityHandle& dual_set, moab::Range& cells )
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
32 moab::ReadUtilIface* iface;
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 }
1.9.1.