hireconst_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 <cmath>

Include dependency graph for hireconst_test.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)
ErrorCode	create_unitsq_tris (Interface *mbImpl, size_t n, std::vector< EntityHandle > &tris)
ErrorCode	create_unitsq_quads (Interface *mbImpl, size_t n, std::vector< EntityHandle > &quads)
ErrorCode	test_unitsq_tris ()
ErrorCode	test_unitsq_quads ()
ErrorCode	test_unitsphere ()
ErrorCode	test_unitcircle ()
ErrorCode	exact_error_torus (const double R, const double r, const double center[3], int npnts, double *pnt, double &error_l1, double &error_l2, double &error_li)
ErrorCode	project_exact_torus (Interface *mbImpl, EntityHandle meshset, int dim, const double R, const double r, const double center[3])
int	main (int argc, char *argv[])
ErrorCode	project_exact_torus (Interface *mbImpl, EntityHandle meshset, int dim, const double R, const double r, const double center[])
ErrorCode	exact_error_torus (const double R, const double r, const double center[], int npnts, double *pnts, double &error_l1, double &error_l2, double &error_li)

Function Documentation

compute_linear_coords()

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

Definition at line 418 of file hireconst_test.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(), test_unitcircle(), test_unitsphere(), test_unitsq_quads(), and test_unitsq_tris().

create_unitsq_quads()

ErrorCode create_unitsq_quads	(	Interface *	mbImpl,
		size_t	n,
		std::vector< EntityHandle > &	quads
	)

Definition at line 318 of file hireconst_test.cpp.

{
 if( n < 2 )
 {
 MB_SET_ERR( MB_FAILURE, "n must be at least 2" );
 }
 
 ErrorCode error;
 std::vector< EntityHandle > verts( n * n );
 size_t istr = quads.size();
 quads.resize( istr + ( n - 1 ) * ( n - 1 ) );
 double istep = 1.0 / (double)( n - 1 );
 
 for( size_t j = 0; j < n; ++j )
 {
 for( size_t i = 0; i < n; ++i )
 {
 double coord[3] = { i * istep, j * istep, 0 };
 error = mbImpl->create_vertex( coord, verts[j * n + i] );MB_CHK_ERR( error );
 }
 }
 
 for( size_t jj = 0; jj < n - 1; ++jj )
 {
 for( size_t ii = 0; ii < n - 1; ++ii )
 {
 EntityHandle conn[4] = { verts[jj * n + ii], verts[jj * n + ii + 1], verts[( jj + 1 ) * n + ii + 1],
 verts[( jj + 1 ) * n + ii] };
 error = mbImpl->create_element( MBQUAD, conn, 4, quads[istr + jj * ( n - 1 ) + ii] );MB_CHK_ERR( error );
 }
 }
 
 return error;
}

References moab::Interface::create_element(), moab::Interface::create_vertex(), moab::error(), ErrorCode, MB_CHK_ERR, MB_SET_ERR, and MBQUAD.

Referenced by test_unitsq_quads().

create_unitsq_tris()

ErrorCode create_unitsq_tris	(	Interface *	mbImpl,
		size_t	n,
		std::vector< EntityHandle > &	tris
	)

Definition at line 278 of file hireconst_test.cpp.

{
 if( n < 2 )
 {
 MB_SET_ERR( MB_FAILURE, "n must be at least 2" );
 }
 
 ErrorCode error;
 std::vector< EntityHandle > verts( n * n );
 size_t istr = tris.size();
 tris.resize( istr + 2 * ( n - 1 ) * ( n - 1 ) );
 double istep = 1.0 / (double)( n - 1 );
 
 for( size_t j = 0; j < n; ++j )
 {
 for( size_t i = 0; i < n; ++i )
 {
 double coord[3] = { i * istep, j * istep, 0 };
 EntityHandle temp;
 error = mbImpl->create_vertex( coord, temp );MB_CHK_ERR( error );
 verts[j * n + i] = temp;
 }
 }
 
 for( size_t jj = 0; jj < n - 1; ++jj )
 {
 for( size_t ii = 0; ii < n - 1; ++ii )
 {
 EntityHandle conn[3] = { verts[jj * n + ii], verts[( jj + 1 ) * n + ii + 1], verts[( jj + 1 ) * n + ii] };
 error = mbImpl->create_element( MBTRI, conn, 3, tris[istr + 2 * jj * ( n - 1 ) + 2 * ii] );MB_CHK_ERR( error );
 conn[0] = verts[jj * n + ii];
 conn[1] = verts[jj * n + ii + 1];
 conn[2] = verts[( jj + 1 ) * n + ii + 1];
 error = mbImpl->create_element( MBTRI, conn, 3, tris[istr + 2 * jj * ( n - 1 ) + 2 * ii + 1] );MB_CHK_ERR( error );
 }
 }
 
 return error;
}

References moab::Interface::create_element(), moab::Interface::create_vertex(), moab::error(), ErrorCode, MB_CHK_ERR, MB_SET_ERR, and MBTRI.

Referenced by test_unitsq_tris().

exact_error_torus() [1/2]

ErrorCode exact_error_torus	(	const double	R,
		const double	r,
		const double	center[3],
		int	npnts,
		double *	pnt,
		double &	error_l1,
		double &	error_l2,
		double &	error_li
	)

Referenced by test_mesh().

exact_error_torus() [2/2]

ErrorCode exact_error_torus	(	const double	R,
		const double	r,
		const double	center[],
		int	npnts,
		double *	pnts,
		double &	error_l1,
		double &	error_l2,
		double &	error_li
	)

Definition at line 691 of file hireconst_test.cpp.

{
 error_l1 = 0;
 error_l2 = 0;
 error_li = 0;
 double x = 0, y = 0, z = 0;
 double error_single = 0;
 
 for( int i = 0; i < npnts; i++ )
 {
 x = pnts[3 * i] - center[0];
 y = pnts[3 * i + 1] - center[1];
 z = pnts[3 * i + 2] - center[2];
 error_single = fabs( sqrt( pow( sqrt( pow( x, 2 ) + pow( y, 2 ) ) - R, 2 ) + pow( z, 2 ) ) - r );
 error_l1 = error_l1 + error_single;
 error_l2 = error_l2 + error_single * error_single;
 
 if( error_li < error_single )
 {
 error_li = error_single;
 }
 }
 
 error_l1 = error_l1 / (double)npnts;
 error_l2 = sqrt( error_l2 / (double)npnts );
 return MB_SUCCESS;
}

References center(), MB_SUCCESS, and moab::R.

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 153 of file hireconst_test.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;
 
 if( nghlayers )
 {
 // get ghost layers
 if( dim == 2 )
 {
 read_options = "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARALLEL_RESOLVE_"
 "SHARED_ENTS;PARALLEL_GHOST=2.0.";
 }
 else if( dim == 1 )
 {
 read_options = "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARALLEL_RESOLVE_"
 "SHARED_ENTS;PARALLEL_GHOST=1.0.";
 }
 else
 {
 MB_SET_ERR( MB_FAILURE, "unsupported dimension" );
 }
 
 read_options += ( '0' + nghlayers );
 }
 else
 {
 read_options = "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;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(), test_unitcircle(), and test_unitsphere().

main()

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

Definition at line 70 of file hireconst_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;
 int degree = 2, dim = 0;
 bool interp = false;
 ErrorCode error;
 
 if( argc == 1 )
 {
 error = test_unitsq_tris();MB_CHK_ERR( error );
 error = test_unitsq_quads();MB_CHK_ERR( error );
 error = test_unitsphere();MB_CHK_ERR( error );
 error = test_unitcircle();MB_CHK_ERR( error );
#ifdef MOAB_HAVE_MPI
 MPI_Finalize();
#endif
 return 0;
 }
 else
 {
 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
 {
 std::cout << argv[i] << std::endl;
 std::cout << "usage: " << argv[0]
 << " <mesh file> -degree <degree> -interp <0=least square, "
 "1=interpolation> -dim <mesh dimension>"
 << std::endl;
 return 0;
 }
 }
 }
 
 if( !hasdim )
 {
 std::cout << "Dimension of input mesh should be provided, positive and less than 3" << std::endl;
 return 0;
 }
 
 if( degree <= 0 || dim > 2 || dim <= 0 )
 {
 std::cout << "Input degree should be positive number;\n";
 std::cout << "Input dimesion should be positive and less than 3;" << std::endl;
 return -1;
 }
 
 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;
 }
 
 error = test_mesh( infile.c_str(), degree, interp, dim );MB_CHK_ERR( error );
#ifdef MOAB_HAVE_MPI
 MPI_Finalize();
#endif
 return 0;
}

References dim, moab::error(), ErrorCode, MB_CHK_ERR, MPI_COMM_WORLD, rank, test_mesh(), test_unitcircle(), test_unitsphere(), test_unitsq_quads(), and test_unitsq_tris().

project_exact_torus() [1/2]

ErrorCode project_exact_torus	(	Interface *	mbImpl,
		EntityHandle	meshset,
		int	dim,
		const double	R,
		const double	r,
		const double	center[3]
	)

Referenced by test_mesh().

project_exact_torus() [2/2]

ErrorCode project_exact_torus	(	Interface *	mbImpl,
		EntityHandle	meshset,
		int	dim,
		const double	R,
		const double	r,
		const double	center[]
	)

Definition at line 653 of file hireconst_test.cpp.

{
 ErrorCode error;
 Range elems, verts;
 error = mbImpl->get_entities_by_dimension( meshset, dim, elems );CHECK_ERR( error );
 error = mbImpl->get_connectivity( elems, verts );CHECK_ERR( error );
 double pnts[3] = { 0, 0, 0 }, cnt[3] = { 0, 0, 0 }, nrms[3] = { 0, 0, 0 };
 double x = 0, y = 0, z = 0, d1 = 0, d2 = 0;
 
 for( int i = 0; i < (int)verts.size(); i++ )
 {
 EntityHandle v = verts[i];
 error = mbImpl->get_coords( &v, 1, &pnts[0] );CHECK_ERR( error );
 x = pnts[0] - center[0];
 y = pnts[1] - center[0];
 z = pnts[2] - center[2];
 d1 = sqrt( x * x + y * y );
 cnt[0] = R * x / d1;
 cnt[1] = R * y / d1;
 cnt[2] = 0;
 d2 = sqrt( ( x - cnt[0] ) * ( x - cnt[0] ) + ( y - cnt[1] ) * ( y - cnt[1] ) + z * z );
 nrms[0] = ( x - cnt[0] ) / d2;
 nrms[1] = ( y - cnt[1] ) / d2;
 nrms[2] = z / d2;
 pnts[0] = cnt[0] + r * nrms[0];
 pnts[1] = cnt[1] + r * nrms[1];
 pnts[2] = cnt[2] + r * nrms[2];
 error = mbImpl->set_coords( &v, 1, &pnts[0] );CHECK_ERR( error );
 }
 
 return MB_SUCCESS;
}

References center(), CHECK_ERR, dim, moab::error(), ErrorCode, moab::Interface::get_connectivity(), moab::Interface::get_coords(), moab::Interface::get_entities_by_dimension(), MB_SUCCESS, moab::R, moab::Interface::set_coords(), and moab::Range::size().

test_mesh()

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

Definition at line 218 of file hireconst_test.cpp.

{
 Core moab;
 Interface* mbimpl = &moab;
 ParallelComm* pc = NULL;
 EntityHandle meshset;
 // load mesh file
 ErrorCode error;
 error = load_meshset_hirec( infile, mbimpl, meshset, pc, degree, dim );MB_CHK_ERR( error );
 // project to exact surface: torus
 double center[3] = { 0, 0, 0 };
 double R = 1, r = 0.3;
 error = project_exact_torus( mbimpl, meshset, dim, R, r, center );CHECK_ERR( error );
 // initialize
 HiReconstruction hirec( dynamic_cast< Core* >( mbimpl ), pc, meshset );
 Range elems;
 error = mbimpl->get_entities_by_dimension( meshset, dim, elems );MB_CHK_ERR( error );
 
 // reconstruction
 if( dim == 2 )
 {
 // error = hirec.reconstruct3D_surf_geom(degree, interp, false); MB_CHK_ERR(error);
 error = hirec.reconstruct3D_surf_geom( degree, interp, true );MB_CHK_ERR( error );
 }
 else if( dim == 1 )
 {
 error = hirec.reconstruct3D_curve_geom( degree, interp, true );MB_CHK_ERR( error );
 }
 
 // fitting
 double mxdist = 0, errl1 = 0, errl2 = 0, errli = 0;
 double* pnts_proj = new double[elems.size() * 3];
 
 for( Range::iterator ielem = elems.begin(); ielem != elems.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 );
 double newcoords[3], linearcoords[3];
 error = hirec.hiproj_walf_in_element( *ielem, nvpe, 1, &( naturalcoords2fit[0] ), newcoords );MB_CHK_ERR( error );
 pnts_proj[3 * ( *ielem - *elems.begin() )] = newcoords[0];
 pnts_proj[3 * ( *ielem - *elems.begin() ) + 1] = newcoords[1];
 pnts_proj[3 * ( *ielem - *elems.begin() ) + 2] = newcoords[2];
 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 );
 mxdist = std::max( mxdist, DGMSolver::vec_distance( 3, newcoords, linearcoords ) );
 }
 
 delete[] pnts_proj;
 // compute error for torus
 error = exact_error_torus( R, r, center, (int)elems.size(), pnts_proj, errl1, errl2, errli );MB_CHK_ERR( error );
 std::cout << "Errors using exact torus for degree " << degree << " fit : L1 = " << errl1 << ", L2 = " << errl2
 << ", Linf = " << errli << std::endl;
 std::cout << "Maximum projection lift is " << mxdist << std::endl;
 return error;
}

References moab::Range::begin(), center(), CHECK_ERR, compute_linear_coords(), dim, moab::Range::end(), moab::error(), ErrorCode, exact_error_torus(), moab::Interface::get_connectivity(), moab::Interface::get_coords(), moab::Interface::get_entities_by_dimension(), moab::HiReconstruction::hiproj_walf_in_element(), load_meshset_hirec(), MB_CHK_ERR, project_exact_torus(), moab::R, moab::HiReconstruction::reconstruct3D_curve_geom(), moab::HiReconstruction::reconstruct3D_surf_geom(), moab::Range::size(), and moab::DGMSolver::vec_distance().

Referenced by main().

test_unitcircle()

ErrorCode test_unitcircle

(

)

Definition at line 574 of file hireconst_test.cpp.

{
 // path to test files
 int nfiles = 4;
 std::string filenames[] = { TestDir + "unittest/circle_3.vtk", TestDir + "unittest/circle_4.vtk",
 TestDir + "unittest/circle_10.vtk", TestDir + "unittest/circle_20.vtk" };
 ErrorCode error;
 int maxdeg = 6;
 
 for( int ifile = 0; ifile < nfiles; ++ifile )
 {
 Core moab;
 Interface* mbimpl = &moab;
 ParallelComm* pc = 0;
 EntityHandle meshset;
 int dim = 1;
 // load file
 error = load_meshset_hirec( filenames[ifile].c_str(), mbimpl, meshset, pc, maxdeg, dim );MB_CHK_ERR( error );
 // initialize
 HiReconstruction hirec( &moab, pc, meshset );
 Range edges;
 error = mbimpl->get_entities_by_dimension( meshset, dim, edges );MB_CHK_ERR( error );
 
 // reconstruction
 for( int degree = 1; degree <= maxdeg; ++degree )
 {
 hirec.reconstruct3D_curve_geom( degree, true, false, true );
 // fitting
 double mxdist = 0, mxerr = 0;
 
 for( Range::iterator iedge = edges.begin(); iedge != edges.end(); ++iedge )
 {
 int nvpe;
 const EntityHandle* conn;
 error = mbimpl->get_connectivity( *iedge, conn, nvpe );MB_CHK_ERR( error );
 double w = 1.0 / (double)nvpe;
 std::vector< double > naturalcoords2fit( nvpe, w );
 double newcoords[3], linearcoords[3];
 error = hirec.hiproj_walf_in_element( *iedge, nvpe, 1, &( naturalcoords2fit[0] ), newcoords );MB_CHK_ERR( error );
 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 );
 mxdist = std::max( mxdist, DGMSolver::vec_distance( 3, newcoords, linearcoords ) );
 mxerr = std::max( mxerr, fabs( DGMSolver::vec_2norm( 3, newcoords ) - 1 ) );
 }
 
 std::cout << filenames[ifile] << ": unit circle"
 << " degree= " << degree << " interpolation:\n";
 std::cout << "maximum projection lift is " << mxdist << ", maximum error is " << mxerr << std::endl;
 mxdist = 0;
 mxerr = 0;
 hirec.reconstruct3D_curve_geom( degree, false, false, true );
 
 // fitting
 for( Range::iterator iedge = edges.begin(); iedge != edges.end(); ++iedge )
 {
 int nvpe;
 const EntityHandle* conn;
 error = mbimpl->get_connectivity( *iedge, conn, nvpe );MB_CHK_ERR( error );
 double w = 1.0 / (double)nvpe;
 std::vector< double > naturalcoords2fit( nvpe, w );
 double newcoords[3], linearcoords[3];
 error = hirec.hiproj_walf_in_element( *iedge, nvpe, 1, &( naturalcoords2fit[0] ), newcoords );MB_CHK_ERR( error );
 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 );
 mxdist = std::max( mxdist, DGMSolver::vec_distance( 3, newcoords, linearcoords ) );
 mxerr = std::max( mxerr, fabs( DGMSolver::vec_2norm( 3, newcoords ) - 1 ) );
 }
 
 std::cout << filenames[ifile] << ": unit circle"
 << " degree= " << degree << " least square:\n";
 std::cout << "maximum projection lift is " << mxdist << ", maximum error is " << mxerr << std::endl;
 }
 }
 
 return error;
}

References moab::Range::begin(), compute_linear_coords(), dim, moab::Range::end(), moab::error(), ErrorCode, moab::Interface::get_connectivity(), moab::Interface::get_coords(), moab::Interface::get_entities_by_dimension(), moab::HiReconstruction::hiproj_walf_in_element(), load_meshset_hirec(), MB_CHK_ERR, moab::HiReconstruction::reconstruct3D_curve_geom(), moab::DGMSolver::vec_2norm(), and moab::DGMSolver::vec_distance().

Referenced by main().

test_unitsphere()

ErrorCode test_unitsphere

(

)

Definition at line 496 of file hireconst_test.cpp.

{
 // path to test files
 int nfiles = 4;
 std::string filenames[] = { TestDir + "unittest/sphere_tris_5.vtk", TestDir + "unittest/sphere_tris_20.vtk",
 TestDir + "unittest/sphere_quads_5.vtk", TestDir + "unittest/sphere_quads_20.vtk" };
 ErrorCode error;
 int maxdeg = 6;
 
 for( int ifile = 0; ifile < nfiles; ++ifile )
 {
 Core moab;
 Interface* mbimpl = &moab;
 ParallelComm* pc = NULL;
 EntityHandle meshset;
 // load file
 error = load_meshset_hirec( filenames[ifile].c_str(), mbimpl, meshset, pc, maxdeg );MB_CHK_ERR( error );
 // initialize
 HiReconstruction hirec( &moab, pc, meshset );
 Range elems;
 error = mbimpl->get_entities_by_dimension( meshset, 2, elems );MB_CHK_ERR( error );
 
 // reconstruction
 for( int degree = 1; degree <= maxdeg; ++degree )
 {
 hirec.reconstruct3D_surf_geom( degree, true, false, true );
 // fitting
 double mxdist = 0, mxerr = 0;
 
 for( Range::iterator ielem = elems.begin(); ielem != elems.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 );
 double newcoords[3], linearcoords[3];
 error = hirec.hiproj_walf_in_element( *ielem, nvpe, 1, &( naturalcoords2fit[0] ), newcoords );MB_CHK_ERR( error );
 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 );
 mxdist = std::max( mxdist, DGMSolver::vec_distance( 3, newcoords, linearcoords ) );
 mxerr = std::max( mxerr, fabs( DGMSolver::vec_2norm( 3, newcoords ) - 1 ) );
 }
 
 std::cout << filenames[ifile] << ": unit sphere"
 << " degree= " << degree << " interpolation:\n";
 std::cout << "maximum projection lift is " << mxdist << ", maximum error is " << mxerr << std::endl;
 mxdist = 0;
 mxerr = 0;
 hirec.reconstruct3D_surf_geom( degree, false, false, true );
 
 // fitting
 for( Range::iterator ielem = elems.begin(); ielem != elems.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 );
 double newcoords[3], linearcoords[3];
 error = hirec.hiproj_walf_in_element( *ielem, nvpe, 1, &( naturalcoords2fit[0] ), newcoords );MB_CHK_ERR( error );
 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 );
 mxdist = std::max( mxdist, DGMSolver::vec_distance( 3, newcoords, linearcoords ) );
 mxerr = std::max( mxerr, fabs( DGMSolver::vec_2norm( 3, newcoords ) - 1 ) );
 }
 
 std::cout << filenames[ifile] << ": unit sphere"
 << " degree= " << degree << " least square:\n";
 std::cout << "maximum projection lift is " << mxdist << ", maximum error is " << mxerr << std::endl;
 }
 }
 
 return error;
}

References moab::Range::begin(), compute_linear_coords(), moab::Range::end(), moab::error(), ErrorCode, moab::Interface::get_connectivity(), moab::Interface::get_coords(), moab::Interface::get_entities_by_dimension(), moab::HiReconstruction::hiproj_walf_in_element(), load_meshset_hirec(), MB_CHK_ERR, moab::HiReconstruction::reconstruct3D_surf_geom(), moab::DGMSolver::vec_2norm(), and moab::DGMSolver::vec_distance().

Referenced by main().

test_unitsq_quads()

ErrorCode test_unitsq_quads

(

)

Definition at line 433 of file hireconst_test.cpp.

{
 ErrorCode error;
 
 for( size_t n = 2; n <= 8; ++n )
 {
 Core moab;
 Interface* mbImpl = &moab;
 std::vector< EntityHandle > quads;
 error = create_unitsq_quads( mbImpl, n, quads );MB_CHK_ERR( error );
 EntityHandle meshIn = 0;
 HiReconstruction hirec( dynamic_cast< Core* >( mbImpl ), 0, meshIn );
 
 for( int degree = 1; degree <= 6; ++degree )
 {
 // reconstruct geometry, interpolation
 hirec.reconstruct3D_surf_geom( degree, true, false, true );
 // test fitting result
 double mxdist = 0;
 
 for( size_t iquad = 0; iquad < quads.size(); ++iquad )
 {
 const int nvpe = 4;
 double w = 1.0 / (double)nvpe;
 double naturalcoords2fit[nvpe] = { w, w, w, w }, newcoords[3];
 error = hirec.hiproj_walf_in_element( quads[iquad], nvpe, 1, naturalcoords2fit, newcoords );MB_CHK_ERR( error );
 std::vector< EntityHandle > conn;
 error = mbImpl->get_connectivity( &( quads[iquad] ), 1, conn );MB_CHK_ERR( error );
 double coords[3 * nvpe], linearcoords[3];
 error = mbImpl->get_coords( &( conn[0] ), nvpe, coords );MB_CHK_ERR( error );
 compute_linear_coords( nvpe, coords, naturalcoords2fit, linearcoords );
 mxdist = std::max( mxdist, DGMSolver::vec_distance( 3, newcoords, linearcoords ) );
 }
 
 std::cout << "quadrilateral unit square n= " << n << " degree= " << degree << " interpolation:\n";
 std::cout << "maximum projection lift is " << mxdist << std::endl;
 mxdist = 0;
 // reconstruct geometry, least square fitting
 hirec.reconstruct3D_surf_geom( degree, false, false, true );
 
 // test fitting result
 for( size_t iquad = 0; iquad < quads.size(); ++iquad )
 {
 const int nvpe = 4;
 double w = 1.0 / (double)nvpe;
 double naturalcoords2fit[nvpe] = { w, w, w, w }, newcoords[3];
 error = hirec.hiproj_walf_in_element( quads[iquad], nvpe, 1, naturalcoords2fit, newcoords );MB_CHK_ERR( error );
 std::vector< EntityHandle > conn;
 error = mbImpl->get_connectivity( &( quads[iquad] ), 1, conn );MB_CHK_ERR( error );
 double coords[3 * nvpe], linearcoords[3];
 error = mbImpl->get_coords( &( conn[0] ), nvpe, coords );MB_CHK_ERR( error );
 compute_linear_coords( nvpe, coords, naturalcoords2fit, linearcoords );
 mxdist = std::max( mxdist, DGMSolver::vec_distance( 3, newcoords, linearcoords ) );
 }
 
 std::cout << "quadrilateral unit square n= " << n << " degree= " << degree << " least square:\n";
 std::cout << "maximum projection lift is " << mxdist << std::endl;
 }
 }
 
 return error;
}

References compute_linear_coords(), create_unitsq_quads(), moab::error(), ErrorCode, moab::Interface::get_connectivity(), moab::Interface::get_coords(), moab::HiReconstruction::hiproj_walf_in_element(), MB_CHK_ERR, moab::HiReconstruction::reconstruct3D_surf_geom(), and moab::DGMSolver::vec_distance().

Referenced by main().

test_unitsq_tris()

ErrorCode test_unitsq_tris

(

)

Definition at line 353 of file hireconst_test.cpp.

{
 ErrorCode error;
 
 for( size_t n = 2; n <= 2; ++n )
 {
 Core moab;
 Interface* mbImpl = &moab;
 std::vector< EntityHandle > tris;
 error = create_unitsq_tris( mbImpl, n, tris );MB_CHK_ERR( error );
 EntityHandle meshIn = 0;
 HiReconstruction hirec( dynamic_cast< Core* >( mbImpl ), 0, meshIn );
 
 for( int degree = 2; degree <= 2; ++degree )
 {
 // reconstruct geometry, interpolation
 hirec.reconstruct3D_surf_geom( degree, true, true, true );
 // test fitting result
 double mxdist = 0;
 
 for( size_t itri = 0; itri < tris.size(); ++itri )
 {
 const int nvpe = 3;
 double naturalcoords2fit[nvpe] = { 1.0 / (double)nvpe, 1.0 / (double)nvpe, 1.0 / (double)nvpe },
 newcoords[3];
 error = hirec.hiproj_walf_in_element( tris[itri], nvpe, 1, naturalcoords2fit, newcoords );MB_CHK_ERR( error );
 std::vector< EntityHandle > conn;
 error = mbImpl->get_connectivity( &( tris[itri] ), 1, conn );MB_CHK_ERR( error );
 double coords[3 * nvpe], linearcoords[3];
 error = mbImpl->get_coords( &( conn[0] ), nvpe, coords );MB_CHK_ERR( error );
 compute_linear_coords( nvpe, coords, naturalcoords2fit, linearcoords );
 mxdist = std::max( mxdist, DGMSolver::vec_distance( 3, newcoords, linearcoords ) );
 }
 
 std::cout << "triangulated unit square n= " << n << " degree= " << degree << " interpolation:\n";
 std::cout << "maximum projection lift is " << mxdist << std::endl;
 // for debug
 // return error;
 mxdist = 0;
 // reconstruct geometry, least square fitting
 hirec.reconstruct3D_surf_geom( degree, false, false, true );
 
 // test fitting result
 for( size_t itri = 0; itri < tris.size(); ++itri )
 {
 const int nvpe = 3;
 double naturalcoords2fit[nvpe] = { 1.0 / (double)nvpe, 1.0 / (double)nvpe, 1.0 / (double)nvpe },
 newcoords[3];
 error = hirec.hiproj_walf_in_element( tris[itri], nvpe, 1, naturalcoords2fit, newcoords );MB_CHK_ERR( error );
 std::vector< EntityHandle > conn;
 error = mbImpl->get_connectivity( &( tris[itri] ), 1, conn );MB_CHK_ERR( error );
 double coords[3 * nvpe], linearcoords[3];
 error = mbImpl->get_coords( &( conn[0] ), nvpe, coords );MB_CHK_ERR( error );
 compute_linear_coords( nvpe, coords, naturalcoords2fit, linearcoords );
 mxdist = std::max( mxdist, DGMSolver::vec_distance( 3, newcoords, linearcoords ) );
 }
 
 std::cout << "unit square n= " << n << " degree= " << degree << " least square:\n";
 std::cout << "maximum projection lift is " << mxdist << std::endl;
 }
 }
 
 return error;
}

References compute_linear_coords(), create_unitsq_tris(), moab::error(), ErrorCode, moab::Interface::get_connectivity(), moab::Interface::get_coords(), moab::HiReconstruction::hiproj_walf_in_element(), MB_CHK_ERR, moab::HiReconstruction::reconstruct3D_surf_geom(), and moab::DGMSolver::vec_distance().

Referenced by main().