MOAB: Mesh Oriented datABase  (version 5.2.1)
cslam_par_test.cpp File Reference
#include <iostream>
#include <sstream>
#include <time.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include "moab/Core.hpp"
#include "moab/Interface.hpp"
#include "moab/IntxMesh/Intx2MeshOnSphere.hpp"
#include <math.h>
#include "TestUtil.hpp"
#include "moab/ParallelComm.hpp"
#include "moab/ProgOptions.hpp"
#include "MBParallelConventions.h"
#include "moab/ReadUtilIface.hpp"
#include "MBTagConventions.hpp"
#include "moab/IntxMesh/IntxUtils.hpp"
#include "IntxUtilsCSLAM.hpp"
+ Include dependency graph for cslam_par_test.cpp:

Go to the source code of this file.

Functions

std::string input_mesh_file ("VELO00_16p.h5m")
void test_intx_in_parallel_elem_based ()
int main (int argc, char **argv)
ErrorCode compute_lagrange_mesh_on_sphere (Interface *mb, EntityHandle euler_set)

Variables

double EPS1 = 0.2
double Radius = 1.0
double deltaT = 1.e-6

Function Documentation

ErrorCode compute_lagrange_mesh_on_sphere ( Interface *  mb,
EntityHandle  euler_set 
)

Definition at line 75 of file cslam_par_test.cpp.

References moab::Range::begin(), deltaT, moab::Range::end(), moab::CartVect::length(), MB_CHK_ERR, MB_TAG_CREAT, MB_TAG_DENSE, MB_TYPE_DOUBLE, MBQUAD, Radius, and moab::Range::size().

Referenced by test_intx_in_parallel_elem_based().

{
    /*
     * get all quads first, then vertices, then move them on the surface of the sphere
     *  radius is 1, usually
     *  pos (t-dt) = pos(t) -Velo(t)*dt; this will be lagrange mesh, on each processor
     */
    Range quads;
    ErrorCode rval = mb->get_entities_by_type( euler_set, MBQUAD, quads );MB_CHK_ERR( rval );

    Range connecVerts;
    rval = mb->get_connectivity( quads, connecVerts );MB_CHK_ERR( rval );

    // the LOC tag, should be provided by the user?
    Tag tagh = 0;
    std::string tag_name( "DP" );
    rval = mb->tag_get_handle( tag_name.c_str(), 3, MB_TYPE_DOUBLE, tagh, MB_TAG_DENSE | MB_TAG_CREAT );MB_CHK_ERR( rval );
    void* data;  // pointer to the DP in memory, for each vertex
    int count;

    rval = mb->tag_iterate( tagh, connecVerts.begin(), connecVerts.end(), count, data );MB_CHK_ERR( rval );
    // here we are checking contiguity
    assert( count == (int)connecVerts.size() );
    double* ptr_DP = (double*)data;
    // get the coordinates of the old mesh, and move it around using velocity tag

    Tag tagv = 0;
    std::string velo_tag_name( "VELO" );
    rval = mb->tag_get_handle( velo_tag_name.c_str(), 3, MB_TYPE_DOUBLE, tagv, MB_TAG_DENSE );MB_CHK_ERR( rval );

    /*void *datavelo; // pointer to the VELO in memory, for each vertex

    rval = mb->tag_iterate(tagv, connecVerts.begin(), connecVerts.end(), count, datavelo);MB_CHK_ERR(rval);*/
    // here we are checking contiguity
    assert( count == (int)connecVerts.size() );
    // now put the vertices in the right place....
    // int vix=0; // vertex index in new array

    for( Range::iterator vit = connecVerts.begin(); vit != connecVerts.end(); ++vit )
    {
        EntityHandle oldV = *vit;
        CartVect posi;
        rval = mb->get_coords( &oldV, 1, &( posi[0] ) );MB_CHK_ERR( rval );
        CartVect velo;
        rval = mb->tag_get_data( tagv, &oldV, 1, (void*)&( velo[0] ) );MB_CHK_ERR( rval );
        // do some mumbo jumbo, as in python script
        CartVect newPos = posi - deltaT * velo;
        double len1     = newPos.length();
        newPos          = Radius * newPos / len1;

        ptr_DP[0] = newPos[0];
        ptr_DP[1] = newPos[1];
        ptr_DP[2] = newPos[2];
        ptr_DP += 3;  // increment to the next node
    }

    return rval;
}
std::string input_mesh_file ( "VELO00_16p.h5m"  )
int main ( int  argc,
char **  argv 
)

Definition at line 47 of file cslam_par_test.cpp.

References deltaT, EPS1, input_mesh_file(), Radius, RUN_TEST, and test_intx_in_parallel_elem_based().

{
    MPI_Init( &argc, &argv );
    EPS1       = 0.000002;
    int result = 0;

    if( argc > 1 )
    {
        int index = 1;
        while( index < argc )
        {
            if( !strcmp( argv[index], "-eps" ) )  // this is for box error
            { EPS1 = atof( argv[++index] ); }
            if( !strcmp( argv[index], "-input" ) ) { input_mesh_file = argv[++index]; }
            if( !strcmp( argv[index], "-radius" ) ) { Radius = atof( argv[++index] ); }
            if( !strcmp( argv[index], "-deltaT" ) ) { deltaT = atof( argv[++index] ); }
            index++;
        }
    }
    std::cout << " run: -input " << input_mesh_file << "  -eps " << EPS1 << " -radius " << Radius << " -deltaT "
              << deltaT << "\n";

    result += RUN_TEST( test_intx_in_parallel_elem_based );

    MPI_Finalize();
    return result;
}

Definition at line 134 of file cslam_par_test.cpp.

References moab::IntxAreaUtils::area_on_sphere(), moab::ParallelComm::check_all_shared_handles(), compute_lagrange_mesh_on_sphere(), moab::Intx2Mesh::create_departure_mesh_2nd_alg(), moab::Core::create_meshset(), EPS1, example, moab::Intx2Mesh::FindMaxEdges(), moab::ParallelComm::get_pcomm(), input_mesh_file(), moab::Intx2Mesh::intersect_meshes(), mb, MB_CHK_ERR_RET, MESHSET_SET, moab::ParallelComm::proc_config(), moab::ProcConfig::proc_rank(), Radius, rank, moab::Intx2Mesh::set_box_error(), moab::Intx2Mesh::set_error_tolerance(), moab::Intx2MeshOnSphere::set_radius_destination_mesh(), and moab::Intx2MeshOnSphere::set_radius_source_mesh().

Referenced by main().

{
    std::string opts = std::string( "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION" ) +
                       std::string( ";PARALLEL_RESOLVE_SHARED_ENTS" );
    Core moab;
    Interface& mb = moab;
    EntityHandle euler_set;
    ErrorCode rval;
    rval = mb.create_meshset( MESHSET_SET, euler_set );MB_CHK_ERR_RET( rval );
    std::string example( TestDir + "/" + input_mesh_file );

    rval = mb.load_file( example.c_str(), &euler_set, opts.c_str() );

    ParallelComm* pcomm = ParallelComm::get_pcomm( &mb, 0 );MB_CHK_ERR_RET( rval );

    rval = pcomm->check_all_shared_handles();MB_CHK_ERR_RET( rval );

    // everybody will get a DP tag, including the non owned entities; so exchange tags is not
    // required for LOC (here)
    rval = compute_lagrange_mesh_on_sphere( &mb, euler_set );MB_CHK_ERR_RET( rval );

    int rank = pcomm->proc_config().proc_rank();

    std::stringstream ste;
    ste << "initial" << rank << ".vtk";
    mb.write_file( ste.str().c_str(), 0, 0, &euler_set, 1 );

    Intx2MeshOnSphere worker( &mb );

    worker.set_radius_source_mesh( Radius );
    worker.set_radius_destination_mesh( Radius );
    worker.set_box_error( EPS1 );  //
    // worker.SetEntityType(MBQUAD);

    worker.set_error_tolerance( Radius * 1.e-8 );
    std::cout << "error tolerance epsilon_1=" << Radius * 1.e-8 << "\n";
    //  worker.locate_departure_points(euler_set);

    rval = worker.FindMaxEdges( euler_set, euler_set );  // departure will be the same max_edges
    // we need to make sure the covering set is bigger than the euler mesh
    EntityHandle covering_lagr_set;
    rval = mb.create_meshset( MESHSET_SET, covering_lagr_set );MB_CHK_ERR_RET( rval );

    rval = worker.create_departure_mesh_2nd_alg( euler_set, covering_lagr_set );MB_CHK_ERR_RET( rval );

    std::stringstream ss;
    ss << "partial" << rank << ".vtk";
    mb.write_file( ss.str().c_str(), 0, 0, &covering_lagr_set, 1 );
    EntityHandle outputSet;
    rval = mb.create_meshset( MESHSET_SET, outputSet );MB_CHK_ERR_RET( rval );
    rval = worker.intersect_meshes( covering_lagr_set, euler_set, outputSet );MB_CHK_ERR_RET( rval );

    // std::string opts_write("PARALLEL=WRITE_PART");
    // rval = mb.write_file("manuf.h5m", 0, opts_write.c_str(), &outputSet, 1);
    // std::string opts_write("");
    std::stringstream outf;
    outf << "intersect" << rank << ".h5m";
    rval = mb.write_file( outf.str().c_str(), 0, 0, &outputSet, 1 );MB_CHK_ERR_RET( rval );

    moab::IntxAreaUtils sphAreaUtils;
    double intx_area    = sphAreaUtils.area_on_sphere( &mb, outputSet, Radius );
    double arrival_area = sphAreaUtils.area_on_sphere( &mb, euler_set, Radius );
    std::cout << "On rank : " << rank << " arrival area: " << arrival_area << "  intersection area:" << intx_area
              << " rel error: " << fabs( ( intx_area - arrival_area ) / arrival_area ) << "\n";
}

Variable Documentation

double deltaT = 1.e-6

Definition at line 44 of file cslam_par_test.cpp.

Referenced by compute_lagrange_mesh_on_sphere(), and main().

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines