MOAB: Mesh Oriented datABase  (version 5.4.1)
ComputeTriDual.cpp
Go to the documentation of this file.
00001 /** @example ComputeTriDual.cpp \n
00002  * \brief Compute the polygonal (dual) mesh of the primal simplex grid in 2-D. \n
00003  * <b>To run</b>: ComputeTriDual  input_file output_file \n
00004  * An example to demonstrate computing the dual polygonal grid of a triangulation
00005  * on a spherical surface.
00006  */
00007 
00008 #include <iostream>
00009 #include <vector>
00010 
00011 // Include header for MOAB instance and tag conventions for
00012 #include "moab/Interface.hpp"
00013 #include "moab/Range.hpp"
00014 #include "moab/ProgOptions.hpp"
00015 #include "moab/CartVect.hpp"
00016 #include "moab/Core.hpp"
00017 #include "moab/ReadUtilIface.hpp"
00018 
00019 #ifdef MOAB_HAVE_MPI
00020 #include "moab/ParallelComm.hpp"
00021 #endif
00022 
00023 ///////////////////////////////////////////////////////////////////////////////
00024 
00025 moab::ErrorCode compute_dual_mesh( moab::Interface* mb, moab::EntityHandle& dual_set, moab::Range& cells )
00026 {
00027     moab::ErrorCode rval;
00028 
00029     moab::Range verts;
00030     rval = mb->get_connectivity( cells, verts );MB_CHK_SET_ERR( rval, "Failed to get connectivity" );
00031 
00032     moab::ReadUtilIface* iface;
00033     rval = mb->query_interface( iface );MB_CHK_SET_ERR( rval, "Can't get reader interface" );
00034 
00035     moab::EntityHandle startv;
00036     std::vector< double* > ccenters;
00037     rval = iface->get_node_coords( 3, cells.size(), 0, startv, ccenters );MB_CHK_SET_ERR( rval, "Can't get node coords" );
00038 
00039     moab::Tag gidTag = mb->globalId_tag();
00040 
00041     std::vector< int > gids( cells.size() );
00042     rval = mb->get_coords( cells, ccenters[0], ccenters[1], ccenters[2] );MB_CHK_SET_ERR( rval, "Failed to get coordinates" );
00043     rval = mb->tag_get_data( gidTag, cells, &gids[0] );MB_CHK_SET_ERR( rval, "Can't set global_id tag" );
00044 
00045     moab::Range dualverts( startv, startv + cells.size() - 1 );
00046     rval = mb->add_entities( dual_set, dualverts );MB_CHK_SET_ERR( rval, "Can't add entities" );
00047     rval = mb->tag_set_data( gidTag, dualverts, &gids[0] );MB_CHK_SET_ERR( rval, "Can't set global_id tag" );
00048 
00049 #define CC( ind ) moab::CartVect( ccenters[0][ind], ccenters[1][ind], ccenters[2][ind] )
00050 #define CCXMY( ind, cv ) \
00051     moab::CartVect( ccenters[0][ind] - ( cv )[0], ccenters[1][ind] - ( cv )[1], ccenters[2][ind] - ( cv )[2] )
00052 
00053     const moab::EntityHandle svtx = *( cells.begin() );
00054 
00055     // maintain block-wise polygon elements
00056     moab::Range dualcells;
00057     std::map< size_t, std::vector< moab::EntityHandle > > polygonConn;
00058     // Generate new Face array
00059     for( moab::Range::iterator vit = verts.begin(); vit != verts.end(); ++vit )
00060     {
00061         moab::EntityHandle vtx = *vit;
00062 
00063         std::vector< moab::EntityHandle > adjs;
00064         // generate all required adjacencies
00065         rval = mb->get_adjacencies( &vtx, 1, 2, true, adjs );MB_CHK_SET_ERR( rval, "Failed to get adjacencies" );
00066         const size_t nEdges = adjs.size();
00067 
00068         double vxyz[3];
00069         rval = mb->get_coords( &vtx, 1, vxyz );MB_CHK_SET_ERR( rval, "Failed to get coordinates" );
00070 
00071         // Reorient Faces
00072         moab::CartVect nodeCentral( vxyz );
00073         moab::CartVect noderef = CC( adjs[0] - svtx );
00074         moab::CartVect node0   = noderef - nodeCentral;
00075         double dNode0Mag       = node0.length();
00076 
00077         moab::CartVect nodeCross = noderef * nodeCentral;
00078 
00079         // Determine the angles about the central Node of each Face Node
00080         std::vector< double > dAngles;
00081         dAngles.resize( nEdges );
00082         dAngles[0] = 0.0;
00083 
00084         for( size_t j = 1; j < nEdges; j++ )
00085         {
00086             moab::CartVect nodeDiff = CCXMY( adjs[j] - svtx, nodeCentral );
00087             double dNodeDiffMag     = nodeDiff.length();
00088 
00089             double dSide    = nodeCross % nodeDiff;
00090             double dDotNorm = ( node0 % nodeDiff ) / ( dNode0Mag * dNodeDiffMag );
00091 
00092             // double dAngle;
00093             if( dDotNorm > 1.0 )
00094             {
00095                 dDotNorm = 1.0;
00096             }
00097 
00098             dAngles[j] = acos( dDotNorm );
00099 
00100             if( dSide > 0.0 )
00101             {
00102                 dAngles[j] = -dAngles[j] + 2.0 * M_PI;
00103             }
00104         }
00105 
00106         std::vector< moab::EntityHandle >& face = polygonConn[nEdges];
00107         if( face.size() == 0 ) face.reserve( nEdges * verts.size() / 4 );
00108         face.push_back( dualverts[adjs[0] - svtx] );
00109 
00110         // Orient each Face by putting Nodes in order of increasing angle
00111         double dCurrentAngle = 0.0;
00112         for( size_t j = 1; j < nEdges; j++ )
00113         {
00114             int ixNextNode    = 1;
00115             double dNextAngle = 2.0 * M_PI;
00116 
00117             for( size_t k = 1; k < nEdges; k++ )
00118             {
00119                 if( ( dAngles[k] > dCurrentAngle ) && ( dAngles[k] < dNextAngle ) )
00120                 {
00121                     ixNextNode = k;
00122                     dNextAngle = dAngles[k];
00123                 }
00124             }
00125 
00126             face.push_back( dualverts[adjs[ixNextNode] - svtx] );
00127             dCurrentAngle = dNextAngle;
00128         }
00129     }
00130 
00131     std::map< size_t, std::vector< moab::EntityHandle > >::iterator eit;
00132     for( eit = polygonConn.begin(); eit != polygonConn.end(); ++eit )
00133     {
00134         const size_t nVPerE      = eit->first;
00135         const size_t nElePerType = static_cast< int >( eit->second.size() / nVPerE );
00136 
00137         moab::EntityHandle starte;  // Connectivity
00138         moab::EntityHandle* conn;
00139 
00140         rval = iface->get_element_connect( nElePerType, eit->first, moab::MBPOLYGON, 0, starte, conn );MB_CHK_SET_ERR( rval, "Can't get element connectivity" );
00141 
00142         // copy the connectivity that we have accumulated
00143         std::copy( eit->second.begin(), eit->second.end(), conn );
00144         eit->second.clear();
00145 
00146         // add this polygon sequence to the aggregate data
00147         moab::Range mbcells( starte, starte + nElePerType - 1 );
00148         dualcells.merge( mbcells );
00149     }
00150 
00151     // add the computed dual cells to mesh
00152     rval = mb->add_entities( dual_set, dualcells );MB_CHK_SET_ERR( rval, "Can't add polygonal entities" );
00153 
00154     // Assign global IDs to all the dual cells - same as original vertices
00155     assert( dualcells.size() == verts.size() );
00156     gids.resize( verts.size() );
00157     rval = mb->tag_get_data( gidTag, verts, &gids[0] );MB_CHK_SET_ERR( rval, "Can't set global_id tag" );
00158     if( gids[0] == gids[1] && gids[0] < 0 )
00159     {
00160 #ifdef MOAB_HAVE_MPI
00161         moab::ParallelComm pcomm( mb, MPI_COMM_WORLD );
00162 
00163         rval = pcomm.assign_global_ids( dual_set, 2, 1, false, true, true );MB_CHK_SET_ERR( rval, "Can't assign global_ids" );
00164 #else
00165         // No global ID assigned to input vertices.
00166         // Can we use std::iota ??
00167         for( size_t ix = 0; ix < gids.size(); ++ix )
00168             gids[ix] = ix + 1;
00169 
00170         // set GID tag
00171         rval = mb->tag_set_data( gidTag, dualcells, &gids[0] );MB_CHK_SET_ERR( rval, "Can't set global_id tag" );
00172         std::cout << "GIDs: " << gids[0] << " " << gids[1] << " " << gids[2] << " " << gids[3] << " " << gids[4] << " "
00173                   << gids[5] << "\n";
00174 #endif
00175     }
00176     gids.clear();
00177 
00178     // delete the original entities from the mesh
00179     rval = mb->delete_entities( cells );MB_CHK_SET_ERR( rval, "Can't remove entities" );
00180     rval = mb->delete_entities( verts );MB_CHK_SET_ERR( rval, "Can't remove entities" );
00181 
00182     return moab::MB_SUCCESS;
00183 }
00184 
00185 ///////////////////////////////////////////////////////////////////////////////
00186 
00187 int main( int argc, char** argv )
00188 {
00189     moab::ErrorCode rval;
00190     ProgOptions opts;
00191 
00192     // Get MOAB instance
00193     moab::Interface* mb = new( std::nothrow ) moab::Core;
00194     if( NULL == mb ) return 1;
00195 
00196     std::string inputFile = std::string( MESH_DIR ) + "/globalmpas_deltri.h5m";
00197     opts.addOpt< std::string >( "inFile,i", "Specify the input file name string ", &inputFile );
00198 #ifdef MOAB_HAVE_HDF5
00199     std::string outputFile = "outDualMesh.h5m";
00200 #else
00201     std::string outputFile = "outDualMesh.vtk";
00202 #endif
00203     opts.addOpt< std::string >( "outFile,o", "Specify the input file name string ", &outputFile );
00204 
00205     opts.parseCommandLine( argc, argv );
00206 
00207     moab::EntityHandle triangle_set, dual_set;
00208     rval = mb->create_meshset( moab::MESHSET_SET, triangle_set );MB_CHK_SET_ERR( rval, "Can't create new set" );
00209     rval = mb->create_meshset( moab::MESHSET_SET, dual_set );MB_CHK_SET_ERR( rval, "Can't create new set" );
00210 
00211     // This file is in the mesh files directory
00212     const char* readopts = "";
00213     rval                 = mb->load_file( inputFile.c_str(), &triangle_set, readopts );MB_CHK_SET_ERR( rval, "Failed to read" );
00214 
00215     // get all cells of dimension 2;
00216     moab::Range cells;
00217     rval = mb->get_entities_by_dimension( triangle_set, 2, cells );MB_CHK_SET_ERR( rval, "Failed to get cells" );
00218 
00219     std::cout << "Original number of triangular cells : " << cells.size() << "\n";
00220 
00221     // call the routine to compute the dual grid
00222     rval = compute_dual_mesh( mb, dual_set, cells );MB_CHK_SET_ERR( rval, "Failed to compute dual mesh" );
00223 
00224     // write the mesh to disk
00225     rval = mb->write_file( outputFile.c_str() );MB_CHK_SET_ERR( rval, "Failed to write new file" );
00226 
00227     cells.clear();
00228     rval = mb->get_entities_by_dimension( dual_set, 2, cells );MB_CHK_SET_ERR( rval, "Failed to get cells" );
00229     std::cout << "Wrote the dual mesh: " << outputFile << " containing " << cells.size() << " polygonal cells\n";
00230 
00231     delete mb;
00232 
00233     return 0;
00234 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines