MOAB: Mesh Oriented datABase
(version 5.4.1)
|
Compute the polygonal (dual) mesh of the primal simplex grid in 2-D.
To run: ComputeTriDual input_file output_file
An example to demonstrate computing the dual polygonal grid of a triangulation on a spherical surface.
/** @example ComputeTriDual.cpp \n * \brief Compute the polygonal (dual) mesh of the primal simplex grid in 2-D. \n * <b>To run</b>: ComputeTriDual input_file output_file \n * An example to demonstrate computing the dual polygonal grid of a triangulation * on a spherical surface. */ #include <iostream> #include <vector> // Include header for MOAB instance and tag conventions for #include "moab/Interface.hpp" #include "moab/Range.hpp" #include "moab/ProgOptions.hpp" #include "moab/CartVect.hpp" #include "moab/Core.hpp" #include "moab/ReadUtilIface.hpp" #ifdef MOAB_HAVE_MPI #include "moab/ParallelComm.hpp" #endif /////////////////////////////////////////////////////////////////////////////// moab::ErrorCode compute_dual_mesh( moab::Interface* mb, moab::EntityHandle& dual_set, moab::Range& cells ) { moab::ErrorCode rval; moab::Range verts; rval = mb->get_connectivity( cells, verts );MB_CHK_SET_ERR( rval, "Failed to get connectivity" ); moab::ReadUtilIface* iface; rval = mb->query_interface( iface );MB_CHK_SET_ERR( rval, "Can't get reader interface" ); moab::EntityHandle startv; std::vector< double* > ccenters; rval = iface->get_node_coords( 3, cells.size(), 0, startv, ccenters );MB_CHK_SET_ERR( rval, "Can't get node coords" ); moab::Tag gidTag = mb->globalId_tag(); std::vector< int > gids( cells.size() ); rval = mb->get_coords( cells, ccenters[0], ccenters[1], ccenters[2] );MB_CHK_SET_ERR( rval, "Failed to get coordinates" ); rval = mb->tag_get_data( gidTag, cells, &gids[0] );MB_CHK_SET_ERR( rval, "Can't set global_id tag" ); moab::Range dualverts( startv, startv + cells.size() - 1 ); rval = mb->add_entities( dual_set, dualverts );MB_CHK_SET_ERR( rval, "Can't add entities" ); rval = mb->tag_set_data( gidTag, dualverts, &gids[0] );MB_CHK_SET_ERR( rval, "Can't set global_id tag" ); #define CC( ind ) moab::CartVect( ccenters[0][ind], ccenters[1][ind], ccenters[2][ind] ) #define CCXMY( ind, cv ) \ moab::CartVect( ccenters[0][ind] - ( cv )[0], ccenters[1][ind] - ( cv )[1], ccenters[2][ind] - ( cv )[2] ) const moab::EntityHandle svtx = *( cells.begin() ); // maintain block-wise polygon elements moab::Range dualcells; std::map< size_t, std::vector< moab::EntityHandle > > polygonConn; // Generate new Face array for( moab::Range::iterator vit = verts.begin(); vit != verts.end(); ++vit ) { moab::EntityHandle vtx = *vit; std::vector< moab::EntityHandle > adjs; // generate all required adjacencies rval = mb->get_adjacencies( &vtx, 1, 2, true, adjs );MB_CHK_SET_ERR( rval, "Failed to get adjacencies" ); const size_t nEdges = adjs.size(); double vxyz[3]; rval = mb->get_coords( &vtx, 1, vxyz );MB_CHK_SET_ERR( rval, "Failed to get coordinates" ); // Reorient Faces moab::CartVect nodeCentral( vxyz ); moab::CartVect noderef = CC( adjs[0] - svtx ); moab::CartVect node0 = noderef - nodeCentral; double dNode0Mag = node0.length(); moab::CartVect nodeCross = noderef * nodeCentral; // Determine the angles about the central Node of each Face Node std::vector< double > dAngles; dAngles.resize( nEdges ); dAngles[0] = 0.0; for( size_t j = 1; j < nEdges; j++ ) { moab::CartVect nodeDiff = CCXMY( adjs[j] - svtx, nodeCentral ); double dNodeDiffMag = nodeDiff.length(); double dSide = nodeCross % nodeDiff; double dDotNorm = ( node0 % nodeDiff ) / ( dNode0Mag * dNodeDiffMag ); // double dAngle; if( dDotNorm > 1.0 ) { dDotNorm = 1.0; } dAngles[j] = acos( dDotNorm ); if( dSide > 0.0 ) { dAngles[j] = -dAngles[j] + 2.0 * M_PI; } } std::vector< moab::EntityHandle >& face = polygonConn[nEdges]; if( face.size() == 0 ) face.reserve( nEdges * verts.size() / 4 ); face.push_back( dualverts[adjs[0] - svtx] ); // Orient each Face by putting Nodes in order of increasing angle double dCurrentAngle = 0.0; for( size_t j = 1; j < nEdges; j++ ) { int ixNextNode = 1; double dNextAngle = 2.0 * M_PI; for( size_t k = 1; k < nEdges; k++ ) { if( ( dAngles[k] > dCurrentAngle ) && ( dAngles[k] < dNextAngle ) ) { ixNextNode = k; dNextAngle = dAngles[k]; } } face.push_back( dualverts[adjs[ixNextNode] - svtx] ); dCurrentAngle = dNextAngle; } } std::map< size_t, std::vector< moab::EntityHandle > >::iterator eit; for( eit = polygonConn.begin(); eit != polygonConn.end(); ++eit ) { const size_t nVPerE = eit->first; const size_t nElePerType = static_cast< int >( eit->second.size() / nVPerE ); moab::EntityHandle starte; // Connectivity moab::EntityHandle* conn; rval = iface->get_element_connect( nElePerType, eit->first, moab::MBPOLYGON, 0, starte, conn );MB_CHK_SET_ERR( rval, "Can't get element connectivity" ); // copy the connectivity that we have accumulated std::copy( eit->second.begin(), eit->second.end(), conn ); eit->second.clear(); // add this polygon sequence to the aggregate data moab::Range mbcells( starte, starte + nElePerType - 1 ); dualcells.merge( mbcells ); } // add the computed dual cells to mesh rval = mb->add_entities( dual_set, dualcells );MB_CHK_SET_ERR( rval, "Can't add polygonal entities" ); // Assign global IDs to all the dual cells - same as original vertices assert( dualcells.size() == verts.size() ); gids.resize( verts.size() ); rval = mb->tag_get_data( gidTag, verts, &gids[0] );MB_CHK_SET_ERR( rval, "Can't set global_id tag" ); if( gids[0] == gids[1] && gids[0] < 0 ) { #ifdef MOAB_HAVE_MPI moab::ParallelComm pcomm( mb, MPI_COMM_WORLD ); rval = pcomm.assign_global_ids( dual_set, 2, 1, false, true, true );MB_CHK_SET_ERR( rval, "Can't assign global_ids" ); #else // No global ID assigned to input vertices. // Can we use std::iota ?? for( size_t ix = 0; ix < gids.size(); ++ix ) gids[ix] = ix + 1; // set GID tag rval = mb->tag_set_data( gidTag, dualcells, &gids[0] );MB_CHK_SET_ERR( rval, "Can't set global_id tag" ); std::cout << "GIDs: " << gids[0] << " " << gids[1] << " " << gids[2] << " " << gids[3] << " " << gids[4] << " " << gids[5] << "\n"; #endif } gids.clear(); // delete the original entities from the mesh rval = mb->delete_entities( cells );MB_CHK_SET_ERR( rval, "Can't remove entities" ); rval = mb->delete_entities( verts );MB_CHK_SET_ERR( rval, "Can't remove entities" ); return moab::MB_SUCCESS; } /////////////////////////////////////////////////////////////////////////////// int main( int argc, char** argv ) { moab::ErrorCode rval; ProgOptions opts; // Get MOAB instance moab::Interface* mb = new( std::nothrow ) moab::Core; if( NULL == mb ) return 1; std::string inputFile = std::string( MESH_DIR ) + "/globalmpas_deltri.h5m"; opts.addOpt< std::string >( "inFile,i", "Specify the input file name string ", &inputFile ); #ifdef MOAB_HAVE_HDF5 std::string outputFile = "outDualMesh.h5m"; #else std::string outputFile = "outDualMesh.vtk"; #endif opts.addOpt< std::string >( "outFile,o", "Specify the input file name string ", &outputFile ); opts.parseCommandLine( argc, argv ); moab::EntityHandle triangle_set, dual_set; rval = mb->create_meshset( moab::MESHSET_SET, triangle_set );MB_CHK_SET_ERR( rval, "Can't create new set" ); rval = mb->create_meshset( moab::MESHSET_SET, dual_set );MB_CHK_SET_ERR( rval, "Can't create new set" ); // This file is in the mesh files directory const char* readopts = ""; rval = mb->load_file( inputFile.c_str(), &triangle_set, readopts );MB_CHK_SET_ERR( rval, "Failed to read" ); // get all cells of dimension 2; moab::Range cells; rval = mb->get_entities_by_dimension( triangle_set, 2, cells );MB_CHK_SET_ERR( rval, "Failed to get cells" ); std::cout << "Original number of triangular cells : " << cells.size() << "\n"; // call the routine to compute the dual grid rval = compute_dual_mesh( mb, dual_set, cells );MB_CHK_SET_ERR( rval, "Failed to compute dual mesh" ); // write the mesh to disk rval = mb->write_file( outputFile.c_str() );MB_CHK_SET_ERR( rval, "Failed to write new file" ); cells.clear(); rval = mb->get_entities_by_dimension( dual_set, 2, cells );MB_CHK_SET_ERR( rval, "Failed to get cells" ); std::cout << "Wrote the dual mesh: " << outputFile << " containing " << cells.size() << " polygonal cells\n"; delete mb; return 0; }