intx_imesh.cpp File Reference

#include <string>
#include <iostream>
#include "moab/Core.hpp"
#include "moab/IntxMesh/Intx2MeshOnSphere.hpp"
#include "moab/IntxMesh/IntxUtils.hpp"
#include "moab_mpi.h"
#include "iMeshP.h"
#include "MBiMesh.hpp"

Include dependency graph for intx_imesh.cpp:

Go to the source code of this file.

Defines
#define	IMESH_ASSERT(ierr)   if( ( ierr ) != 0 ) std::cout << "iMesh Assert: \n";
#define	IMESH_NULL   0
Functions
void	update_tracer (iMesh_Instance instance, iBase_EntitySetHandle imesh_euler_set, int *ierr)
int	main (int argc, char *argv[])

---

Define Documentation

#define IMESH_ASSERT

(

ierr

)

if( ( ierr ) != 0 ) std::cout << "iMesh Assert: \n";

Definition at line 16 of file intx_imesh.cpp.

Referenced by main().

#define IMESH_NULL   0

Definition at line 18 of file intx_imesh.cpp.

Referenced by main().

---

Function Documentation

int main	(	int	argc,
		char *	argv[]
	)

Definition at line 22 of file intx_imesh.cpp.

References filename, iBase_FACE, ierr, iMesh_addEntArrToSet, iMesh_ALL_TOPOLOGIES, IMESH_ASSERT, iMesh_createEntSet, iMesh_getEntities, iMesh_getNumEntSets, iMesh_getRootSet, iMesh_newMesh, IMESH_NULL, iMeshP_createPartitionAll, iMeshP_loadAll, iMeshP_saveAll, MPI_COMM_WORLD, rank, size, and update_tracer().

{
    MPI_Init( &argc, &argv );

    iMesh_Instance imesh;
    iMeshP_PartitionHandle partn;
    int ierr, num_sets;

    iBase_EntitySetHandle root;
    imesh = IMESH_NULL;
    iMesh_newMesh( 0, &imesh, &ierr, 0 );
    IMESH_ASSERT( ierr );
    iMesh_getRootSet( imesh, &root, &ierr );
    IMESH_ASSERT( ierr );

    iMeshP_createPartitionAll( imesh, MPI_COMM_WORLD, &partn, &ierr );
    int rank, size;
    MPI_Comm_rank( MPI_COMM_WORLD, &rank );
    MPI_Comm_size( MPI_COMM_WORLD, &size );
    IMESH_ASSERT( ierr );

    const char options[] = " moab:PARALLEL=READ_PART "
                           " moab:PARTITION=PARALLEL_PARTITION "
                           " moab:PARALLEL_RESOLVE_SHARED_ENTS "
                           " moab:PARTITION_DISTRIBUTE ";
    const char* filename = "HN16DP.h5m";  // the file should have the dp tag already

    if( 0 == rank ) std::cout << "Load in parallel the file: " << filename << std::endl;
    iMeshP_loadAll( imesh, partn, root, filename, options, &ierr, strlen( filename ), strlen( options ) );
    IMESH_ASSERT( ierr );

    iMesh_getNumEntSets( imesh, IMESH_NULL, 1, &num_sets, &ierr );
    IMESH_ASSERT( ierr );
    std::cout << "There's " << num_sets << " entity sets here on process rank " << rank << std::endl;

    iBase_EntitySetHandle euler_set;

    iMesh_createEntSet( imesh, 0, &euler_set, &ierr );
    IMESH_ASSERT( ierr );

    iBase_EntityHandle* cells = NULL;
    int ents_alloc            = 0;
    int ents_size             = 0;

    iMesh_getEntities( imesh, root, iBase_FACE, iMesh_ALL_TOPOLOGIES, &cells, &ents_alloc, &ents_size, &ierr );
    IMESH_ASSERT( ierr );

    iMesh_addEntArrToSet( imesh, cells, ents_size, euler_set, &ierr );
    IMESH_ASSERT( ierr );

    update_tracer( imesh, euler_set, &ierr );
    IMESH_ASSERT( ierr );

    // write everything
    const char* out_name      = "out.h5m";
    const char optionswrite[] = " moab:PARALLEL=WRITE_PART ";
    iMeshP_saveAll( imesh, partn, euler_set, out_name, optionswrite, &ierr, strlen( out_name ),
                    strlen( optionswrite ) );
    IMESH_ASSERT( ierr );

    if( 0 == rank ) std::cout << "Done\n";
    MPI_Finalize();

    return 0;
}

void update_tracer	(	iMesh_Instance	instance,
		iBase_EntitySetHandle	imesh_euler_set,
		int *	ierr
	)

Definition at line 150 of file wrap_intx.cpp.

References moab::Intx2Mesh::create_departure_mesh_2nd_alg(), moab::Core::create_meshset(), moab::debug, moab::IntxUtils::enforce_convexity(), ErrorCode, ERRORV, gtol, moab::Intx2Mesh::intersect_meshes(), mb, MB_CHK_SET_ERR_RET, MB_TAG_CREAT, MB_TAG_DENSE, MB_TYPE_DOUBLE, MESHSET_SET, MOABI, radius, moab::Intx2Mesh::set_error_tolerance(), moab::Intx2MeshOnSphere::set_radius_destination_mesh(), moab::Intx2MeshOnSphere::set_radius_source_mesh(), moab::Core::tag_get_handle(), moab::Intx2MeshOnSphere::update_tracer_data(), and moab::Core::write_file().

Referenced by main().

{
    using namespace moab;
    const double radius = 1.;
    const double gtol   = 1.e-9;
    const bool debug    = false;

    Range ents;
    moab::Interface* mb = MOABI;
    *ierr               = 1;

    EntityHandle euler_set = (EntityHandle)imesh_euler_set;

    Intx2MeshOnSphere worker( mb );
    worker.SetRadius( radius );

    worker.SetErrorTolerance( gtol );

    EntityHandle covering_lagr_set;

    ErrorCode rval = mb->create_meshset( MESHSET_SET, covering_lagr_set );
    ERRORV( rval, "can't create covering set " );

    // we need to update the correlation tag and remote tuples
    rval = worker.create_departure_mesh_2nd_alg( euler_set, covering_lagr_set );
    ERRORV( rval, "can't populate covering set " );

    if( debug )
    {
        rval = mb->write_file( "lagr.h5m", 0, 0, &covering_lagr_set, 1 );
        ERRORV( rval, "can't write covering set " );
    }

    //
    rval = enforce_convexity( mb, covering_lagr_set );
    ERRORV( rval, "can't write covering set " );

    EntityHandle outputSet;
    rval = mb->create_meshset( MESHSET_SET, outputSet );
    ERRORV( rval, "can't create output set " );

    rval = worker.intersect_meshes( covering_lagr_set, euler_set, outputSet );
    ERRORV( rval, "can't intersect " );

    if( debug )
    {
        rval = mb->write_file( "output.vtk", 0, 0, &outputSet, 1 );
        ERRORV( rval, "can't write covering set " );
    }

    // tagElem is the average computed at each element, from nodal values
    Tag tagElem = 0;
    std::string tag_name2( "TracerAverage" );
    rval = mb->tag_get_handle( tag_name2.c_str(), 1, MB_TYPE_DOUBLE, tagElem, MB_TAG_DENSE | MB_TAG_CREAT );
    ERRORV( rval, "can't get tracer tag " );

    // area of the euler element is fixed, store it; it is used to recompute the averages at each
    // time step
    Tag tagArea = 0;
    std::string tag_name4( "Area" );
    rval = mb->tag_get_handle( tag_name4.c_str(), 1, MB_TYPE_DOUBLE, tagArea, MB_TAG_DENSE | MB_TAG_CREAT );
    ERRORV( rval, "can't get area tag" );

    rval = worker.update_tracer_data( outputSet, tagElem, tagArea );
    ERRORV( rval, "can't update tracer " );

    // everything can be deleted now from intx data; polygons, etc.

    *ierr = 0;
    return;
}