hireconst_test_parallel.cpp File Reference

#include <iostream>
#include <string>
#include <sstream>
#include <ctime>
#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>

Include dependency graph for hireconst_test_parallel.cpp:

Go to the source code of this file.

Functions
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 ()
int	main (int argc, char *argv[])

Function Documentation

compute_linear_coords()

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

Definition at line 338 of file hireconst_test_parallel.cpp.

{
 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];
 }
 }
}

Referenced by test_mesh().

load_meshset_hirec()

ErrorCode load_meshset_hirec	(	const char *	infile,
		Interface *	mbimpl,
		EntityHandle &	meshset,
		ParallelComm *&	pc,
		const int	degree = 0,
		const int	dim = 2
	)

Definition at line 177 of file hireconst_test_parallel.cpp.

{
 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;
}

References moab::Interface::create_meshset(), dim, moab::error(), ErrorCode, moab::HiReconstruction::estimate_num_ghost_layers(), moab::ParallelComm::get_pcomm(), moab::Interface::load_file(), MB_CHK_ERR, MB_SET_ERR, MESHSET_SET, MPI_COMM_WORLD, rank, and read_options.

Referenced by test_mesh().

main()

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

Definition at line 58 of file hireconst_test_parallel.cpp.

{
#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
}

References dim, moab::error(), ErrorCode, MB_CHK_ERR, MPI_COMM_WORLD, rank, test_mesh(), and usage().

test_mesh()

ErrorCode test_mesh	(	const char *	infile,
		const int	degree,
		const bool	interp,
		const int	dim
	)

Definition at line 244 of file hireconst_test_parallel.cpp.

{
 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;
}

References moab::CartVect::array(), moab::Range::begin(), compute_linear_coords(), dim, moab::Range::end(), moab::error(), ErrorCode, moab::ParallelComm::filter_pstatus(), moab::Interface::get_connectivity(), moab::Interface::get_coords(), moab::Interface::get_entities_by_dimension(), moab::HiReconstruction::hiproj_walf_in_element(), moab::CartVect::length(), load_meshset_hirec(), MB_CHK_ERR, MPI_COMM_WORLD, PSTATUS_GHOST, PSTATUS_NOT, rank, moab::HiReconstruction::reconstruct3D_curve_geom(), moab::HiReconstruction::reconstruct3D_surf_geom(), and moab::Range::size().

Referenced by main().

usage()

void usage

(

)

Definition at line 51 of file hireconst_test_parallel.cpp.

{
 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;
}

Referenced by main().