Cppcheck report - MOAB: tools/mbtempest.cpp

/*
 * Usage: MOAB-Tempest tool
 *
 * Generate a Cubed-Sphere mesh: ./mbtempest -t 0 -res 25 -f cubed_sphere_mesh.exo
 * Generate a RLL mesh: ./mbtempest -t 1 -res 25 -f rll_mesh.exo
 * Generate a Icosahedral-Sphere mesh: ./mbtempest -t 2 -res 25 <-dual> -f icosa_mesh.exo
 *
 * Now you can compute the intersections between the meshes too!
 *
 * Generate the overlap mesh: ./mbtempest -t 3 -l cubed_sphere_mesh.exo -l rll_mesh.exo -f
 * overlap_mesh.exo
 *
 */

#include <iostream>
#include <iomanip>
#include <cstdlib>
#include <vector>
#include <string>
#include <sstream>
#include <cassert>

#include "moab/Core.hpp"
#include "moab/IntxMesh/IntxUtils.hpp"
#include "moab/Remapping/TempestRemapper.hpp"
#include "moab/Remapping/TempestOnlineMap.hpp"
#include "moab/ProgOptions.hpp"
#include "moab/CpuTimer.hpp"
#include "DebugOutput.hpp"

//#ifndef MOAB_HAVE_MPI
//    #error mbtempest tool requires MPI configuration
//#endif

#ifdef MOAB_HAVE_MPI
// MPI includes
#include "moab_mpi.h"
#include "moab/ParallelComm.hpp"
#include "MBParallelConventions.h"
#endif

struct ToolContext
{
    moab::Interface* mbcore;
#ifdef MOAB_HAVE_MPI
    moab::ParallelComm* pcomm;
#endif
    const int proc_id, n_procs;
    moab::DebugOutput outputFormatter;
    int blockSize;
    std::vector< std::string > inFilenames;
    std::vector< Mesh* > meshes;
    std::vector< moab::EntityHandle > meshsets;
    std::vector< int > disc_orders;
    std::vector< std::string > disc_methods;
    std::vector< std::string > doftag_names;
    std::string outFilename;
    std::string intxFilename;
    std::string baselineFile;
    moab::TempestRemapper::TempestMeshType meshType;
    bool computeDual;
    bool computeWeights;
    bool verifyWeights;
    bool enforceConvexity;
    int ensureMonotonicity;
    bool rrmGrids;
    bool kdtreeSearch;
    bool fCheck;
    bool fVolumetric;
    bool fInverseDistanceMap;
    GenerateOfflineMapAlgorithmOptions mapOptions;

#ifdef MOAB_HAVE_MPI
    ToolContext( moab::Interface* icore, moab::ParallelComm* p_pcomm )<--- Member variable 'ToolContext::timer_ops' is not initialized in the constructor.
        : mbcore( icore ), pcomm( p_pcomm ), proc_id( pcomm->rank() ), n_procs( pcomm->size() ),
          outputFormatter( std::cout, pcomm->rank(), 0 ),
#else
    ToolContext( moab::Interface* icore )<--- Member variable 'ToolContext::timer_ops' is not initialized in the constructor.
<--- Struct 'ToolContext' has a constructor with 1 argument that is not explicit. [+]
Struct 'ToolContext' has a constructor with 1 argument that is not explicit. Such constructors should in general be explicit for type safety reasons. Using the explicit keyword in the constructor means some mistakes when using the class can be avoided.
        : mbcore( icore ), proc_id( 0 ), n_procs( 1 ), outputFormatter( std::cout, 0, 0 ),
#endif
          blockSize( 5 ), outFilename( "outputFile.nc" ), intxFilename( "intxFile.h5m" ), baselineFile( "" ),
          meshType( moab::TempestRemapper::DEFAULT ), computeDual( false ), computeWeights( false ),
          verifyWeights( false ), enforceConvexity( false ), ensureMonotonicity( 0 ), rrmGrids( false ),
          kdtreeSearch( true ), fCheck( n_procs > 1 ? false : true ), fVolumetric( false ), fInverseDistanceMap( false )
    {
        inFilenames.resize( 2 );
        doftag_names.resize( 2 );
        timer = new moab::CpuTimer();<--- Struct 'ToolContext' does not have a copy constructor which is recommended since it has dynamic memory/resource allocation(s).<--- Struct 'ToolContext' does not have a operator= which is recommended since it has dynamic memory/resource allocation(s).

        outputFormatter.set_prefix( "[MBTempest]: " );
    }

    ~ToolContext()
    {
        // for (unsigned i=0; i < meshes.size(); ++i) delete meshes[i];
        meshes.clear();
        inFilenames.clear();
        disc_orders.clear();
        disc_methods.clear();
        doftag_names.clear();
        outFilename.clear();
        intxFilename.clear();
        baselineFile.clear();
        meshsets.clear();
        delete timer;
    }

    void timer_push( std::string operation )
<--- Function parameter 'operation' should be passed by const reference. [+]
Parameter 'operation' is passed by value. It could be passed as a const reference which is usually faster and recommended in C++.
    {
        timer_ops = timer->time_since_birth();
        opName    = operation;
    }

    void timer_pop()
    {
        double locElapsed = timer->time_since_birth() - timer_ops, avgElapsed = 0, maxElapsed = 0;
#ifdef MOAB_HAVE_MPI
        MPI_Reduce( &locElapsed, &maxElapsed, 1, MPI_DOUBLE, MPI_MAX, 0, pcomm->comm() );
        MPI_Reduce( &locElapsed, &avgElapsed, 1, MPI_DOUBLE, MPI_SUM, 0, pcomm->comm() );
#else
        maxElapsed = locElapsed;
        avgElapsed = locElapsed;
#endif
        if( !proc_id )
        {
            avgElapsed /= n_procs;
            std::cout << "[LOG] Time taken to " << opName.c_str() << ": max = " << maxElapsed
                      << ", avg = " << avgElapsed << "\n";
        }
        // std::cout << "\n[LOG" << proc_id << "] Time taken to " << opName << " = " <<
        // timer->time_since_birth() - timer_ops << std::endl;
        opName.clear();
    }

    void ParseCLOptions( int argc, char* argv[] )
    {
        ProgOptions opts;
        int imeshType                  = 0;
        std::string expectedFName      = "output.exo";
        std::string expectedMethod     = "fv";
        std::string expectedDofTagName = "GLOBAL_ID";
        int expectedOrder              = 1;

        if( !proc_id )
        {
            std::cout << "Command line options provided to mbtempest:\n  ";
            for( int iarg = 0; iarg < argc; ++iarg )
                std::cout << argv[iarg] << " ";
            std::cout << std::endl << std::endl;
        }

        opts.addOpt< int >( "type,t",
                            "Type of mesh (default=CS; Choose from [CS=0, RLL=1, ICO=2, OVERLAP_FILES=3, "
                            "OVERLAP_MEMORY=4, OVERLAP_MOAB=5])",
                            &imeshType );

        opts.addOpt< int >( "res,r", "Resolution of the mesh (default=5)", &blockSize );

        opts.addOpt< void >( "dual,d", "Output the dual of the mesh (relevant only for ICO mesh type)", &computeDual );

        opts.addOpt< std::string >( "file,f", "Output computed mesh or remapping weights to specified filename",
                                    &outFilename );
        opts.addOpt< std::string >(
            "load,l", "Input mesh filenames for source and target meshes. (relevant only when computing weights)",
            &expectedFName );

        opts.addOpt< void >( "advfront,a",
                             "Use the advancing front intersection instead of the Kd-tree based algorithm "
                             "to compute mesh intersections. (relevant only when computing weights)" );

        opts.addOpt< std::string >( "intx,i", "Output TempestRemap intersection mesh filename", &intxFilename );

        opts.addOpt< void >( "weights,w",
                             "Compute and output the weights using the overlap mesh (generally "
                             "relevant only for OVERLAP mesh)",
                             &computeWeights );

        opts.addOpt< std::string >( "method,m", "Discretization method for the source and target solution fields",
                                    &expectedMethod );

        opts.addOpt< int >( "order,o", "Discretization orders for the source and target solution fields",
                            &expectedOrder );

        opts.addOpt< std::string >( "global_id,g",
                                    "Tag name that contains the global DoF IDs for source and target solution fields",
                                    &expectedDofTagName );

        opts.addOpt< void >( "noconserve",
                             "Do not apply conservation to the resultant weights (relevant only "
                             "when computing weights)",
                             &mapOptions.fNoConservation );

        opts.addOpt< void >( "volumetric",
                             "Apply a volumetric projection to compute the weights (relevant only "
                             "when computing weights)",
                             &fVolumetric );

        opts.addOpt< void >( "invdist",
                             "Apply a inverse distance weight projection to compute the weights (relevant only "
                             "when computing FV-FV weights)",
                             &fInverseDistanceMap );

        opts.addOpt< int >( "monotonicity", "Ensure monotonicity in the weight generation. Options=[0,1,2,3]",
                            &ensureMonotonicity );

        opts.addOpt< void >( "enforce_convexity", "check convexity of input meshes to compute mesh intersections",
                             &enforceConvexity );

        opts.addOpt< void >( "nobubble", "do not use bubble on interior of spectral element nodes",
                             &mapOptions.fNoBubble );

        opts.addOpt< void >( "sparseconstraints", "do not use bubble on interior of spectral element nodes",
                             &mapOptions.fSparseConstraints );

        opts.addOpt< void >( "rrmgrids",
                             "At least one of the meshes is a regionally refined grid (relevant to "
                             "accelerate intersection computation)",
                             &rrmGrids );

        opts.addOpt< void >( "checkmap", "Check the generated map for conservation and consistency", &fCheck );

        opts.addOpt< void >( "verify",
                             "Verify the accuracy of the maps by projecting analytical functions "
                             "from source to target "
                             "grid by applying the maps",
                             &verifyWeights );

        opts.addOpt< std::string >( "baseline", "Output baseline file", &baselineFile );

        opts.parseCommandLine( argc, argv );

        // By default - use Kd-tree based search; if user asks for advancing front, disable Kd-tree
        // algorithm
        kdtreeSearch = opts.numOptSet( "advfront,a" ) == 0;

        switch( imeshType )
        {
            case 0:
                meshType = moab::TempestRemapper::CS;
                break;

            case 1:
                meshType = moab::TempestRemapper::RLL;
                break;

            case 2:
                meshType = moab::TempestRemapper::ICO;
                break;

            case 3:
                meshType = moab::TempestRemapper::OVERLAP_FILES;
                break;

            case 4:
                meshType = moab::TempestRemapper::OVERLAP_MEMORY;
                break;

            case 5:
                meshType = moab::TempestRemapper::OVERLAP_MOAB;
                break;

            default:
                meshType = moab::TempestRemapper::DEFAULT;
                break;
        }

        if( meshType > moab::TempestRemapper::ICO ) // compute overlap mesh and maps possibly
        {
            opts.getOptAllArgs( "load,l", inFilenames );
            opts.getOptAllArgs( "order,o", disc_orders );
            opts.getOptAllArgs( "method,m", disc_methods );
            opts.getOptAllArgs( "global_id,i", doftag_names );

            assert( inFilenames.size() == 2 );
            assert( disc_orders.size() == 2 );
            assert( disc_methods.size() == 2 );
            assert( ensureMonotonicity >= 0 && ensureMonotonicity <= 3 );

            // get discretization order parameters
            if( disc_orders.size() == 0 ) disc_orders.resize( 2, 1 );
            if( disc_orders.size() == 1 ) disc_orders.push_back( 1 );

            // get discretization method parameters
            if( disc_methods.size() == 0 ) disc_methods.resize( 2, "fv" );
            if( disc_methods.size() == 1 ) disc_methods.push_back( "fv" );

            // get DoF tagname parameters
            if( doftag_names.size() == 0 ) doftag_names.resize( 2, "GLOBAL_ID" );
            if( doftag_names.size() == 1 ) doftag_names.push_back( "GLOBAL_ID" );

            // for computing maps and overlaps, set discretization orders
            mapOptions.nPin           = disc_orders[0];
            mapOptions.nPout          = disc_orders[1];
            mapOptions.fSourceConcave = false;
            mapOptions.fTargetConcave = false;

            mapOptions.strMethod = "";
            switch( ensureMonotonicity )
            {
                case 0:
                    mapOptions.fMonotone = false;
                    break;
                case 3:
                    mapOptions.strMethod += "mono3;";
                    break;
                case 2:
                    mapOptions.strMethod += "mono2;";
                    break;
                default:
                    mapOptions.fMonotone = true;
            }
            mapOptions.fNoCorrectAreas = false;
            mapOptions.fNoCheck        = !fCheck;

            //assert( fVolumetric && fInverseDistanceMap == false );  // both options cannot be active
            if( fVolumetric ) mapOptions.strMethod += "volumetric;";
            if( fInverseDistanceMap ) mapOptions.strMethod += "invdist;";
        }

        // clear temporary string name
        expectedFName.clear();

        mapOptions.strOutputMapFile = outFilename;
        mapOptions.strOutputFormat  = "Netcdf4";
    }

  private:
    moab::CpuTimer* timer;
    double timer_ops;
    std::string opName;
};

// Forward declare some methods
static moab::ErrorCode CreateTempestMesh( ToolContext&, moab::TempestRemapper& remapper, Mesh* );
inline double sample_slow_harmonic( double dLon, double dLat );
inline double sample_fast_harmonic( double dLon, double dLat );
inline double sample_constant( double dLon, double dLat );
inline double sample_stationary_vortex( double dLon, double dLat );

int main( int argc, char* argv[] )
{
    moab::ErrorCode rval;
    NcError error( NcError::verbose_nonfatal );
    std::stringstream sstr;

    int proc_id = 0, nprocs = 1;<--- 'nprocs' is assigned value '1' here.<--- 'nprocs' is assigned value '1' here.
#ifdef MOAB_HAVE_MPI
    MPI_Init( &argc, &argv );
    MPI_Comm_rank( MPI_COMM_WORLD, &proc_id );
    MPI_Comm_size( MPI_COMM_WORLD, &nprocs );
#endif

    moab::Interface* mbCore = new( std::nothrow ) moab::Core;

    if( NULL == mbCore )
    {
        return 1;
    }

    ToolContext* runCtx;
#ifdef MOAB_HAVE_MPI
    moab::ParallelComm* pcomm = new moab::ParallelComm( mbCore, MPI_COMM_WORLD, 0 );

    runCtx                   = new ToolContext( mbCore, pcomm );
    const char* writeOptions = ( nprocs > 1 ? "PARALLEL=WRITE_PART" : "" );
#else
    runCtx                   = new ToolContext( mbCore );
    const char* writeOptions = "";
#endif
    runCtx->ParseCLOptions( argc, argv );

    const double radius_src  = 1.0 /*2.0*acos(-1.0)*/;
    const double radius_dest = 1.0 /*2.0*acos(-1.0)*/;

    moab::DebugOutput& outputFormatter = runCtx->outputFormatter;

#ifdef MOAB_HAVE_MPI
    moab::TempestRemapper remapper( mbCore, pcomm );
#else
    moab::TempestRemapper remapper( mbCore );
#endif
    remapper.meshValidate     = true;
    remapper.constructEdgeMap = false;
    remapper.initialize();

    // Default area_method = lHuiller; Options: Girard, GaussQuadrature (if TR is available)
#ifdef MOAB_HAVE_TEMPESTREMAP
    moab::IntxAreaUtils areaAdaptor( moab::IntxAreaUtils::GaussQuadrature );
#else
    moab::IntxAreaUtils areaAdaptor( moab::IntxAreaUtils::lHuiller );
#endif

    Mesh* tempest_mesh = new Mesh();
    runCtx->timer_push( "create Tempest mesh" );
    rval = CreateTempestMesh( *runCtx, remapper, tempest_mesh );MB_CHK_ERR( rval );
    runCtx->timer_pop();

    const double epsrel = ReferenceTolerance;  // ReferenceTolerance is defined in Defines.h in tempestremap;
                                               // Defines.h is included in SparseMatrix.h
                                               // SparseMatrix.h is included in OfflineMap.h
                                               // OfflineMap.h is included in TempestOnlineMap.hpp
                                               // TempestOnlineMap.hpp is included in this file, and is part of MOAB
    // Some constant parameters

    const double boxeps = 1e-6;

    if( runCtx->meshType == moab::TempestRemapper::OVERLAP_MEMORY )
    {
        // Compute intersections with MOAB
        // For the overlap method, choose between: "fuzzy", "exact" or "mixed"
        assert( runCtx->meshes.size() == 3 );

#ifdef MOAB_HAVE_MPI
        rval = pcomm->check_all_shared_handles();MB_CHK_ERR( rval );
#endif

        // Load the meshes and validate
        rval = remapper.ConvertTempestMesh( moab::Remapper::SourceMesh );MB_CHK_ERR( rval );
        rval = remapper.ConvertTempestMesh( moab::Remapper::TargetMesh );MB_CHK_ERR( rval );
        remapper.SetMeshType( moab::Remapper::OverlapMesh, moab::TempestRemapper::OVERLAP_FILES );
        rval = remapper.ConvertTempestMesh( moab::Remapper::OverlapMesh );MB_CHK_ERR( rval );
        rval = mbCore->write_mesh( "tempest_intersection.h5m", &runCtx->meshsets[2], 1 );MB_CHK_ERR( rval );

        // print verbosely about the problem setting
        {
            moab::Range rintxverts, rintxelems;
            rval = mbCore->get_entities_by_dimension( runCtx->meshsets[0], 0, rintxverts );MB_CHK_ERR( rval );
            rval = mbCore->get_entities_by_dimension( runCtx->meshsets[0], 2, rintxelems );MB_CHK_ERR( rval );
            outputFormatter.printf( 0, "The source set contains %lu vertices and %lu elements \n", rintxverts.size(),
                                    rintxelems.size() );

            moab::Range bintxverts, bintxelems;
            rval = mbCore->get_entities_by_dimension( runCtx->meshsets[1], 0, bintxverts );MB_CHK_ERR( rval );
            rval = mbCore->get_entities_by_dimension( runCtx->meshsets[1], 2, bintxelems );MB_CHK_ERR( rval );
            outputFormatter.printf( 0, "The target set contains %lu vertices and %lu elements \n", bintxverts.size(),
                                    bintxelems.size() );
        }

        moab::EntityHandle intxset;  // == remapper.GetMeshSet(moab::Remapper::OverlapMesh);

        // Compute intersections with MOAB
        {
            // Create the intersection on the sphere object
            runCtx->timer_push( "setup the intersector" );

            moab::Intx2MeshOnSphere* mbintx = new moab::Intx2MeshOnSphere( mbCore );
            mbintx->set_error_tolerance( epsrel );
            mbintx->set_box_error( boxeps );
            mbintx->set_radius_source_mesh( radius_src );
            mbintx->set_radius_destination_mesh( radius_dest );
#ifdef MOAB_HAVE_MPI
            mbintx->set_parallel_comm( pcomm );
#endif
            rval = mbintx->FindMaxEdges( runCtx->meshsets[0], runCtx->meshsets[1] );MB_CHK_ERR( rval );

#ifdef MOAB_HAVE_MPI
            moab::Range local_verts;
            rval = mbintx->build_processor_euler_boxes( runCtx->meshsets[1], local_verts );MB_CHK_ERR( rval );

            runCtx->timer_pop();

            moab::EntityHandle covering_set;
            runCtx->timer_push( "communicate the mesh" );
            rval = mbintx->construct_covering_set( runCtx->meshsets[0], covering_set );MB_CHK_ERR( rval );  // lots of communication if mesh is distributed very differently
            runCtx->timer_pop();
#else
            moab::EntityHandle covering_set = runCtx->meshsets[0];
#endif
            // Now let's invoke the MOAB intersection algorithm in parallel with a
            // source and target mesh set representing two different decompositions
            runCtx->timer_push( "compute intersections with MOAB" );
            rval = mbCore->create_meshset( moab::MESHSET_SET, intxset );MB_CHK_SET_ERR( rval, "Can't create new set" );
            rval = mbintx->intersect_meshes( covering_set, runCtx->meshsets[1], intxset );MB_CHK_SET_ERR( rval, "Can't compute the intersection of meshes on the sphere" );
            runCtx->timer_pop();

            // free the memory
            delete mbintx;
        }

        {
            moab::Range intxelems, intxverts;
            rval = mbCore->get_entities_by_dimension( intxset, 2, intxelems );MB_CHK_ERR( rval );
            rval = mbCore->get_entities_by_dimension( intxset, 0, intxverts, true );MB_CHK_ERR( rval );
            outputFormatter.printf( 0, "The intersection set contains %lu elements and %lu vertices \n",
                                    intxelems.size(), intxverts.size() );

            double initial_sarea =
                areaAdaptor.area_on_sphere( mbCore, runCtx->meshsets[0],
                                            radius_src );  // use the target to compute the initial area
            double initial_tarea =
                areaAdaptor.area_on_sphere( mbCore, runCtx->meshsets[1],
                                            radius_dest );  // use the target to compute the initial area
            double intx_area = areaAdaptor.area_on_sphere( mbCore, intxset, radius_src );

            outputFormatter.printf( 0, "mesh areas: source = %12.10f, target = %12.10f, intersection = %12.10f \n",
                                    initial_sarea, initial_tarea, intx_area );
            outputFormatter.printf( 0, "relative error w.r.t source = %12.10e, target = %12.10e \n",
                                    fabs( intx_area - initial_sarea ) / initial_sarea,
                                    fabs( intx_area - initial_tarea ) / initial_tarea );
        }

        // Write out our computed intersection file
        rval = mbCore->write_mesh( "moab_intersection.h5m", &intxset, 1 );MB_CHK_ERR( rval );

        if( runCtx->computeWeights )
        {
            runCtx->timer_push( "compute weights with the Tempest meshes" );
            // Call to generate an offline map with the tempest meshes
            OfflineMap weightMap;
            int err = GenerateOfflineMapWithMeshes( *runCtx->meshes[0], *runCtx->meshes[1], *runCtx->meshes[2],
                                                    runCtx->disc_methods[0],  // std::string strInputType
                                                    runCtx->disc_methods[1],  // std::string strOutputType,
                                                    runCtx->mapOptions, weightMap );
            runCtx->timer_pop();

            std::map< std::string, std::string > mapAttributes;
            if( err )
            {
                rval = moab::MB_FAILURE;
            }
            else
            {
                weightMap.Write( "outWeights.nc", mapAttributes );
            }
        }
    }
    else if( runCtx->meshType == moab::TempestRemapper::OVERLAP_MOAB )
    {
        // Usage: mpiexec -n 2 tools/mbtempest -t 5 -l mycs_2.h5m -l myico_2.h5m -f myoverlap_2.h5m
#ifdef MOAB_HAVE_MPI
        rval = pcomm->check_all_shared_handles();MB_CHK_ERR( rval );
#endif
        // print verbosely about the problem setting
        {
            moab::Range rintxverts, rintxelems;
            rval = mbCore->get_entities_by_dimension( runCtx->meshsets[0], 0, rintxverts );MB_CHK_ERR( rval );
            rval = mbCore->get_entities_by_dimension( runCtx->meshsets[0], 2, rintxelems );MB_CHK_ERR( rval );
            rval = moab::IntxUtils::fix_degenerate_quads( mbCore, runCtx->meshsets[0] );MB_CHK_ERR( rval );
            if( runCtx->enforceConvexity )
            {
                rval = moab::IntxUtils::enforce_convexity( mbCore, runCtx->meshsets[0], proc_id );MB_CHK_ERR( rval );
            }
            rval = areaAdaptor.positive_orientation( mbCore, runCtx->meshsets[0], radius_src );MB_CHK_ERR( rval );
            if( !proc_id )
                outputFormatter.printf( 0, "The source set contains %lu vertices and %lu elements \n",
                                        rintxverts.size(), rintxelems.size() );

            moab::Range bintxverts, bintxelems;
            rval = mbCore->get_entities_by_dimension( runCtx->meshsets[1], 0, bintxverts );MB_CHK_ERR( rval );
            rval = mbCore->get_entities_by_dimension( runCtx->meshsets[1], 2, bintxelems );MB_CHK_ERR( rval );
            rval = moab::IntxUtils::fix_degenerate_quads( mbCore, runCtx->meshsets[1] );MB_CHK_ERR( rval );
            if( runCtx->enforceConvexity )
            {
                rval = moab::IntxUtils::enforce_convexity( mbCore, runCtx->meshsets[1], proc_id );MB_CHK_ERR( rval );
            }
            rval = areaAdaptor.positive_orientation( mbCore, runCtx->meshsets[1], radius_dest );MB_CHK_ERR( rval );
            if( !proc_id )
                outputFormatter.printf( 0, "The target set contains %lu vertices and %lu elements \n",
                                        bintxverts.size(), bintxelems.size() );
        }

        // First compute the covering set such that the target elements are fully covered by the
        // lcoal source grid
        runCtx->timer_push( "construct covering set for intersection" );
        rval = remapper.ConstructCoveringSet( epsrel, 1.0, 1.0, boxeps, runCtx->rrmGrids );MB_CHK_ERR( rval );
        runCtx->timer_pop();

        // Compute intersections with MOAB with either the Kd-tree or the advancing front algorithm
        runCtx->timer_push( "setup and compute mesh intersections" );
        rval = remapper.ComputeOverlapMesh( runCtx->kdtreeSearch, false );MB_CHK_ERR( rval );
        runCtx->timer_pop();

        // print some diagnostic checks to see if the overlap grid resolved the input meshes
        // correctly
        {
            double local_areas[3],
                global_areas[3];  // Array for Initial area, and through Method 1 and Method 2
            // local_areas[0] = area_on_sphere_lHuiller ( mbCore, runCtx->meshsets[1], radius_src );
            local_areas[0] = areaAdaptor.area_on_sphere( mbCore, runCtx->meshsets[0], radius_src );
            local_areas[1] = areaAdaptor.area_on_sphere( mbCore, runCtx->meshsets[1], radius_dest );
            local_areas[2] = areaAdaptor.area_on_sphere( mbCore, runCtx->meshsets[2], radius_src );

#ifdef MOAB_HAVE_MPI
            MPI_Allreduce( &local_areas[0], &global_areas[0], 3, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD );
#else
            global_areas[0]                 = local_areas[0];
            global_areas[1]                 = local_areas[1];
            global_areas[2]                 = local_areas[2];
#endif
            if( !proc_id )
            {
                outputFormatter.printf( 0,
                                        "initial area: source mesh = %12.14f, target mesh = "
                                        "%12.14f, overlap mesh = %12.14f\n",
                                        global_areas[0], global_areas[1], global_areas[2] );
                // outputFormatter.printf ( 0, " area with l'Huiller: %12.14f with Girard:
                // %12.14f\n", global_areas[2], global_areas[3] ); outputFormatter.printf ( 0, "
                // relative difference areas = %12.10e\n", fabs ( global_areas[2] - global_areas[3]
                // ) / global_areas[2] );
                outputFormatter.printf( 0, "relative error w.r.t source = %12.14e, and target = %12.14e\n",
                                        fabs( global_areas[0] - global_areas[2] ) / global_areas[0],
                                        fabs( global_areas[1] - global_areas[2] ) / global_areas[1] );
            }
        }

        if( runCtx->intxFilename.size() )
        {
            moab::EntityHandle writableOverlapSet;
            rval = mbCore->create_meshset( moab::MESHSET_SET, writableOverlapSet );MB_CHK_SET_ERR( rval, "Can't create new set" );
            moab::EntityHandle meshOverlapSet = remapper.GetMeshSet( moab::Remapper::OverlapMesh );
            moab::Range ovEnts;
            rval = mbCore->get_entities_by_dimension( meshOverlapSet, 2, ovEnts );MB_CHK_SET_ERR( rval, "Can't create new set" );
            rval = mbCore->get_entities_by_dimension( meshOverlapSet, 0, ovEnts );MB_CHK_SET_ERR( rval, "Can't create new set" );

#ifdef MOAB_HAVE_MPI
            // Do not remove ghosted entities if we still haven't computed weights
            // Remove ghosted entities from overlap set before writing the new mesh set to file
            if( nprocs > 1 )
            {
                moab::Range ghostedEnts;
                rval = remapper.GetOverlapAugmentedEntities( ghostedEnts );MB_CHK_ERR( rval );
                ovEnts = moab::subtract( ovEnts, ghostedEnts );
            }
#endif
            rval = mbCore->add_entities( writableOverlapSet, ovEnts );MB_CHK_SET_ERR( rval, "Deleting ghosted entities failed" );

            size_t lastindex = runCtx->intxFilename.find_last_of( "." );
            sstr.str( "" );
            sstr << runCtx->intxFilename.substr( 0, lastindex ) << ".h5m";
            if( !runCtx->proc_id )
                std::cout << "Writing out the MOAB intersection mesh file to " << sstr.str() << std::endl;

            // Write out our computed intersection file
            rval = mbCore->write_file( sstr.str().c_str(), NULL, writeOptions, &writableOverlapSet, 1 );MB_CHK_ERR( rval );
        }

        if( runCtx->computeWeights )
        {
            runCtx->meshes[2] = remapper.GetMesh( moab::Remapper::OverlapMesh );
            if( !runCtx->proc_id ) std::cout << std::endl;

            runCtx->timer_push( "setup computation of weights" );
            // Call to generate the remapping weights with the tempest meshes
            moab::TempestOnlineMap* weightMap = new moab::TempestOnlineMap( &remapper );
            runCtx->timer_pop();

            runCtx->timer_push( "compute weights with TempestRemap" );
            rval = weightMap->GenerateRemappingWeights(
                runCtx->disc_methods[0],  // std::string strInputType
                runCtx->disc_methods[1],  // std::string strOutputType,
                runCtx->mapOptions,       // const GenerateOfflineMapAlgorithmOptions& options
                runCtx->doftag_names[0],  // const std::string& source_tag_name
                runCtx->doftag_names[1]   // const std::string& target_tag_name
            );MB_CHK_ERR( rval );
            runCtx->timer_pop();

            // Invoke the CheckMap routine on the TempestRemap serial interface directly, if running
            // on a single process
            if( nprocs == 1 )
<--- The comparison 'nprocs == 1' is always true. [+]
Finding the same expression on both sides of an operator is suspicious and might indicate a cut and paste or logic error. Please examine this code carefully to determine if it is correct.
            {
                const double dNormalTolerance = 1.0E-8;
                const double dStrictTolerance = 1.0E-12;
                weightMap->CheckMap( runCtx->fCheck, runCtx->fCheck, runCtx->fCheck && ( runCtx->ensureMonotonicity ),
                                     dNormalTolerance, dStrictTolerance );
            }

            if( runCtx->outFilename.size() )
            {
                // Write the map file to disk in parallel using either HDF5 or SCRIP interface
                rval = weightMap->WriteParallelMap( runCtx->outFilename.c_str() );MB_CHK_ERR( rval );

                // Write out the metadata information for the map file
                if( proc_id == 0 )
                {
                    size_t lastindex = runCtx->outFilename.find_last_of( "." );
                    sstr.str( "" );
                    sstr << runCtx->outFilename.substr( 0, lastindex ) << ".meta";

                    std::ofstream metafile( sstr.str() );
                    metafile << "Generator = MOAB-TempestRemap (mbtempest) Offline Regridding "
                                "Weight Generator"
                             << std::endl;
                    metafile << "domain_a = " << runCtx->inFilenames[0] << std::endl;
                    metafile << "domain_b = " << runCtx->inFilenames[1] << std::endl;
                    metafile << "domain_aNb = "
                             << ( runCtx->intxFilename.size() ? runCtx->intxFilename : "outOverlap.h5m" ) << std::endl;
                    metafile << "map_aPb = " << runCtx->outFilename << std::endl;
                    metafile << "type_src = " << runCtx->disc_methods[0] << std::endl;
                    metafile << "np_src = " << runCtx->disc_orders[0] << std::endl;
                    metafile << "concave_src = " << ( runCtx->mapOptions.fSourceConcave ? "true" : "false" )
                             << std::endl;
                    metafile << "type_dst = " << runCtx->disc_methods[1] << std::endl;
                    metafile << "np_dst = " << runCtx->disc_orders[1] << std::endl;
                    metafile << "concave_dst = " << ( runCtx->mapOptions.fTargetConcave ? "true" : "false" )
                             << std::endl;
                    metafile << "mono_type = " << runCtx->ensureMonotonicity << std::endl;
                    metafile << "bubble = " << ( runCtx->mapOptions.fNoBubble ? "false" : "true" ) << std::endl;
                    metafile << "version = "
                             << "MOAB v5.1.0+" << std::endl;
                    metafile.close();
                }
            }

            if( runCtx->verifyWeights )
            {
                // Let us pick a sampling test function for solution evaluation
                moab::TempestOnlineMap::sample_function testFunction =
                    &sample_stationary_vortex;  // &sample_slow_harmonic;

                runCtx->timer_push( "describe a solution on source grid" );
                moab::Tag srcAnalyticalFunction;
                rval = weightMap->DefineAnalyticalSolution( srcAnalyticalFunction, "AnalyticalSolnSrcExact",
                                                            moab::Remapper::SourceMesh, testFunction );MB_CHK_ERR( rval );
                runCtx->timer_pop();
                // rval = mbCore->write_file ( "srcWithSolnTag.h5m", NULL, writeOptions,
                // &runCtx->meshsets[0], 1 ); MB_CHK_ERR ( rval );

                runCtx->timer_push( "describe a solution on target grid" );
                moab::Tag tgtAnalyticalFunction;
                moab::Tag tgtProjectedFunction;
                rval = weightMap->DefineAnalyticalSolution( tgtAnalyticalFunction, "AnalyticalSolnTgtExact",
                                                            moab::Remapper::TargetMesh, testFunction,
                                                            &tgtProjectedFunction, "ProjectedSolnTgt" );MB_CHK_ERR( rval );
                // rval = mbCore->write_file ( "tgtWithSolnTag.h5m", NULL, writeOptions,
                // &runCtx->meshsets[1], 1 ); MB_CHK_ERR ( rval );
                runCtx->timer_pop();

                runCtx->timer_push( "compute solution projection on target grid" );
                rval = weightMap->ApplyWeights( srcAnalyticalFunction, tgtProjectedFunction );MB_CHK_ERR( rval );
                runCtx->timer_pop();
                rval = mbCore->write_file( "tgtWithSolnTag2.h5m", NULL, writeOptions, &runCtx->meshsets[1], 1 );MB_CHK_ERR( rval );

                if( nprocs == 1 && runCtx->baselineFile.size() )
<--- The comparison 'nprocs == 1' is always true. [+]
Finding the same expression on both sides of an operator is suspicious and might indicate a cut and paste or logic error. Please examine this code carefully to determine if it is correct.
                {
                    // save the field from tgtWithSolnTag2 in a text file, and global ids for cells
                    moab::Range tgtCells;
                    rval = mbCore->get_entities_by_dimension( runCtx->meshsets[1], 2, tgtCells );MB_CHK_ERR( rval );
                    std::vector< int > globIds;
                    globIds.resize( tgtCells.size() );
                    std::vector< double > vals;
                    vals.resize( tgtCells.size() );
                    moab::Tag projTag;
                    rval = mbCore->tag_get_handle( "ProjectedSolnTgt", projTag );MB_CHK_ERR( rval );
                    moab::Tag gid = mbCore->globalId_tag();
                    rval          = mbCore->tag_get_data( gid, tgtCells, &globIds[0] );MB_CHK_ERR( rval );
                    rval = mbCore->tag_get_data( projTag, tgtCells, &vals[0] );MB_CHK_ERR( rval );
                    std::fstream fs;
                    fs.open( runCtx->baselineFile.c_str(), std::fstream::out );
                    fs << std::setprecision( 15 );  // maximum precision for doubles
                    for( size_t i = 0; i < tgtCells.size(); i++ )
                        fs << globIds[i] << " " << vals[i] << "\n";
                    fs.close();
                    // for good measure, save the source file too, with the tag AnalyticalSolnSrcExact
                    // it will be used later to test, along with a target file
                    rval = mbCore->write_file( "srcWithSolnTag.h5m", NULL, writeOptions, &runCtx->meshsets[0], 1 );MB_CHK_ERR( rval );
                }
                runCtx->timer_push( "compute error metrics against analytical solution on target grid" );
                std::map< std::string, double > errMetrics;
                rval = weightMap->ComputeMetrics( moab::Remapper::TargetMesh, tgtAnalyticalFunction,
                                                  tgtProjectedFunction, errMetrics, true );MB_CHK_ERR( rval );
                runCtx->timer_pop();
            }

            delete weightMap;
        }
    }

    // Clean up
    remapper.clear();
    delete runCtx;
    delete mbCore;

#ifdef MOAB_HAVE_MPI
    MPI_Finalize();
#endif
    exit( 0 );
}

static moab::ErrorCode CreateTempestMesh( ToolContext& ctx, moab::TempestRemapper& remapper, Mesh* tempest_mesh )
{
    moab::ErrorCode rval = moab::MB_SUCCESS;
    int err;
    moab::DebugOutput& outputFormatter = ctx.outputFormatter;
<--- The scope of the variable 'outputFormatter' can be reduced. [+]
The scope of the variable 'outputFormatter' can be reduced. Warning: Be careful when fixing this message, especially when there are inner loops. Here is an example where cppcheck will write that the scope for 'i' can be reduced:
void f(int x)
{
    int i = 0;
    if (x) {
        // it's safe to move 'int i = 0;' here
        for (int n = 0; n < 10; ++n) {
            // it is possible but not safe to move 'int i = 0;' here
            do_something(&i);
        }
    }
}
When you see this message it is always safe to reduce the variable scope 1 level.

    if( !ctx.proc_id )
    {
        outputFormatter.printf( 0, "Creating TempestRemap Mesh object ...\n" );
    }

    if( ctx.meshType == moab::TempestRemapper::OVERLAP_FILES )
    {
        // For the overlap method, choose between: "fuzzy", "exact" or "mixed"
        err = GenerateOverlapMesh( ctx.inFilenames[0], ctx.inFilenames[1], *tempest_mesh, ctx.outFilename, "NetCDF4",
                                   "exact", true );

        if( err )
        {
            rval = moab::MB_FAILURE;
        }
        else
        {
            ctx.meshes.push_back( tempest_mesh );
        }
    }
    else if( ctx.meshType == moab::TempestRemapper::OVERLAP_MEMORY )
    {
        // Load the meshes and validate
        ctx.meshsets.resize( 3 );
        ctx.meshes.resize( 3 );
        ctx.meshsets[0] = remapper.GetMeshSet( moab::Remapper::SourceMesh );
        ctx.meshsets[1] = remapper.GetMeshSet( moab::Remapper::TargetMesh );
        ctx.meshsets[2] = remapper.GetMeshSet( moab::Remapper::OverlapMesh );

        // First the source
        rval = remapper.LoadMesh( moab::Remapper::SourceMesh, ctx.inFilenames[0], moab::TempestRemapper::DEFAULT );MB_CHK_ERR( rval );
        ctx.meshes[0] = remapper.GetMesh( moab::Remapper::SourceMesh );

        // Next the target
        rval = remapper.LoadMesh( moab::Remapper::TargetMesh, ctx.inFilenames[1], moab::TempestRemapper::DEFAULT );MB_CHK_ERR( rval );
        ctx.meshes[1] = remapper.GetMesh( moab::Remapper::TargetMesh );

        // Now let us construct the overlap mesh, by calling TempestRemap interface directly
        // For the overlap method, choose between: "fuzzy", "exact" or "mixed"
        err = GenerateOverlapWithMeshes( *ctx.meshes[0], *ctx.meshes[1], *tempest_mesh, "" /*ctx.outFilename*/,
                                         "NetCDF4", "exact", false );

        if( err )
        {
            rval = moab::MB_FAILURE;
        }
        else
        {
            remapper.SetMesh( moab::Remapper::OverlapMesh, tempest_mesh );
            ctx.meshes[2] = remapper.GetMesh( moab::Remapper::OverlapMesh );
            // ctx.meshes.push_back(*tempest_mesh);
        }
    }
    else if( ctx.meshType == moab::TempestRemapper::OVERLAP_MOAB )
    {
        ctx.meshsets.resize( 3 );
        ctx.meshes.resize( 3 );
        ctx.meshsets[0] = remapper.GetMeshSet( moab::Remapper::SourceMesh );
        ctx.meshsets[1] = remapper.GetMeshSet( moab::Remapper::TargetMesh );
        ctx.meshsets[2] = remapper.GetMeshSet( moab::Remapper::OverlapMesh );

        const double radius_src  = 1.0 /*2.0*acos(-1.0)*/;
        const double radius_dest = 1.0 /*2.0*acos(-1.0)*/;

        const char* additional_read_opts = ( ctx.n_procs > 1 ? "NO_SET_CONTAINING_PARENTS;" : "" );
        std::vector< int > smetadata, tmetadata;
        // Load the source mesh and validate
        rval = remapper.LoadNativeMesh( ctx.inFilenames[0], ctx.meshsets[0], smetadata, additional_read_opts );MB_CHK_ERR( rval );
        if( smetadata.size() )
        {
            // remapper.SetMeshType( moab::Remapper::SourceMesh,
            //                       static_cast< moab::TempestRemapper::TempestMeshType >( smetadata[0] ), &smetadata
            //                       );
            remapper.SetMeshType( moab::Remapper::SourceMesh,
                                  static_cast< moab::TempestRemapper::TempestMeshType >( smetadata[0] ) );
        }
        // Rescale the radius of both to compute the intersection
        rval = moab::IntxUtils::ScaleToRadius( ctx.mbcore, ctx.meshsets[0], radius_src );MB_CHK_ERR( rval );
        rval = remapper.ConvertMeshToTempest( moab::Remapper::SourceMesh );MB_CHK_ERR( rval );
        ctx.meshes[0] = remapper.GetMesh( moab::Remapper::SourceMesh );

        // Load the target mesh and validate
        rval = remapper.LoadNativeMesh( ctx.inFilenames[1], ctx.meshsets[1], tmetadata, additional_read_opts );MB_CHK_ERR( rval );
        if( tmetadata.size() )
        {
            // remapper.SetMeshType( moab::Remapper::TargetMesh,
            //                       static_cast< moab::TempestRemapper::TempestMeshType >( tmetadata[0] ), &tmetadata
            //                       );
            remapper.SetMeshType( moab::Remapper::TargetMesh,
                                  static_cast< moab::TempestRemapper::TempestMeshType >( tmetadata[0] ) );
        }
        rval = moab::IntxUtils::ScaleToRadius( ctx.mbcore, ctx.meshsets[1], radius_dest );MB_CHK_ERR( rval );
        rval = remapper.ConvertMeshToTempest( moab::Remapper::TargetMesh );MB_CHK_ERR( rval );
        ctx.meshes[1] = remapper.GetMesh( moab::Remapper::TargetMesh );
    }
    else if( ctx.meshType == moab::TempestRemapper::ICO )
    {
        err = GenerateICOMesh( *tempest_mesh, ctx.blockSize, ctx.computeDual, ctx.outFilename, "NetCDF4" );

        if( err )
        {
            rval = moab::MB_FAILURE;
        }
        else
        {
            ctx.meshes.push_back( tempest_mesh );
        }
    }
    else if( ctx.meshType == moab::TempestRemapper::RLL )
    {
        err = GenerateRLLMesh( *tempest_mesh,                     // Mesh& meshOut,
                               ctx.blockSize * 2, ctx.blockSize,  // int nLongitudes, int nLatitudes,
                               0.0, 360.0,                        // double dLonBegin, double dLonEnd,
                               -90.0, 90.0,                       // double dLatBegin, double dLatEnd,
                               false, false, false,  // bool fGlobalCap, bool fFlipLatLon, bool fForceGlobal,
                               "" /*ctx.inFilename*/, "",
                               "",  // std::string strInputFile, std::string strInputFileLonName, std::string
                                    // strInputFileLatName,
                               ctx.outFilename, "NetCDF4",  // std::string strOutputFile, std::string strOutputFormat
                               true                         // bool fVerbose
        );

        if( err )
        {
            rval = moab::MB_FAILURE;
        }
        else
        {
            ctx.meshes.push_back( tempest_mesh );
        }
    }
    else  // default
    {
        err = GenerateCSMesh( *tempest_mesh, ctx.blockSize, ctx.outFilename, "NetCDF4" );

        if( err )
        {
            rval = moab::MB_FAILURE;
        }
        else
        {
            ctx.meshes.push_back( tempest_mesh );
        }
    }

    if( ctx.meshType != moab::TempestRemapper::OVERLAP_MOAB && !tempest_mesh )
    {
        std::cout << "Tempest Mesh is not a complete object; Quitting...";
        exit( -1 );
    }

    return rval;
}

///////////////////////////////////////////////
//         Test functions

double sample_slow_harmonic( double dLon, double dLat )<--- The function 'sample_slow_harmonic' is never used.
{
    return ( 2.0 + cos( dLat ) * cos( dLat ) * cos( 2.0 * dLon ) );
}

double sample_fast_harmonic( double dLon, double dLat )<--- The function 'sample_fast_harmonic' is never used.
{
    return ( 2.0 + pow( sin( 2.0 * dLat ), 16.0 ) * cos( 16.0 * dLon ) );
    // return (2.0 + pow(cos(2.0 * dLat), 16.0) * cos(16.0 * dLon));
}

double sample_constant( double /*dLon*/, double /*dLat*/ )<--- The function 'sample_constant' is never used.
{
    return 1.0;
}

double sample_stationary_vortex( double dLon, double dLat )
{
    const double dLon0 = 0.0;
    const double dLat0 = 0.6;
    const double dR0   = 3.0;
    const double dD    = 5.0;
    const double dT    = 6.0;

    ///		Find the rotated longitude and latitude of a point on a sphere
    ///		with pole at (dLonC, dLatC).
    {
        double dSinC = sin( dLat0 );
        double dCosC = cos( dLat0 );
        double dCosT = cos( dLat );
        double dSinT = sin( dLat );

        double dTrm = dCosT * cos( dLon - dLon0 );
        double dX   = dSinC * dTrm - dCosC * dSinT;
        double dY   = dCosT * sin( dLon - dLon0 );
        double dZ   = dSinC * dSinT + dCosC * dTrm;

        dLon = atan2( dY, dX );
        if( dLon < 0.0 )
        {
            dLon += 2.0 * M_PI;
        }
        dLat = asin( dZ );
    }

    double dRho = dR0 * cos( dLat );
    double dVt  = 3.0 * sqrt( 3.0 ) / 2.0 / cosh( dRho ) / cosh( dRho ) * tanh( dRho );

    double dOmega;
    if( dRho == 0.0 )
    {
        dOmega = 0.0;
    }
    else
    {
        dOmega = dVt / dRho;
    }

    return ( 1.0 - tanh( dRho / dD * sin( dLon - dOmega * dT ) ) );
}

///////////////////////////////////////////////