ComputeTriDual.cpp

Go to the documentation of this file.

/** @example ComputeTriDual.cpp
 * This example demonstrates computation of dual polygonal mesh from primal triangulation.
 * It shows how to load a triangular mesh from a file,
 * compute the dual polygonal mesh of a triangulation,
 * handle spherical surface triangulations,
 * create dual vertices at triangle centroids,
 * generate dual polygons around primal vertices,
 * orient dual polygons correctly on spherical surfaces,
 * preserve global IDs during dual mesh construction,
 * and write the dual mesh to a new file.
 *
 * The dual mesh computation is particularly useful for finite volume
 * methods and other applications requiring dual grid structures.
 *
 * To run: ComputeTriDual input_file output_file
 * 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;
    MB_CHK_SET_ERR( mb->get_connectivity( cells, verts ), "Failed to get connectivity" );
 
    moab::ReadUtilIface* iface;
    MB_CHK_SET_ERR( mb->query_interface( iface ), "Can't get reader interface" );
 
    moab::EntityHandle startv;
    std::vector< double* > ccenters;
    MB_CHK_SET_ERR( iface->get_node_coords( 3, cells.size(), 0, startv, ccenters ), "Can't get node coords" );
 
    moab::Tag gidTag = mb->globalId_tag();
 
    std::vector< int > gids( cells.size() );
    MB_CHK_SET_ERR( mb->get_coords( cells, ccenters[0], ccenters[1], ccenters[2] ), "Failed to get coordinates" );
    MB_CHK_SET_ERR( mb->tag_get_data( gidTag, cells, &gids[0] ), "Can't set global_id tag" );
 
    moab::Range dualverts( startv, startv + cells.size() - 1 );
    MB_CHK_SET_ERR( mb->add_entities( dual_set, dualverts ), "Can't add entities" );
    MB_CHK_SET_ERR( mb->tag_set_data( gidTag, dualverts, &gids[0] ), "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
        MB_CHK_SET_ERR( mb->get_adjacencies( &vtx, 1, 2, true, adjs ), "Failed to get adjacencies" );
        const size_t nEdges = adjs.size();
 
        double vxyz[3];
        MB_CHK_SET_ERR( mb->get_coords( &vtx, 1, vxyz ), "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;
 
        MB_CHK_SET_ERR( iface->get_element_connect( nElePerType, eit->first, moab::MBPOLYGON, 0, starte, conn ),
                        "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
    MB_CHK_SET_ERR( mb->add_entities( dual_set, dualcells ), "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() );
    MB_CHK_SET_ERR( mb->tag_get_data( gidTag, verts, &gids[0] ), "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 );
 
        MB_CHK_SET_ERR( pcomm.assign_global_ids( dual_set, 2, 1, false, true, true ), "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
        MB_CHK_SET_ERR( mb->tag_set_data( gidTag, dualcells, &gids[0] ), "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
    MB_CHK_SET_ERR( mb->delete_entities( cells ), "Can't remove entities" );
    MB_CHK_SET_ERR( mb->delete_entities( verts ), "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;
    MB_CHK_SET_ERR( mb->create_meshset( moab::MESHSET_SET, triangle_set ), "Can't create new set" );
    MB_CHK_SET_ERR( mb->create_meshset( moab::MESHSET_SET, dual_set ), "Can't create new set" );
 
    // This file is in the mesh files directory
    const char* readopts = "";
    MB_CHK_SET_ERR( mb->load_file( inputFile.c_str(), &triangle_set, readopts ), "Failed to read" );
 
    // get all cells of dimension 2;
    moab::Range cells;
    MB_CHK_SET_ERR( mb->get_entities_by_dimension( triangle_set, 2, cells ), "Failed to get cells" );
 
    std::cout << "Original number of triangular cells : " << cells.size() << "\n";
 
    // call the routine to compute the dual grid
    MB_CHK_SET_ERR( compute_dual_mesh( mb, dual_set, cells ), "Failed to compute dual mesh" );
 
    // write the mesh to disk
    MB_CHK_SET_ERR( mb->write_file( outputFile.c_str() ), "Failed to write new file" );
 
    cells.clear();
    MB_CHK_SET_ERR( mb->get_entities_by_dimension( dual_set, 2, cells ), "Failed to get cells" );
    std::cout << "Wrote the dual mesh: " << outputFile << " containing " << cells.size() << " polygonal cells\n";
 
    delete mb;
 
    return 0;
}