|
MOAB: Mesh Oriented datABase
(version 5.4.1)
|
|
MOAB: Mesh Oriented datABase
(version 5.4.1)
|
#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) |
| void compute_linear_coords | ( | const int | nvpe, |
| double * | elemcoords, | ||
| double * | naturals, | ||
| double * | linearcoords | ||
| ) |
Definition at line 418 of file hireconst_test.cpp.
Referenced by test_mesh(), test_unitcircle(), test_unitsphere(), test_unitsq_quads(), and test_unitsq_tris().
{
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];
}
}
}
| ErrorCode create_unitsq_quads | ( | Interface * | mbImpl, |
| size_t | n, | ||
| std::vector< EntityHandle > & | quads | ||
| ) |
Definition at line 318 of file hireconst_test.cpp.
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().
{
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;
}
| ErrorCode create_unitsq_tris | ( | Interface * | mbImpl, |
| size_t | n, | ||
| std::vector< EntityHandle > & | tris | ||
| ) |
Definition at line 278 of file hireconst_test.cpp.
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().
{
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;
}
| 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().
| 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.
References MB_SUCCESS.
{
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;
}
| 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.
References moab::Interface::create_meshset(), 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 closedsurface_uref_hirec_convergence_study(), test_mesh(), test_unitcircle(), and test_unitsphere().
{
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;
}
| int main | ( | int | argc, |
| char * | argv[] | ||
| ) |
Definition at line 70 of file hireconst_test.cpp.
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().
{
#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;
}
| ErrorCode project_exact_torus | ( | Interface * | mbImpl, |
| EntityHandle | meshset, | ||
| int | dim, | ||
| const double | R, | ||
| const double | r, | ||
| const double | center[3] | ||
| ) |
Referenced by test_mesh().
| 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.
References CHECK_ERR, moab::error(), ErrorCode, moab::Interface::get_connectivity(), moab::Interface::get_coords(), moab::Interface::get_entities_by_dimension(), MB_SUCCESS, moab::Interface::set_coords(), and moab::Range::size().
{
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;
}
Definition at line 218 of file hireconst_test.cpp.
References moab::Range::begin(), center(), CHECK_ERR, compute_linear_coords(), 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_1D(), test_2D(), and test_3D().
{
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;
}
Definition at line 574 of file hireconst_test.cpp.
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().
{
// 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;
}
Definition at line 496 of file hireconst_test.cpp.
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().
{
// 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;
}
Definition at line 433 of file hireconst_test.cpp.
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().
{
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;
}
Definition at line 353 of file hireconst_test.cpp.
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().
{
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;
}
1.7.6.1