hireconst_test_parallel.cpp

Go to the documentation of this file.

/*This unit test is for high order reconstruction capability, which could be used for mesh
 * refinement*/
#include <iostream>
#include <string>
#include <sstream>
#if defined( __MINGW32__ )
#include <sys/time.h>
#else
#include <ctime>
#endif

#include <vector>
#include <algorithm>
#include "moab/Core.hpp"
#include "moab/Range.hpp"
#include "moab/CartVect.hpp"
#include "moab/MeshTopoUtil.hpp"
#include "moab/NestedRefine.hpp"
#include "moab/DiscreteGeometry/HiReconstruction.hpp"
#include "TestUtil.hpp"
#include <cmath>

#ifdef MOAB_HAVE_MPI
#include "moab/ParallelComm.hpp"
#include "MBParallelConventions.h"
#include "ReadParallel.hpp"
#include "moab/FileOptions.hpp"
#include "MBTagConventions.hpp"
#include "moab_mpi.h"
#endif

using namespace moab;

#ifdef MOAB_HAVE_MPI
std::string read_options;
#endif

#ifdef MOAB_HAVE_HDF5
#undef MOAB_HAVE_HDF5
#endif
ErrorCode load_meshset_hirec( const char* infile,
                              Interface* mbimpl,
                              EntityHandle& meshset,
                              ParallelComm*& pc,
                              const int degree = 0,
                              const int dim    = 2 );
ErrorCode test_mesh( const char* infile, const int degree, const bool interp, const int dim );

void compute_linear_coords( const int nvpe, double* elemcoords, double* naturals, double* linearcoords );

void usage()
{
    std::cout << "usage: mpirun -np <number of processors> ./hireconst_test_parallel <mesh file> "
                 "-degree <degree> -interp <0=least square, 1=interpolation> -dim <mesh dimension>"
              << std::endl;
}

int main( int argc, char* argv[] )
{
#ifdef MOAB_HAVE_MPI
    MPI_Init( &argc, &argv );
    int nprocs, rank;
    MPI_Comm_size( MPI_COMM_WORLD, &nprocs );
    MPI_Comm_rank( MPI_COMM_WORLD, &rank );
#endif
    int degree = 3, dim = 2;
    bool interp = false;
    ErrorCode error;

#ifdef MOAB_HAVE_HDF5
    std::string infile = TestDir + "unittest/mbcslam/fine4.h5m";
#else
    std::string infile = TestDir + "unittest/sphere_quads_20.vtk";
#endif

    if( argc == 1 )
    {
        usage();
        std::cout << "Using default arguments: ./hireconst_test_parallel " << infile << " -degree 3 -interp 0 -dim 2"
                  << std::endl;
    }
    else
    {
        infile      = argv[1];
        bool hasdim = false;

        for( int i = 2; i < argc; ++i )
        {
            if( i + 1 != argc )
            {
                if( std::string( argv[i] ) == "-degree" )
                {
                    degree = atoi( argv[++i] );
                }
                else if( std::string( argv[i] ) == "-interp" )
                {
                    interp = atoi( argv[++i] );
                }
                else if( std::string( argv[i] ) == "-dim" )
                {
                    dim    = atoi( argv[++i] );
                    hasdim = true;
                }
                else
                {
#ifdef MOAB_HAVE_MPI

                    if( 0 == rank )
                    {
                        usage();
                    }

                    MPI_Finalize();
#else
                    usage();
#endif
                    return 0;
                }
            }
        }

        if( !hasdim )
        {
#ifdef MOAB_HAVE_MPI

            if( 0 == rank )
            {
                std::cout << "Dimension of input mesh should be provided, positive and less than 3" << std::endl;
            }

#else
            std::cout << "Dimension of input mesh should be provided, positive and less than 3" << std::endl;
#endif
            return 0;
        }

        if( degree <= 0 || dim > 2 || dim <= 0 )
        {
#ifdef MOAB_HAVE_MPI

            if( 0 == rank )
            {
                std::cout << "Input degree should be positive number;\n";
                std::cout << "Input dimesion should be positive and less than 3;" << std::endl;
            }

#else
            std::cout << "Input degree should be positive number;\n";
            std::cout << "Input dimesion should be positive and less than 3;" << std::endl;
#endif
            return 0;
        }

#ifdef MOAB_HAVE_MPI

        if( 0 == rank )
        {
            std::cout << "Testing on " << infile << " with dimension " << dim << "\n";
            std::string opts = interp ? "interpolation" : "least square fitting";
            std::cout << "High order reconstruction with degree " << degree << " " << opts << std::endl;
        }

#else
        std::cout << "Testing on " << infile << " with dimension " << dim << "\n";
        std::string opts = interp ? "interpolation" : "least square fitting";
        std::cout << "High order reconstruction with degree " << degree << " " << opts << std::endl;
#endif
    }

    error = test_mesh( infile.c_str(), degree, interp, dim );MB_CHK_ERR( error );

#ifdef MOAB_HAVE_MPI
    MPI_Finalize();
#endif
}

ErrorCode load_meshset_hirec( const char* infile,
                              Interface* mbimpl,
                              EntityHandle& meshset,
                              ParallelComm*& pc,
                              const int degree,
                              const int dim )
{
    ErrorCode error;
    error = mbimpl->create_meshset( moab::MESHSET_SET, meshset );MB_CHK_ERR( error );
#ifdef MOAB_HAVE_MPI
    int nprocs, rank;
    MPI_Comm comm = MPI_COMM_WORLD;
    MPI_Comm_size( comm, &nprocs );
    MPI_Comm_rank( comm, &rank );
    EntityHandle partnset;
    error = mbimpl->create_meshset( moab::MESHSET_SET, partnset );MB_CHK_ERR( error );

    if( nprocs > 1 )
    {
        pc = moab::ParallelComm::get_pcomm( mbimpl, partnset, &comm );
    }

    if( nprocs > 1 )
    {
        int nghlayers           = degree > 0 ? HiReconstruction::estimate_num_ghost_layers( degree, true ) : 0;
        std::string part_method = "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;";
#ifndef MOAB_HAVE_HDF5
        part_method = "PARALLEL=BCAST_DELETE;PARTITION=TRIVIAL;";
#endif

        if( nghlayers )
        {
            // get ghost layers
            if( dim == 2 )
            {
                read_options = part_method + ";PARALLEL_RESOLVE_SHARED_ENTS;PARALLEL_GHOSTS=2.0.";
            }
            else if( dim == 1 )
            {
                read_options = part_method + ";PARALLEL_RESOLVE_SHARED_ENTS;PARALLEL_GHOSTS=1.0.";
            }
            else
            {
                MB_SET_ERR( MB_FAILURE, "unsupported dimension" );
            }

            read_options += (char)( '0' + nghlayers );
        }
        else
        {
            read_options = part_method + ";PARALLEL_RESOLVE_SHARED_ENTS;";
        }

        error = mbimpl->load_file( infile, &meshset, read_options.c_str() );MB_CHK_ERR( error );
    }
    else
    {
        error = mbimpl->load_file( infile, &meshset );MB_CHK_ERR( error );
    }

#else
    assert( !pc && degree && dim );
    error = mbimpl->load_file( infile, &meshset );MB_CHK_ERR( error );
#endif
    return error;
}

ErrorCode test_mesh( const char* infile, const int degree, const bool interp, const int dim )
{
    Core moab;
    Interface* mbimpl = &moab;
    ParallelComm* pc  = NULL;
    EntityHandle meshset;
#ifdef MOAB_HAVE_MPI
    int nprocs, rank;
    MPI_Comm comm = MPI_COMM_WORLD;
    MPI_Comm_size( comm, &nprocs );
    MPI_Comm_rank( comm, &rank );
#endif

    ErrorCode error;
    // mesh will be loaded and communicator pc will be updated
    error = load_meshset_hirec( infile, mbimpl, meshset, pc, degree, dim );MB_CHK_ERR( error );
    // initialize
    HiReconstruction hirec( dynamic_cast< Core* >( mbimpl ), pc, meshset );
    Range elems, elems_owned;
    error = mbimpl->get_entities_by_dimension( meshset, dim, elems );MB_CHK_ERR( error );
    int nelems = elems.size();

#ifdef MOAB_HAVE_MPI

    if( pc )
    {
        error = pc->filter_pstatus( elems, PSTATUS_GHOST, PSTATUS_NOT, -1, &elems_owned );MB_CHK_ERR( error );
    }
    else
    {
        elems_owned = elems;
    }

#endif

#ifdef MOAB_HAVE_MPI
    std::cout << "Mesh has " << nelems << " elements on Processor " << rank << " in total;";
    std::cout << elems_owned.size() << " of which are locally owned elements" << std::endl;
#else
    std::cout << "Mesh has " << nelems << " elements" << std::endl;
#endif

    // reconstruction
    if( dim == 2 )
    {
        error = hirec.reconstruct3D_surf_geom( degree, interp, false );MB_CHK_ERR( error );
    }
    else if( dim == 1 )
    {
        error = hirec.reconstruct3D_curve_geom( degree, interp, false );MB_CHK_ERR( error );
    }

#ifdef MOAB_HAVE_MPI
    std::cout << "HiRec has been done on Processor " << rank << std::endl;
#else
    std::cout << "HiRec has been done " << std::endl;
#endif
    // fitting
    double mxdist = 0;

    for( Range::iterator ielem = elems_owned.begin(); ielem != elems_owned.end(); ++ielem )
    {
        int nvpe;
        const EntityHandle* conn;
        error = mbimpl->get_connectivity( *ielem, conn, nvpe );MB_CHK_ERR( error );
        double w = 1.0 / (double)nvpe;
        std::vector< double > naturalcoords2fit( nvpe, w );
        CartVect newcoords, linearcoords;
        error = hirec.hiproj_walf_in_element( *ielem, nvpe, 1, &( naturalcoords2fit[0] ), newcoords.array() );

        if( MB_FAILURE == error )
        {
            continue;
        }

        std::vector< double > coords( 3 * nvpe );
        error = mbimpl->get_coords( conn, nvpe, &( coords[0] ) );MB_CHK_ERR( error );
        compute_linear_coords( nvpe, &( coords[0] ), &( naturalcoords2fit[0] ), linearcoords.array() );
        CartVect nlcoords = newcoords - linearcoords;
        mxdist            = std::max( mxdist, nlcoords.length() );
        /*#ifdef MOAB_HAVE_MPI
            std::cout << "Error on element " << *ielem << " is " << nlcoords.length() << "on
        Processor " << rank << std::endl; #else std::cout << "Error on element " << *ielem << " is "
        << nlcoords.length() << std::endl; #endif*/
    }

#ifdef MOAB_HAVE_MPI
    std::cout << "Maximum projection lift is " << mxdist << " on Processor " << rank << std::endl;
#else
    std::cout << "Maximum projection lift is " << mxdist << std::endl;
#endif
    return error;
}

void compute_linear_coords( const int nvpe, double* elemcoords, double* naturals, double* linearcoords )
{
    assert( elemcoords && linearcoords );

    for( int i = 0; i < 3; ++i )
    {
        linearcoords[i] = 0;

        for( int j = 0; j < nvpe; ++j )
        {
            linearcoords[i] += naturals[j] * elemcoords[3 * j + i];
        }
    }
}