uref_hirec_test.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/MeshTopoUtil.hpp"
#include "moab/NestedRefine.hpp"
#include "moab/DiscreteGeometry/DGMSolver.hpp"
#include "moab/DiscreteGeometry/HiReconstruction.hpp"
#include "TestUtil.hpp"
#include "geomObject.cpp"
#include <cmath>
#include <cstdlib>

Include dependency graph for uref_hirec_test.cpp:

Go to the source code of this file.

Macros
#define	nsamples   10

Functions
ErrorCode	test_closedsurface_mesh (const char *filename, int *level_degrees, int num_levels, int degree, bool interp, int dim, geomObject *obj)
ErrorCode	closedsurface_uref_hirec_convergence_study (const char *filename, int *level_degrees, int num_levels, std::vector< int > &degs2fit, bool interp, geomObject *obj)
void	usage ()
int	main (int argc, char *argv[])

Macro Definition Documentation

nsamples

#define nsamples   10

Definition at line 36 of file uref_hirec_test.cpp.

Function Documentation

closedsurface_uref_hirec_convergence_study()

ErrorCode closedsurface_uref_hirec_convergence_study	(	const char *	filename,
		int *	level_degrees,
		int	num_levels,
		std::vector< int > &	degs2fit,
		bool	interp,
		geomObject *	obj
	)

Definition at line 141 of file uref_hirec_test.cpp.

{
 Core moab;
 Interface* mbImpl = &moab;
 ParallelComm* pc = NULL;
 EntityHandle fileset;
 ErrorCode error;
 error = mbImpl->create_meshset( moab::MESHSET_SET, fileset );MB_CHK_ERR( error );
 // load mesh from file
#ifdef MOAB_HAVE_MPI
 MPI_Comm comm = MPI_COMM_WORLD;
 EntityHandle partnset;
 error = mbImpl->create_meshset( moab::MESHSET_SET, partnset );MB_CHK_ERR( error );
 pc = moab::ParallelComm::get_pcomm( mbImpl, partnset, &comm );
 int procs = 1;
 MPI_Comm_size( comm, &procs );
 
 if( procs > 1 )
 {
 read_options = "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARALLEL_RESOLVE_SHARED_ENTS;";
 error = mbImpl->load_file( filename, &fileset, read_options.c_str() );MB_CHK_ERR( error );
 }
 else
 {
 error = mbImpl->load_file( filename, &fileset );MB_CHK_ERR( error );
 }
 
#else
 error = mbImpl->load_file( filename, &fileset );MB_CHK_ERR( error );
#endif
 // Generate hierarchy
 NestedRefine uref( &moab, pc, fileset );
 std::vector< EntityHandle > meshes;
 std::vector< int > meshsizes;
 error = uref.generate_mesh_hierarchy( num_levels, level_degrees, meshes );MB_CHK_ERR( error );
 std::vector< std::vector< double > > geoml1errs( 1 + degs2fit.size() ), geoml2errs( 1 + degs2fit.size() ),
 geomlinferrs( 1 + degs2fit.size() );
 
 // Perform high order reconstruction on each level of mesh (projected onto exact geometry) and
 // estimate geometric error with various degrees of fitting
 for( size_t i = 0; i < meshes.size(); ++i )
 {
 EntityHandle& mesh = meshes[i];
 // project onto exact geometry since each level with uref has only linear coordinates
 Range verts;
 error = mbImpl->get_entities_by_dimension( mesh, 0, verts );MB_CHK_ERR( error );
 
 for( Range::iterator ivert = verts.begin(); ivert != verts.end(); ++ivert )
 {
 EntityHandle currvert = *ivert;
 double currcoords[3], exactcoords[3];
 error = mbImpl->get_coords( &currvert, 1, currcoords );
 obj->project_points2geom( 3, currcoords, exactcoords, NULL );
 error = mbImpl->set_coords( &currvert, 1, exactcoords );
 }
 
 // generate random points on each elements, assument 3D coordinates
 Range elems;
 error = mbImpl->get_entities_by_dimension( mesh, 2, elems );MB_CHK_ERR( error );
 meshsizes.push_back( elems.size() );
 int nvpe = TYPE_FROM_HANDLE( *elems.begin() ) == MBTRI ? 3 : 4;
 std::vector< double > testpnts, testnaturalcoords;
 testpnts.reserve( 3 * elems.size() * nsamples );
 testnaturalcoords.reserve( nvpe * elems.size() * nsamples );
 
 for( Range::iterator ielem = elems.begin(); ielem != elems.end(); ++ielem )
 {
 EntityHandle currelem = *ielem;
 std::vector< EntityHandle > conn;
 error = mbImpl->get_connectivity( &currelem, 1, conn );MB_CHK_ERR( error );
 nvpe = conn.size();
 std::vector< double > elemcoords( 3 * conn.size() );
 error = mbImpl->get_coords( &( conn[0] ), conn.size(), &( elemcoords[0] ) );MB_CHK_ERR( error );
 EntityType type = TYPE_FROM_HANDLE( currelem );
 
 for( int s = 0; s < nsamples; ++s )
 {
 if( type == MBTRI )
 {
 double a = (double)rand() / RAND_MAX, b = (double)rand() / RAND_MAX, c = (double)rand() / RAND_MAX,
 sum;
 sum = a + b + c;
 
 if( sum < 1e-12 )
 {
 --s;
 continue;
 }
 else
 {
 a /= sum, b /= sum, c /= sum;
 }
 
 testpnts.push_back( a * elemcoords[0] + b * elemcoords[3] + c * elemcoords[6] );
 testpnts.push_back( a * elemcoords[1] + b * elemcoords[4] + c * elemcoords[7] );
 testpnts.push_back( a * elemcoords[2] + b * elemcoords[5] + c * elemcoords[8] );
 testnaturalcoords.push_back( a );
 testnaturalcoords.push_back( b );
 testnaturalcoords.push_back( c );
 }
 else if( type == MBQUAD )
 {
 double xi = (double)rand() / RAND_MAX, eta = (double)rand() / RAND_MAX, N[4];
 xi = 2 * xi - 1;
 eta = 2 * eta - 1;
 N[0] = ( 1 - xi ) * ( 1 - eta ) / 4, N[1] = ( 1 + xi ) * ( 1 - eta ) / 4,
 N[2] = ( 1 + xi ) * ( 1 + eta ) / 4, N[3] = ( 1 - xi ) * ( 1 + eta ) / 4;
 testpnts.push_back( N[0] * elemcoords[0] + N[1] * elemcoords[3] + N[2] * elemcoords[6] +
 N[3] * elemcoords[9] );
 testpnts.push_back( N[0] * elemcoords[1] + N[1] * elemcoords[4] + N[2] * elemcoords[7] +
 N[3] * elemcoords[10] );
 testpnts.push_back( N[0] * elemcoords[2] + N[1] * elemcoords[5] + N[2] * elemcoords[8] +
 N[3] * elemcoords[11] );
 testnaturalcoords.push_back( N[0] );
 testnaturalcoords.push_back( N[1] );
 testnaturalcoords.push_back( N[2] );
 testnaturalcoords.push_back( N[3] );
 }
 else
 {
 throw std::invalid_argument( "NOT SUPPORTED ELEMENT TYPE" );
 }
 }
 }
 
 // Compute error of linear interpolation
 double l1err, l2err, linferr;
 obj->compute_projecterror( 3, elems.size() * nsamples, &( testpnts[0] ), l1err, l2err, linferr );
 geoml1errs[0].push_back( l1err );
 geoml2errs[0].push_back( l2err );
 geomlinferrs[0].push_back( linferr );
 // Perform high order projection and compute error
 HiReconstruction hirec( &moab, pc, mesh );
 
 for( size_t ideg = 0; ideg < degs2fit.size(); ++ideg )
 {
 // High order reconstruction
 error = hirec.reconstruct3D_surf_geom( degs2fit[ideg], interp, false, true );MB_CHK_ERR( error );
 int index = 0;
 
 for( Range::iterator ielem = elems.begin(); ielem != elems.end(); ++ielem, ++index )
 {
 // Projection
 error = hirec.hiproj_walf_in_element( *ielem, nvpe, nsamples,
 &( testnaturalcoords[nvpe * nsamples * index] ),
 &( testpnts[3 * nsamples * index] ) );MB_CHK_ERR( error );
 }
 
 // Estimate error
 obj->compute_projecterror( 3, elems.size() * nsamples, &( testpnts[0] ), l1err, l2err, linferr );
 geoml1errs[ideg + 1].push_back( l1err );
 geoml2errs[ideg + 1].push_back( l2err );
 geomlinferrs[ideg + 1].push_back( linferr );
 }
 }
 
 std::cout << "Mesh Size: ";
 
 for( size_t i = 0; i < meshsizes.size(); ++i )
 {
 std::cout << meshsizes[i] << " ";
 }
 
 std::cout << std::endl;
 std::cout << "Degrees: 0 ";
 
 for( size_t ideg = 0; ideg < degs2fit.size(); ++ideg )
 {
 std::cout << degs2fit[ideg] << " ";
 }
 
 std::cout << std::endl;
 std::cout << "L1-norm error: \n";
 
 for( size_t i = 0; i < geoml1errs.size(); ++i )
 {
 for( size_t j = 0; j < geoml1errs[i].size(); ++j )
 {
 std::cout << geoml1errs[i][j] << " ";
 }
 
 std::cout << std::endl;
 }
 
 std::cout << "L2-norm error: \n";
 
 for( size_t i = 0; i < geoml2errs.size(); ++i )
 {
 for( size_t j = 0; j < geoml2errs[i].size(); ++j )
 {
 std::cout << geoml2errs[i][j] << " ";
 }
 
 std::cout << std::endl;
 }
 
 std::cout << "Linf-norm error: \n";
 
 for( size_t i = 0; i < geomlinferrs.size(); ++i )
 {
 for( size_t j = 0; j < geomlinferrs[i].size(); ++j )
 {
 std::cout << geomlinferrs[i][j] << " ";
 }
 
 std::cout << std::endl;
 }
 
 return error;
}

References moab::Range::begin(), geomObject::compute_projecterror(), moab::Interface::create_meshset(), moab::Range::end(), moab::error(), ErrorCode, filename, moab::NestedRefine::generate_mesh_hierarchy(), moab::Interface::get_connectivity(), moab::Interface::get_coords(), moab::Interface::get_entities_by_dimension(), moab::ParallelComm::get_pcomm(), moab::HiReconstruction::hiproj_walf_in_element(), moab::Interface::load_file(), MB_CHK_ERR, MBQUAD, MBTRI, mesh, MESHSET_SET, MPI_COMM_WORLD, nsamples, geomObject::project_points2geom(), read_options, moab::HiReconstruction::reconstruct3D_surf_geom(), moab::Interface::set_coords(), moab::Range::size(), moab::sum(), and moab::TYPE_FROM_HANDLE().

Referenced by main().

main()

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

Definition at line 367 of file uref_hirec_test.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
 
 std::string infile = TestDir + "unittest/sphere_tris_5.vtk";
 
 int degree = 2, dim = 2, geom = 0;
 bool interp = false;
 ErrorCode error;
 
 if( argc == 10 )
 {
 infile = std::string( 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 if( std::string( argv[i] ) == "-geom" )
 {
 geom = std::string( argv[++i] ) == "s" ? 0 : 1;
 }
 else
 {
#ifdef MOAB_HAVE_MPI
 
 if( 0 == rank )
 {
 usage();
 }
 MPI_Finalize();
 
#else
 usage();
#endif
 return -1;
 }
 }
 }
 
 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 -1;
 }
 
 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 -1;
 }
 
#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
 }
 else
 {
 if( argc > 1 )
 {
 usage();
 return -1;
 }
 }
 
 {
 int level_degrees[3] = { 2, 2, 2 };
 int num_levels = 3;
 
 // create the geometry object
 geomObject* obj;
 if( geom )
 obj = new torus();
 else
 obj = new sphere();
 
 error = test_closedsurface_mesh( infile.c_str(), level_degrees, num_levels, degree, interp, dim, obj );MB_CHK_ERR( error );
 
 std::vector< int > degs2fit( 6 );
 for( int d = 1; d <= 6; ++d )
 {
 degs2fit[d - 1] = d;
 }
 
 // Call the higher order reconstruction routines and compute error convergence
 error = closedsurface_uref_hirec_convergence_study( infile.c_str(), level_degrees, num_levels, degs2fit, interp,
 obj );MB_CHK_ERR( error );
 // cleanup memory
 delete obj;
 }
 
#ifdef MOAB_HAVE_MPI
 MPI_Finalize();
#endif
 return 0;
}

References closedsurface_uref_hirec_convergence_study(), dim, moab::error(), ErrorCode, geom, MB_CHK_ERR, MPI_COMM_WORLD, rank, test_closedsurface_mesh(), and usage().

test_closedsurface_mesh()

ErrorCode test_closedsurface_mesh	(	const char *	filename,
		int *	level_degrees,
		int	num_levels,
		int	degree,
		bool	interp,
		int	dim,
		geomObject *	obj
	)

Definition at line 42 of file uref_hirec_test.cpp.

{
 Core moab;
 Interface* mbImpl = &moab;
 ParallelComm* pc = NULL;
 EntityHandle fileset;
 ErrorCode error;
 error = mbImpl->create_meshset( moab::MESHSET_SET, fileset );MB_CHK_ERR( error );
 // load mesh from file
#ifdef MOAB_HAVE_MPI
 MPI_Comm comm = MPI_COMM_WORLD;
 EntityHandle partnset;
 error = mbImpl->create_meshset( moab::MESHSET_SET, partnset );MB_CHK_ERR( error );
 pc = moab::ParallelComm::get_pcomm( mbImpl, partnset, &comm );
 int procs = 1;
 MPI_Comm_size( comm, &procs );
 
 if( procs > 1 )
 {
 read_options = "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARALLEL_RESOLVE_SHARED_ENTS;";
 error = mbImpl->load_file( filename, &fileset, read_options.c_str() );MB_CHK_ERR( error );
 }
 else
 {
 error = mbImpl->load_file( filename, &fileset );MB_CHK_ERR( error );
 }
 
#else
 error = mbImpl->load_file( filename, &fileset );MB_CHK_ERR( error );
#endif
 
 // Generate hierarchy
 // error = refine_entities(&moab, pc, fileset, level_degrees, num_levels, true);
 // CHECK_ERR(error);
 NestedRefine uref( &moab, pc, fileset );
 std::vector< EntityHandle > meshes;
 error = uref.generate_mesh_hierarchy( num_levels, level_degrees, meshes );MB_CHK_ERR( error );
 
 // Perform high order reconstruction on level 0 mesh
 HiReconstruction hirec( &moab, pc, fileset );
 assert( dim == 2 );
 error = hirec.reconstruct3D_surf_geom( degree, interp, false );MB_CHK_ERR( error );
 
 // High order projection
 Range verts, vorig;
 error = mbImpl->get_entities_by_dimension( meshes.back(), 0, verts );MB_CHK_ERR( error );
 error = mbImpl->get_entities_by_dimension( meshes.front(), 0, vorig );MB_CHK_ERR( error );
 int nvorig = vorig.size();
 double l1err = 0, l2err = 0, linferr = 0;
 
 for( Range::iterator ivert = verts.begin() + nvorig; ivert != verts.end(); ++ivert )
 {
 // locate the element in level 0 mesh, on which *ivert is lying
 EntityHandle currvert = *ivert;
 
 std::vector< EntityHandle > parentEntities;
 error = uref.get_adjacencies( *ivert, 2, parentEntities );MB_CHK_ERR( error );
 assert( parentEntities.size() );
 
 EntityHandle rootelem;
 error = uref.child_to_parent( parentEntities[0], num_levels, 0, &rootelem );MB_CHK_ERR( error );
 
 // compute the natural coordinates of *ivert in this element
 assert( TYPE_FROM_HANDLE( rootelem ) == MBTRI );
 
 const EntityHandle* conn;
 int nvpe = 3;
 error = mbImpl->get_connectivity( rootelem, conn, nvpe );MB_CHK_ERR( error );
 
 std::vector< double > cornercoords( 3 * nvpe ), currcoords( 3 );
 error = mbImpl->get_coords( conn, nvpe, &( cornercoords[0] ) );MB_CHK_ERR( error );
 error = mbImpl->get_coords( &currvert, 1, &( currcoords[0] ) );MB_CHK_ERR( error );
 
 std::vector< double > naturalcoords2fit( 3 );
 DGMSolver::get_tri_natural_coords( 3, &( cornercoords[0] ), 1, &( currcoords[0] ), &( naturalcoords2fit[0] ) );
 
 // project onto the estimated geometry
 double hicoords[3];
 error = hirec.hiproj_walf_in_element( rootelem, nvpe, 1, &( naturalcoords2fit[0] ), hicoords );MB_CHK_ERR( error );
 
 // estimate error
 double err = obj->compute_projecterror( 3, hicoords );
 l1err += err;
 l2err += err * err;
 linferr = std::max( linferr, err );
 }
 
 l1err /= verts.size() - nvorig;
 l2err = sqrt( l2err / ( verts.size() - nvorig ) );
 std::cout << "L1 error " << l1err << " L2 error " << l2err << " Linf error " << linferr << std::endl;
 return error;
}

References moab::Range::begin(), moab::NestedRefine::child_to_parent(), geomObject::compute_projecterror(), moab::Interface::create_meshset(), dim, moab::Range::end(), moab::error(), ErrorCode, filename, moab::NestedRefine::generate_mesh_hierarchy(), moab::NestedRefine::get_adjacencies(), moab::Interface::get_connectivity(), moab::Interface::get_coords(), moab::Interface::get_entities_by_dimension(), moab::ParallelComm::get_pcomm(), moab::DGMSolver::get_tri_natural_coords(), moab::HiReconstruction::hiproj_walf_in_element(), moab::Interface::load_file(), MB_CHK_ERR, MBTRI, MESHSET_SET, MPI_COMM_WORLD, read_options, moab::HiReconstruction::reconstruct3D_surf_geom(), moab::Range::size(), and moab::TYPE_FROM_HANDLE().

Referenced by main().

usage()

void usage

(

)

Definition at line 357 of file uref_hirec_test.cpp.

{
 std::cout << "usage: ./uref_hirec_test <mesh file> -degree <degree> -interp <0=least square, "
 "1=interpolation> -dim <mesh dimension> -geom <s=sphere,t=torus>"
 << std::endl;
 std::cout << "Example: ./uref_hirec_test $MOAB_SOURCE_DIR/MeshFiles/unittest/sphere_tris_5.vtk "
 "-degree 2 -interp 0 -dim 2 -geom s"
 << std::endl;
}

Referenced by main().