test_remapping.cpp

Go to the documentation of this file.

/**
 * MOAB, a Mesh-Oriented datABase, is a software component for creating,
 * storing and accessing finite element mesh data.
 *
 * Copyright 2004 Sandia Corporation.  Under the terms of Contract
 * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
 * retains certain rights in this software.
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 *
 */

#ifdef WIN32
#ifdef _DEBUG
// turn off warnings that say they debugging identifier has been truncated
// this warning comes up when using some STL containers
#pragma warning( disable : 4786 )
#endif
#endif

#include <iostream>
#include "moab/Core.hpp"
#include "moab/Remapping/TempestRemapper.hpp"

#ifdef MOAB_HAVE_MPI
#include "moab/ParallelComm.hpp"
#endif

#ifndef IS_BUILDING_MB
#define IS_BUILDING_MB
#endif

#include "Internals.hpp"
#include "TestRunner.hpp"

using namespace moab;

const static double radius               = 1.0;
const double MOAB_PI                     = 3.1415926535897932384626433832795028841971693993751058209749445923;
const static double surface_area         = 4.0 * MOAB_PI * radius * radius;
const static std::string outFilenames[5] = { "outTempestCS.g", "outTempestRLL.g", "outTempestICO.g", "outTempestICOD.g",
                                             "outTempestOV.g" };

void test_tempest_cs_create();
void test_tempest_rll_create();
void test_tempest_ico_create();
void test_tempest_mpas_create();
void test_tempest_overlap_combinations();
void test_tempest_to_moab_convert();

int main( int argc, char** argv )
{
    REGISTER_TEST( test_tempest_cs_create );
    REGISTER_TEST( test_tempest_rll_create );
    REGISTER_TEST( test_tempest_ico_create );
    REGISTER_TEST( test_tempest_mpas_create );
    REGISTER_TEST( test_tempest_overlap_combinations );
    REGISTER_TEST( test_tempest_to_moab_convert );

#ifdef MOAB_HAVE_MPI
    MPI_Init( &argc, &argv );
#endif
    int result = RUN_TESTS( argc, argv );
#ifdef MOAB_HAVE_MPI
    MPI_Finalize();
#endif
    return result;
}

void test_tempest_cs_create()
{
    NcError error( NcError::verbose_nonfatal );
    const int blockSize           = 30;
    const std::string outFilename = outFilenames[0];

    std::cout << "Creating TempestRemap Cubed-Sphere Mesh ...\n";
    Mesh tempest_mesh;
    int ierr = GenerateCSMesh( tempest_mesh, blockSize, outFilename, "NetCDF4" );
    CHECK_EQUAL( ierr, 0 );

    // Compute the surface area of CS mesh
    const double sphere_area = tempest_mesh.CalculateFaceAreas( false );
    CHECK_REAL_EQUAL( sphere_area, surface_area, 1e-10 );
}

void test_tempest_rll_create()
{
    NcError error( NcError::verbose_nonfatal );
    const int blockSize           = 30;
    const std::string outFilename = outFilenames[1];

    std::cout << "Creating TempestRemap Latitude-Longitude Mesh ...\n";
    Mesh tempest_mesh;
    int ierr =
        GenerateRLLMesh( tempest_mesh, blockSize * 2, blockSize, 0.0, 360.0, -90.0, 90.0, false, false, true, "", "",
                         "",  // std::string strInputFile, std::string strInputFileLonName, std::string
                              // strInputFileLatName,
                         outFilename, "NetCDF4",  // std::string strOutputFile, std::string strOutputFormat
                         false );
    CHECK_EQUAL( ierr, 0 );

    // Compute the surface area of RLL mesh
    const double sphere_area = tempest_mesh.CalculateFaceAreas( false );
    CHECK_REAL_EQUAL( sphere_area, surface_area, 1e-10 );
}

void test_tempest_ico_create()
{
    NcError error( NcError::verbose_nonfatal );
    const int blockSize           = 30;
    const bool computeDual        = false;
    const std::string outFilename = outFilenames[2];

    std::cout << "Creating TempestRemap Icosahedral Mesh ...\n";
    Mesh tempest_mesh;
    int ierr = GenerateICOMesh( tempest_mesh, blockSize, computeDual, outFilename, "NetCDF4" );
    CHECK_EQUAL( ierr, 0 );

    // Compute the surface area of ICO mesh
    const double sphere_area = tempest_mesh.CalculateFaceAreas( false );
    CHECK_REAL_EQUAL( sphere_area, surface_area, 1e-10 );
}

void test_tempest_mpas_create()
{
    NcError error( NcError::verbose_nonfatal );
    const int blockSize           = 30;
    const bool computeDual        = true;
    const std::string outFilename = outFilenames[3];

    std::cout << "Creating TempestRemap MPAS Mesh (dual of the Icosahedral) ...\n";
    Mesh tempest_mesh;
    int ierr = GenerateICOMesh( tempest_mesh, blockSize, computeDual, outFilename, "NetCDF4" );
    CHECK_EQUAL( ierr, 0 );

    // Compute the surface area of MPAS mesh
    const double sphere_area = tempest_mesh.CalculateFaceAreas( false );
    CHECK_REAL_EQUAL( sphere_area, surface_area, 1e-10 );
}

void test_tempest_overlap_combinations()
{
    NcError error( NcError::verbose_nonfatal );
    const std::string outFilename = outFilenames[4];

    Mesh inpMesh( outFilenames[0] );
    // verify input mesh area first
    const double inpArea = inpMesh.CalculateFaceAreas( false );
    CHECK_REAL_EQUAL( inpArea, surface_area, 1e-10 );

    for( int isrc = 0; isrc < 4; ++isrc )
    {
        for( int jsrc = 0; jsrc < 4; ++jsrc )
        {
            std::cout << "Computing Overlap between " << outFilenames[isrc] << " and " << outFilenames[jsrc]
                      << " ...\n";
            Mesh tempest_mesh;
            int ierr = GenerateOverlapMesh( outFilenames[isrc], outFilenames[jsrc], tempest_mesh, outFilename,
                                            "NetCDF4", "exact", false, false, false, false, false );
            CHECK_EQUAL( ierr, 0 );
            // verify overlap mesh area
            const double ovArea = tempest_mesh.CalculateFaceAreas( false );
            CHECK_REAL_EQUAL( ovArea, surface_area, 1e-10 );
        }
    }
}

void test_tempest_to_moab_convert()
{
    NcError error( NcError::verbose_nonfatal );

    // Allocate and create MOAB Remapper object
    moab::ErrorCode rval;
    moab::Interface* mbCore = new( std::nothrow ) moab::Core;
    CHECK( NULL != mbCore );

#ifdef MOAB_HAVE_MPI
    moab::ParallelComm* pcomm       = new moab::ParallelComm( mbCore, MPI_COMM_WORLD, 0 );
    moab::TempestRemapper* remapper = new moab::TempestRemapper( mbCore, pcomm );
#else
    moab::TempestRemapper* remapper = new moab::TempestRemapper( mbCore );
#endif
    remapper->meshValidate     = true;
    remapper->constructEdgeMap = true;
    remapper->initialize();

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

    rval = remapper->LoadMesh( moab::Remapper::SourceMesh, outFilenames[0], moab::TempestRemapper::CS );CHECK_ERR( rval );

    // Load the meshes and validate
    rval = remapper->ConvertTempestMesh( moab::Remapper::SourceMesh );CHECK_ERR( rval );

    Mesh* srcTempest = remapper->GetMesh( moab::Remapper::SourceMesh );

    moab::EntityHandle srcset = remapper->GetMeshSet( moab::Remapper::SourceMesh );

    moab::EntityHandle& tgtset = remapper->GetMeshSet( moab::Remapper::TargetMesh );

    tgtset = srcset;

    // Load the meshes and validate
    rval = remapper->ConvertMeshToTempest( moab::Remapper::TargetMesh );CHECK_ERR( rval );

    Mesh* tgtTempest = remapper->GetMesh( moab::Remapper::TargetMesh );

    const size_t tempest_nodes_src = srcTempest->nodes.size(), tempest_elems_src = srcTempest->faces.size();
    const size_t tempest_nodes_tgt = tgtTempest->nodes.size(), tempest_elems_tgt = tgtTempest->faces.size();
    CHECK_EQUAL( tempest_nodes_src, tempest_nodes_tgt );
    CHECK_EQUAL( tempest_elems_src, tempest_elems_tgt );

    delete remapper;
#ifdef MOAB_HAVE_MPI
    delete pcomm;
#endif
    delete mbCore;
}