Remapper.hpp

Go to the documentation of this file.

/*
 * =====================================================================================
 *
 *       Filename:  Remapper.hpp
 *
 *    Description:  Interface to the a general remapping capability on arbitrary topology
 *                  that performs both mesh intersection between a source and target grid,
 *                  with arbitrary decompositions. The intersections can then be used to
 *                  either evaluate interpolation weights or to perform high-order
 *                  conservative remapping of solutions defined on the source grid.
 *
 *         Author:  Vijay S. Mahadevan (vijaysm), [email protected]
 *
 * =====================================================================================
 */
 
#ifndef MB_REMAPPER_HPP
#define MB_REMAPPER_HPP
 
#include <string>
 
#include "moab/Interface.hpp"
#ifdef MOAB_HAVE_MPI
#include "moab/ParallelComm.hpp"
#endif
 
// Tempest includes
#ifdef MOAB_HAVE_TEMPESTREMAP
#include "netcdfcpp.h"
#include "TempestRemapAPI.h"
#else
#error "This tool depends on TempestRemap library. Reconfigure using --with-tempestremap"
#endif
 
namespace moab
{
 
class Remapper
{
  public:
#ifdef MOAB_HAVE_MPI
    Remapper( moab::Interface* mbInt, moab::ParallelComm* pcomm = NULL ) : m_interface( mbInt ), m_pcomm( pcomm )
#else
    Remapper( moab::Interface* mbInt ) : m_interface( mbInt )
#endif
    {
    }
 
    virtual ~Remapper()
    {
#ifdef MOAB_HAVE_MPI
        m_pcomm = NULL;
#endif
        m_interface = NULL;
    }
 
    enum IntersectionContext
    {
        DEFAULT              = -1,  // default context
        SourceMesh           = 0,   // source mesh
        TargetMesh           = 1,   // target mesh
        OverlapMesh          = 2,   // overlap/intersection mesh
        CoveringMesh         = 3,   // source mesh covering target mesh
        SourceMeshWithGhosts = 4,   // mesh with extra ghost layers to compute coverage in high order case or bilin
        TargetMeshWithGhosts = 5    // mesh with extra ghost layers to impose target data limiting
    };
 
    moab::Interface* get_interface()
    {
        return m_interface;
    }
 
#ifdef MOAB_HAVE_MPI
    moab::ParallelComm* get_parallel_communicator()
    {
        return m_pcomm;
    }
 
    /// <summary>
    ///     ghost layers
    /// </summary>
    moab::ErrorCode GhostLayers( moab::EntityHandle& meshset, const int ngh_layers, moab::EntityHandle& set_with_ghosts )
    {
        // meshset contains the mesh set distributed already
        //
        moab::ErrorCode rval = m_interface->create_meshset( MESHSET_SET, set_with_ghosts );MB_CHK_ERR( rval );
        // copy original content of mesh set here; we will use it later for local area, for example
        // it will not have any ghosts in it
        moab::Range orgEnts;
        rval = m_interface->get_entities_by_handle( meshset, orgEnts );MB_CHK_ERR( rval );
        rval = m_interface->add_entities( set_with_ghosts, orgEnts );MB_CHK_ERR( rval );
        rval = m_pcomm->exchange_ghost_cells( 2, 0, 1, 0, true, true, &set_with_ghosts );MB_CHK_ERR( rval );
        for( int i = 2; i <= ngh_layers; i++ )
        {
            rval = m_pcomm->correct_thin_ghost_layers();MB_CHK_ERR( rval );
            rval = m_pcomm->exchange_ghost_cells( 2, 0, i, 0, true, true, &set_with_ghosts );MB_CHK_ERR( rval );
        }
 
        // need to set global id tags
        // need also to propagate global id to ghost cells; it is not done by default :(
        moab::Tag gtag = m_interface->globalId_tag();
        moab::Range entities;
        rval = m_interface->get_entities_by_dimension( set_with_ghosts, 2, entities );MB_CHK_ERR( rval );
 
        moab::Tag doftag;
        rval = m_interface->tag_get_handle( "GLOBAL_DOFS", doftag );
        if ( rval == MB_SUCCESS )
        {
            moab::Range quads = entities.subset_by_type(moab::MBQUAD);
            rval = m_pcomm->exchange_tags( doftag, quads );MB_CHK_ERR( rval );
        }
 
        // get all vertices too, need to exchange global ids for vertices too
        moab::Range vertices;
        rval = m_interface->get_connectivity( entities, vertices );MB_CHK_ERR( rval );
        entities.merge( vertices );
        rval = m_pcomm->exchange_tags( gtag, entities );MB_CHK_ERR( rval );
        return rval;
    }
 
#endif
 
    ErrorCode LoadNativeMesh( std::string filename,
                              moab::EntityHandle& meshset,
                              std::vector< int >& metadata,
                              const char* readopts = 0 )
    {
#ifdef MOAB_HAVE_MPI
        std::string opts = "";
        if( readopts )
        {
            if( opts.size() )
                opts = opts + ";" + std::string( readopts );
            else
                opts = std::string( readopts );
        }
 
        if( !m_pcomm->rank() ) std::cout << "Reading file (" << filename << ") with options = [" << opts << "]\n";
#else
        const std::string opts = std::string( ( readopts ? readopts : "" ) );
        std::cout << "Reading file (" << filename << ") with options = [" << opts << "]\n";
#endif
        moab::ErrorCode rval = m_interface->load_file( filename.c_str(), &meshset, opts.c_str() );MB_CHK_ERR( rval );
 
        Tag rectilinearTag;
        rval = m_interface->tag_get_handle( "ClimateMetadata", rectilinearTag );
 
        if( rval != MB_FAILURE && rval != MB_TAG_NOT_FOUND && rval != MB_ALREADY_ALLOCATED &&
            rectilinearTag != nullptr )
        {
            int dimSizes[3];
            rval = m_interface->tag_get_data( rectilinearTag, &meshset, 1,
                                              dimSizes );  // MB_CHK_SET_ERR( rval, "Error geting tag data" );
            metadata.clear();
            metadata.push_back( dimSizes[0] );
            metadata.push_back( dimSizes[1] );
            metadata.push_back( dimSizes[2] );
        }
 
        return MB_SUCCESS;
    }
 
  protected:
    // member data
    Interface* m_interface;
 
#ifdef MOAB_HAVE_MPI
    ParallelComm* m_pcomm;
#endif
};
 
}  // namespace moab
 
#endif /* MB_REMAPPER_HPP */