moab::HiReconstruction Class Reference

#include <HiReconstruction.hpp>

Collaboration diagram for moab::HiReconstruction:

Public Member Functions
	HiReconstruction (Core *impl, ParallelComm *comm=0, EntityHandle meshIn=0, int minpnts=5, bool recwhole=true)
	~HiReconstruction ()
ErrorCode	initialize (bool recwhole)
ErrorCode	reconstruct3D_surf_geom (int degree, bool interp, bool safeguard, bool reset=false)
	Reconstruct a high order surface on given surface mesh. More...
ErrorCode	reconstruct3D_surf_geom (size_t npts, int *degrees, bool *interps, bool safeguard, bool reset=false)
	Reconstruct a high order surface on given surface mesh. More...
ErrorCode	reconstruct3D_curve_geom (int degree, bool interp, bool safeguard, bool reset=false)
	Reconstruct a high order curve on given curve mesh. More...
ErrorCode	reconstruct3D_curve_geom (size_t npts, int *degrees, bool *interps, bool safeguard, bool reset=false)
	Reconstruct a high order curve on given curve mesh. More...
ErrorCode	polyfit3d_walf_surf_vertex (const EntityHandle vid, const bool interp, int degree, int minpnts, const bool safeguard, const int ncoords, double *coords, int *degree_out, const int ncoeffs, double *coeffs)
	Construct vertex based polynomial fitting on a surface mesh. More...
ErrorCode	polyfit3d_walf_curve_vertex (const EntityHandle vid, const bool interp, int degree, int minpnts, const bool safeguard, const int ncoords, double *coords, int *degree_out, const int ncoeffs, double *coeffs)
	Construct vertex based polynomial fitting on a curve mesh. More...
ErrorCode	hiproj_walf_in_element (EntityHandle elem, const int nvpe, const int npts2fit, const double *naturalcoords2fit, double *newcoords)
	Perform high order projection of points in an element, using estimated geometry by reconstruction class. More...
ErrorCode	hiproj_walf_around_vertex (EntityHandle vid, const int npts2fit, const double *coords2fit, double *hiproj_new)
	Perform high order projection of points around a vertex, using estimated geometry by reconstruction class. More...
void	walf3d_surf_vertex_eval (const double *local_origin, const double *local_coords, const int local_deg, const double *local_coeffs, const bool interp, const int npts2fit, const double *coords2fit, double *hiproj_new)
	Perform high order projection of points around a center vertex, assume geometry is surface. More...
void	walf3d_curve_vertex_eval (const double *local_origin, const double *local_coords, const int local_deg, const double *local_coeffs, const bool interp, const int npts2fit, const double *coords2fit, double *hiproj_new)
	Perform high order projection of points around a center vertex, assume geometry is curve. More...
bool	get_fittings_data (EntityHandle vid, GEOMTYPE &geomtype, std::vector< double > &coords, int &degree_out, std::vector< double > &coeffs, bool &interp)
	Get interally stored fitting results. More...

Static Public Member Functions
static int	estimate_num_ghost_layers (int degree, bool interp=false)

Protected Member Functions
	HiReconstruction (const HiReconstruction &source)
HiReconstruction &	operator= (const HiReconstruction &right)
int	estimate_num_rings (int degree, bool interp)
ErrorCode	vertex_get_incident_elements (const EntityHandle &vid, const int elemdim, std::vector< EntityHandle > &adjents)
	Given a vertex, return the incident elements with dimension elemdim. More...
ErrorCode	obtain_nring_ngbvs (const EntityHandle vid, int ring, const int minpnts, Range &ngbvs)
	Get n-ring neighbor vertices, assuming curve/surface mesh, not volume mesh. More...
void	initialize_surf_geom (const int degree)
void	initialize_surf_geom (const size_t npts, const int *degrees)
void	initialize_3Dcurve_geom (const int degree)
void	initialize_3Dcurve_geom (const size_t npts, const int *degrees)
ErrorCode	average_vertex_normal (const EntityHandle vid, double *nrm)
ErrorCode	compute_average_vertex_normals_surf ()
ErrorCode	get_normals_surf (const Range &vertsh, double *nrms)
ErrorCode	average_vertex_tangent (const EntityHandle vid, double *tang)
ErrorCode	compute_average_vertex_tangents_curve ()
ErrorCode	get_tangents_curve (const Range &vertsh, double *tangs)
void	polyfit3d_surf_get_coeff (const int nverts, const double *ngbcors, const double *ngbnrms, int degree, const bool interp, const bool safeguard, const int ncoords, double *coords, const int ncoeffs, double *coeffs, int *degree_out, int *degree_pnt, int *degree_qr)
	Compute local coordinates system of a vertex, and perform vertex based polynomial fittings of local height function. More...
void	eval_vander_bivar_cmf (const int npts2fit, const double *us, const int ndim, double *bs, int degree, const double *ws, const bool interp, const bool safeguard, int *degree_out, int *degree_pnt, int *degree_qr)
	Form and solve Vandermonde system of bi-variables. More...
void	polyfit3d_curve_get_coeff (const int nverts, const double *ngbcors, const double *ngbtangs, int degree, const bool interp, const bool safeguard, const int ncoords, double *coords, const int ncoeffs, double *coeffs, int *degree_out)
	Compute local single variable coordinate system of a vertex, and perform vertex based polynomial fittings of three global coordinates axis. More...
void	eval_vander_univar_cmf (const int npts2fit, const double *us, const int ndim, double *bs, int degree, const double *ws, const bool interp, const bool safeguard, int *degree_out)
	Form and solve Vandermonde system of single-variables. More...
int	compute_weights (const int nrows, const int ncols, const double *us, const int nngbs, const double *ngbnrms, const int degree, const double toler, double *ws)
	Compute weights for points selected in weighted least square fittigns. More...
bool	check_barycentric_coords (const int nws, const double *naturalcoords)
	Check the correctness of barycentric coordination, wi>=0 and sum(wi)=1. More...

Protected Attributes
Core *	mbImpl
ParallelComm *	pcomm
HalfFacetRep *	ahf
EntityHandle	_mesh2rec
Range	_verts2rec
Range	_inverts
Range	_inedges
Range	_infaces
Range	_incells
size_t	_nv2rec
int	_MAXPNTS
int	_MINPNTS
double	_MINEPS
bool	_hasderiv
GEOMTYPE	_geom
int	_dim
bool	_hasfittings
bool	_initfittings
std::vector< double >	_local_coords
std::vector< double >	_local_fit_coeffs
std::vector< size_t >	_vertID2coeffID
std::vector< int >	_degrees_out
std::vector< bool >	_interps

Detailed Description

Definition at line 39 of file HiReconstruction.hpp.

Constructor & Destructor Documentation

HiReconstruction() [1/2]

moab::HiReconstruction::HiReconstruction	(	Core *	impl,
		ParallelComm *	comm = 0,
		EntityHandle	meshIn = 0,
		int	minpnts = 5,
		bool	recwhole = true
	)

Definition at line 18 of file HiReconstruction.cpp.

 : mbImpl( impl ), pcomm( comm ), _mesh2rec( meshIn ), _MINPNTS( minpnts )
{
 assert( NULL != impl );
 ErrorCode error;
 _MINEPS = 1e-12;
 _dim = 0;
 _hasfittings = false;
 _hasderiv = false;
#ifdef MOAB_HAVE_MPI
 if( !pcomm )
 {
 pcomm = moab::ParallelComm::get_pcomm( mbImpl, 0 );
 }
#endif
 
 error = initialize( recwhole );
 if( MB_SUCCESS != error )
 {
 std::cout << "Error initializing HiReconstruction\n" << std::endl;
 exit( 1 );
 }
}

References _dim, _hasderiv, _hasfittings, _MINEPS, moab::error(), ErrorCode, moab::ParallelComm::get_pcomm(), initialize(), MB_SUCCESS, mbImpl, and pcomm.

~HiReconstruction()

moab::HiReconstruction::~HiReconstruction

(

)

Definition at line 42 of file HiReconstruction.cpp.

{
#ifdef MOAB_HAVE_AHF
 ahf = NULL;
#else
 delete ahf;
#endif
}

References ahf.

HiReconstruction() [2/2]

moab::HiReconstruction::HiReconstruction ( const HiReconstruction &  source )

protected

Member Function Documentation

average_vertex_normal()

ErrorCode moab::HiReconstruction::average_vertex_normal ( const EntityHandle  vid, double *  nrm  )

protected

Save fitting results of a vertex into internal storage

Parameters

vid	EntityHandle, a vertex in _mesh2rec, in _verts2rec
coords	Pointer to double array, global coordinates of local uvw coordinate system axis directions
degree_out	Integer, order of polynomial fittings for vid
coeffs	Pointer to double array, coefficients of local polynomial fittings in monomial basis
interp	Boolean, true = interpolation Compute area weighted average vertex normals for given vertex, assuming surface mesh For arbitrary polygon mesh, use incident two edges of each incident polygon of this vertex to form a triangle, then use these incident "triangles" to compute area weighted average vertex normals
vid	EntityHandle, vertex in _mesh2rec, might be ghost vertex
nrm	Pointer to 3-doubles array, preallocated by user

Definition at line 841 of file HiReconstruction.cpp.

{
 ErrorCode error;
 std::vector< EntityHandle > adjfaces;
 error = vertex_get_incident_elements( vid, 2, adjfaces );MB_CHK_ERR( error );
 int npolys = adjfaces.size();
 if( !npolys )
 {
 MB_SET_ERR( MB_FAILURE, "Vertex has no incident 2D entities" );
 }
 else
 {
 double v1[3], v2[3], v3[3], a[3], b[3], c[3];
 nrm[0] = nrm[1] = nrm[2] = 0;
 for( int i = 0; i < npolys; ++i )
 {
 // get incident "triangles"
 std::vector< EntityHandle > elemconn;
 error = mbImpl->get_connectivity( &adjfaces[i], 1, elemconn );MB_CHK_ERR( error );
 EntityHandle pre, nxt;
 int nvpe = elemconn.size();
 for( int j = 0; j < nvpe; ++j )
 {
 if( vid == elemconn[j] )
 {
 pre = j == 0 ? elemconn[nvpe - 1] : elemconn[j - 1];
 nxt = elemconn[( j + 1 ) % nvpe];
 break;
 }
 }
 // compute area weighted normals
 error = mbImpl->get_coords( &pre, 1, a );MB_CHK_ERR( error );
 error = mbImpl->get_coords( &vid, 1, b );MB_CHK_ERR( error );
 error = mbImpl->get_coords( &nxt, 1, c );MB_CHK_ERR( error );
 DGMSolver::vec_linear_operation( 3, 1, c, -1, b, v1 );
 DGMSolver::vec_linear_operation( 3, 1, a, -1, b, v2 );
 DGMSolver::vec_crossprod( v1, v2, v3 );
 DGMSolver::vec_linear_operation( 3, 1, nrm, 1, v3, nrm );
 }
#ifndef NDEBUG
 assert( DGMSolver::vec_normalize( 3, nrm, nrm ) );
#endif
 }
 return error;
}

References moab::error(), ErrorCode, moab::Core::get_connectivity(), moab::Core::get_coords(), MB_CHK_ERR, MB_SET_ERR, mbImpl, moab::DGMSolver::vec_crossprod(), moab::DGMSolver::vec_linear_operation(), moab::DGMSolver::vec_normalize(), and vertex_get_incident_elements().

Referenced by compute_average_vertex_normals_surf(), and get_normals_surf().

average_vertex_tangent()

ErrorCode moab::HiReconstruction::average_vertex_tangent ( const EntityHandle  vid, double *  tang  )

protected

Compute area weighted average vertex tangent vector for given vertex, assuming curve mesh Use incident two edges of vertex as estimatation of tangent vectors, weighted by length

Parameters

vid	EntityHandle, vertex in _mesh2rec, might be ghost vertex
tang	Pointer to 3-doubles array, preallocated by user

Definition at line 943 of file HiReconstruction.cpp.

{
 ErrorCode error;
 std::vector< EntityHandle > adjedges;
 error = vertex_get_incident_elements( vid, 1, adjedges );MB_CHK_ERR( error );
 int nedges = adjedges.size();
 if( !nedges )
 {
 MB_SET_ERR( MB_FAILURE, "Vertex has no incident edges" );
 }
 else
 {
 assert( nedges <= 2 );
 tang[0] = tang[1] = tang[2] = 0;
 for( int i = 0; i < nedges; ++i )
 {
 std::vector< EntityHandle > edgeconn;
 error = mbImpl->get_connectivity( &adjedges[i], 1, edgeconn );
 double istr[3], iend[3], t[3];
 error = mbImpl->get_coords( &( edgeconn[0] ), 1, istr );
 error = mbImpl->get_coords( &( edgeconn[1] ), 1, iend );
 DGMSolver::vec_linear_operation( 3, 1, iend, -1, istr, t );
 DGMSolver::vec_linear_operation( 3, 1, tang, 1, t, tang );
 }
#ifndef NDEBUG
 assert( DGMSolver::vec_normalize( 3, tang, tang ) );
#endif
 }
 return error;
}

References moab::error(), ErrorCode, moab::Core::get_connectivity(), moab::Core::get_coords(), MB_CHK_ERR, MB_SET_ERR, mbImpl, t, moab::DGMSolver::vec_linear_operation(), moab::DGMSolver::vec_normalize(), and vertex_get_incident_elements().

Referenced by compute_average_vertex_tangents_curve(), and get_tangents_curve().

check_barycentric_coords()

bool moab::HiReconstruction::check_barycentric_coords ( const int  nws, const double *  naturalcoords  )

protected

Check the correctness of barycentric coordination, wi>=0 and sum(wi)=1.

Definition at line 1593 of file HiReconstruction.cpp.

{
 double sum = 0;
 for( int i = 0; i < nws; ++i )
 {
 if( naturalcoords[i] < -_MINEPS )
 {
 return false;
 }
 sum += naturalcoords[i];
 }
 if( fabs( 1 - sum ) > _MINEPS )
 {
 return false;
 }
 else
 {
 return true;
 }
}

References _MINEPS, and moab::sum().

Referenced by hiproj_walf_in_element().

compute_average_vertex_normals_surf()

ErrorCode moab::HiReconstruction::compute_average_vertex_normals_surf ( )

protected

Compute weighted average vertex normals for all vertices in _verts2rec, not including ghost vertices, results are stored interally in _local_coords

Definition at line 887 of file HiReconstruction.cpp.

{
 if( _hasderiv )
 {
 return MB_SUCCESS;
 }
 ErrorCode error;
 _local_coords.assign( 9 * _nv2rec, 0 );
 size_t index = 0;
 for( Range::iterator ivert = _verts2rec.begin(); ivert != _verts2rec.end(); ++ivert, ++index )
 {
 error = average_vertex_normal( *ivert, &( _local_coords[9 * index + 6] ) );MB_CHK_ERR( error );
 }
 return error;
}

References _hasderiv, _local_coords, _nv2rec, _verts2rec, average_vertex_normal(), moab::Range::begin(), moab::Range::end(), moab::error(), ErrorCode, MB_CHK_ERR, and MB_SUCCESS.

Referenced by initialize(), and initialize_surf_geom().

compute_average_vertex_tangents_curve()

ErrorCode moab::HiReconstruction::compute_average_vertex_tangents_curve ( )

protected

Compute weighted average vertex tangent vectors for all vertices in _verts2rec, not including ghost vertices, results are stored interally in _local_coords

Definition at line 974 of file HiReconstruction.cpp.

{
 if( _hasderiv )
 {
 return MB_SUCCESS;
 }
 ErrorCode error;
 _local_coords.assign( 3 * _nv2rec, 0 );
 size_t index = 0;
 for( Range::iterator ivert = _verts2rec.begin(); ivert != _verts2rec.end(); ++ivert, ++index )
 {
 error = average_vertex_tangent( *ivert, &( _local_coords[3 * index] ) );MB_CHK_ERR( error );
 }
 return error;
}

References _hasderiv, _local_coords, _nv2rec, _verts2rec, average_vertex_tangent(), moab::Range::begin(), moab::Range::end(), moab::error(), ErrorCode, MB_CHK_ERR, and MB_SUCCESS.

Referenced by initialize(), and initialize_3Dcurve_geom().

compute_weights()

int moab::HiReconstruction::compute_weights ( const int  nrows, const int  ncols, const double *  us, const int  nngbs, const double *  ngbnrms, const int  degree, const double  toler, double *  ws  )

protected

Compute weights for points selected in weighted least square fittigns.

Definition at line 1544 of file HiReconstruction.cpp.

{
 assert( nrows <= _MAXPNTS && ws );
 bool interp = false;
 if( nngbs != nrows )
 {
 assert( nngbs == nrows + 1 );
 interp = true;
 }
 double epsilon = 1e-2;
 
 // First, compute squared distance from each input piont to the center
 for( int i = 0; i < nrows; ++i )
 {
 ws[i] = DGMSolver::vec_innerprod( ncols, us + i * ncols, us + i * ncols );
 }
 
 // Second, compute a small correction termt o guard against zero
 double h = 0;
 for( int i = 0; i < nrows; ++i )
 {
 h += ws[i];
 }
 h /= (double)nrows;
 
 // Finally, compute the weights for each vertex
 int nzeros = 0;
 for( int i = 0; i < nrows; ++i )
 {
 double costheta = DGMSolver::vec_innerprod( 3, ngbnrms, ngbnrms + 3 * ( i + interp ) );
 if( costheta > toler )
 {
 ws[i] = costheta * pow( ws[i] / h + epsilon, -1 * (double)degree / 2.0 );
 }
 else
 {
 ws[i] = 0;
 ++nzeros;
 }
 }
 return nzeros;
}

References _MAXPNTS, and moab::DGMSolver::vec_innerprod().

Referenced by polyfit3d_curve_get_coeff(), and polyfit3d_surf_get_coeff().

estimate_num_ghost_layers()

static int moab::HiReconstruction::estimate_num_ghost_layers ( int  degree, bool  interp = false  )

inlinestatic

Definition at line 277 of file HiReconstruction.hpp.

 {
 return 1 + ( interp ? ( ( degree + 1 ) >> 1 ) + ( ( degree + 1 ) & 1 )
 : ( ( degree + 2 ) >> 1 ) + ( ( degree + 2 ) & 1 ) );
 };

estimate_num_rings()

int moab::HiReconstruction::estimate_num_rings ( int  degree, bool  interp  )

protected

Definition at line 653 of file HiReconstruction.cpp.

{
 return interp ? ( ( degree + 1 ) >> 1 ) + ( ( degree + 1 ) & 1 ) : ( ( degree + 2 ) >> 1 ) + ( ( degree + 2 ) & 1 );
}

Referenced by polyfit3d_walf_curve_vertex(), and polyfit3d_walf_surf_vertex().

eval_vander_bivar_cmf()

void moab::HiReconstruction::eval_vander_bivar_cmf ( const int  npts2fit, const double *  us, const int  ndim, double *  bs, int  degree, const double *  ws, const bool  interp, const bool  safeguard, int *  degree_out, int *  degree_pnt, int *  degree_qr  )

protected

Form and solve Vandermonde system of bi-variables.

Definition at line 1176 of file HiReconstruction.cpp.

{
 // step 1. adjust the degree according to number of points to fit
 int ncols = ( ( ( degree + 2 ) * ( degree + 1 ) ) >> 1 ) - interp;
 while( 1 < degree && npts2fit < ncols )
 {
 --degree;
 ncols = ( ( ( degree + 2 ) * ( degree + 1 ) ) >> 1 ) - interp;
 }
 *degree_pnt = degree;
 
 // std::cout << "degree_pnt: " << *degree_pnt << std::endl;
 
 // step 2. construct Vandermonde matrix, stored in columnwise
 std::vector< double > V; // V(npts2fit*(ncols+interp)); //double *V_init = new double[npts2fit*(ncols+interp)];
 DGMSolver::gen_vander_multivar( npts2fit, 2, us, degree, V );
 // remove the first column of 1s if interpolation
 if( interp )
 {
 V.erase( V.begin(), V.begin() + npts2fit );
 }
 /*double* V;
 if(interp){
 V = new double[npts2fit*ncols];
 std::memcpy(V,V_init+npts2fit,ncols*npts2fit*sizeof(double));
 delete [] V_init; V_init = 0;
 }else{
 V = V_init;
 }*/
 
 // step 3. Scale rows to assign different weights to different points
 for( int i = 0; i < npts2fit; ++i )
 {
 for( int j = 0; j < ncols; ++j )
 {
 V[j * npts2fit + i] *= ws[i];
 }
 for( int k = 0; k < ndim; ++k )
 {
 bs[k * npts2fit + i] *= ws[i];
 }
 }
 
 // step 4. scale columns to reduce condition number
 std::vector< double > ts( ncols ); // double *ts = new double[ncols];
 DGMSolver::rescale_matrix( npts2fit, ncols, &( V[0] ), &( ts[0] ) );
 
 // step 5. Perform Householder QR factorization
 std::vector< double > D( ncols ); // double *D = new double[ncols];
 int rank;
 DGMSolver::qr_polyfit_safeguarded( npts2fit, ncols, &( V[0] ), &( D[0] ), &rank );
 
 // step 6. adjust degree of fitting according to rank of Vandermonde matrix
 int ncols_sub = ncols;
 while( rank < ncols_sub )
 {
 --degree;
 if( degree == 0 )
 {
 // surface is flat, return 0
 *degree_out = *degree_qr = degree;
 for( int i = 0; i < npts2fit; ++i )
 {
 for( int k = 0; k < ndim; ++k )
 {
 bs[k * npts2fit + i] = 0;
 }
 }
 return;
 }
 else
 {
 ncols_sub = ( ( ( degree + 2 ) * ( degree + 1 ) ) >> 1 ) - interp;
 }
 }
 *degree_qr = degree;
 
 // std::cout << "degree_qr: " << *degree_qr << std::endl;
 
 /* DBG
 * std::cout<<"before Qtb"<<std::endl;
 std::cout<<std::endl;
 std::cout<<"bs = "<<std::endl;
 std::cout<<std::endl;
 for (int k=0; k< ndim; k++){
 for (int j=0; j<npts2fit; ++j){
 std::cout<<" "<<bs[npts2fit*k+j]<<std::endl;
 }
 }
 std::cout<<std::endl;*/
 
 // step 7. compute Q'b
 DGMSolver::compute_qtransposeB( npts2fit, ncols_sub, &( V[0] ), ndim, bs );
 
 /* DBG
 * std::cout<<"after Qtb"<<std::endl;
 std::cout<<"bs = "<<std::endl;
 std::cout<<std::endl;
 for (int k=0; k< ndim; k++){
 for (int j=0; j<npts2fit; ++j){
 std::cout<<" "<<bs[npts2fit*k+j]<<std::endl;
 }
 }
 std::cout<<std::endl;*/
 
 // step 8. perform backward substitution and scale the solution
 // assign diagonals of V
 for( int i = 0; i < ncols_sub; ++i )
 {
 V[i * npts2fit + i] = D[i];
 }
 
 // backsolve
 if( safeguard )
 {
 // for debug
 // std::cout << "ts size " << ts.size() << std::endl;
 DGMSolver::backsolve_polyfit_safeguarded( 2, degree, interp, npts2fit, ncols_sub, &( V[0] ), ndim, bs,
 &( ts[0] ), degree_out );
 }
 else
 {
 DGMSolver::backsolve( npts2fit, ncols_sub, &( V[0] ), 1, bs, &( ts[0] ) );
 *degree_out = degree;
 }
 /*if(V_init){
 delete [] V_init;
 }else{
 delete [] V;
 }*/
}

References moab::DGMSolver::backsolve(), moab::DGMSolver::backsolve_polyfit_safeguarded(), moab::DGMSolver::compute_qtransposeB(), moab::DGMSolver::gen_vander_multivar(), moab::DGMSolver::qr_polyfit_safeguarded(), and moab::DGMSolver::rescale_matrix().

Referenced by polyfit3d_surf_get_coeff().

eval_vander_univar_cmf()

void moab::HiReconstruction::eval_vander_univar_cmf ( const int  npts2fit, const double *  us, const int  ndim, double *  bs, int  degree, const double *  ws, const bool  interp, const bool  safeguard, int *  degree_out  )

protected

Form and solve Vandermonde system of single-variables.

Definition at line 1434 of file HiReconstruction.cpp.

{
 // step 1. determine degree of polynomials to fit according to number of points
 int ncols = degree + 1 - interp;
 while( npts2fit < ncols && degree >= 1 )
 {
 --degree;
 ncols = degree + 1 - interp;
 }
 if( !degree )
 {
 if( interp )
 {
 for( int icol = 0; icol < ndim; ++icol )
 {
 bs[icol * npts2fit] = 0;
 }
 }
 for( int irow = 1; irow < npts2fit; ++irow )
 {
 for( int icol = 0; icol < ndim; ++icol )
 {
 bs[icol * npts2fit + irow] = 0;
 }
 }
 *degree_out = 0;
 return;
 }
 // step 2. construct Vandermonde matrix
 std::vector< double > V; // V(npts2fit*(ncols+interp));
 DGMSolver::gen_vander_multivar( npts2fit, 1, us, degree, V );
 
 if( interp )
 {
 V.erase( V.begin(), V.begin() + npts2fit );
 }
 
 // step 3. scale rows with respect to weights
 for( int i = 0; i < npts2fit; ++i )
 {
 for( int j = 0; j < ncols; ++j )
 {
 V[j * npts2fit + i] *= ws[i];
 }
 for( int k = 0; k < ndim; ++k )
 {
 bs[k * npts2fit + i] *= ws[i];
 }
 }
 
 // step 4. scale columns to reduce condition number
 std::vector< double > ts( ncols );
 DGMSolver::rescale_matrix( npts2fit, ncols, &( V[0] ), &( ts[0] ) );
 
 // step 5. perform Householder QR factorization
 std::vector< double > D( ncols );
 int rank;
 DGMSolver::qr_polyfit_safeguarded( npts2fit, ncols, &( V[0] ), &( D[0] ), &rank );
 
 // step 6. adjust degree of fitting
 int ncols_sub = ncols;
 while( rank < ncols_sub )
 {
 --degree;
 if( !degree )
 {
 // matrix is singular, consider curve on flat plane passing center and orthogonal to
 // tangent line
 *degree_out = 0;
 for( int i = 0; i < npts2fit; ++i )
 {
 for( int k = 0; k < ndim; ++k )
 {
 bs[k * npts2fit + i] = 0;
 }
 }
 }
 ncols_sub = degree + 1 - interp;
 }
 
 // step 7. compute Q'*bs
 DGMSolver::compute_qtransposeB( npts2fit, ncols_sub, &( V[0] ), ndim, bs );
 
 // step 8. perform backward substitution and scale solutions
 // assign diagonals of V
 for( int i = 0; i < ncols_sub; ++i )
 {
 V[i * npts2fit + i] = D[i];
 }
 // backsolve
 if( safeguard )
 {
 DGMSolver::backsolve_polyfit_safeguarded( 1, degree, interp, npts2fit, ncols, &( V[0] ), ndim, bs, ws,
 degree_out );
 }
 else
 {
 DGMSolver::backsolve( npts2fit, ncols_sub, &( V[0] ), ndim, bs, &( ts[0] ) );
 *degree_out = degree;
 }
}

References moab::DGMSolver::backsolve(), moab::DGMSolver::backsolve_polyfit_safeguarded(), moab::DGMSolver::compute_qtransposeB(), moab::DGMSolver::gen_vander_multivar(), moab::DGMSolver::qr_polyfit_safeguarded(), and moab::DGMSolver::rescale_matrix().

Referenced by polyfit3d_curve_get_coeff().

get_fittings_data()

bool moab::HiReconstruction::get_fittings_data	(	EntityHandle	vid,
		GEOMTYPE &	geomtype,
		std::vector< double > &	coords,
		int &	degree_out,
		std::vector< double > &	coeffs,
		bool &	interp
	)

Get interally stored fitting results.

Get fittings results of a vertex, stored internally, results will be writtend to user provided memory

Parameters

vid	EntityHandle, a vertex in _verts2rec
geomtype	GEOMTYPE, one of HISURFACE,HI3DCURVE,HI2DCURVE
coords	vector, global coordinates of local uvw coordinate system axis directions will be appended to the end of coords
degree_out	Reference to Integer, order of polynomial fittings for vid
coeffs	vector, coefficients of local polynomial fittings in monomial basis will be appended to the end of coeffs
interp	Reference to Boolean, true = interpolation

Definition at line 607 of file HiReconstruction.cpp.

{
 if( !_hasfittings )
 {
 return false;
 }
 else
 {
 int index = _verts2rec.index( vid );
 if( -1 == index )
 {
 std::cout << "Input vertex is not locally hosted vertex in this mesh set" << std::endl;
 return false;
 }
 geomtype = _geom;
 if( HISURFACE == _geom )
 {
 coords.insert( coords.end(), _local_coords.begin() + 9 * index, _local_coords.begin() + 9 * index + 9 );
 degree_out = _degrees_out[index];
 interp = _interps[index];
 int ncoeffs = ( degree_out + 2 ) * ( degree_out + 1 ) >> 1;
 size_t istr = _vertID2coeffID[index];
 coeffs.insert( coeffs.end(), _local_fit_coeffs.begin() + istr, _local_fit_coeffs.begin() + istr + ncoeffs );
 }
 else if( HI3DCURVE == _geom )
 {
 coords.insert( coords.end(), _local_coords.begin() + 3 * index, _local_coords.begin() + 3 * index + 3 );
 degree_out = _degrees_out[index];
 interp = _interps[index];
 int ncoeffs = 3 * ( degree_out + 1 );
 size_t istr = _vertID2coeffID[index];
 coeffs.insert( coeffs.end(), _local_fit_coeffs.begin() + istr, _local_fit_coeffs.begin() + istr + ncoeffs );
 }
 return true;
 }
}

References _degrees_out, _geom, _hasfittings, _interps, _local_coords, _local_fit_coeffs, _vertID2coeffID, _verts2rec, moab::HI3DCURVE, moab::HISURFACE, and moab::Range::index().

get_normals_surf()

ErrorCode moab::HiReconstruction::get_normals_surf ( const Range &  vertsh, double *  nrms  )

protected

Return the normals of given vertices in a Range, writing to preallocated memory If normals have been computed and stored, just access them If not, compute on the fly

Parameters

vertsh	Range, EntityHandles of vertices
nrms	Pointer of array of doubles, size = 3*vertsh.size()

Definition at line 903 of file HiReconstruction.cpp.

{
 ErrorCode error = MB_SUCCESS;
 if( _hasderiv )
 {
 size_t id = 0;
 for( Range::iterator ivert = vertsh.begin(); ivert != vertsh.end(); ++ivert, ++id )
 {
 int index = _verts2rec.index( *ivert );
#ifdef MOAB_HAVE_MPI
 if( -1 == index )
 {
 // ghost vertex
 error = average_vertex_normal( *ivert, nrms + 3 * id );MB_CHK_ERR( error );
 }
 else
 {
 nrms[3 * id] = _local_coords[9 * index + 6];
 nrms[3 * id + 1] = _local_coords[9 * index + 7];
 nrms[3 * id + 2] = _local_coords[9 * index + 8];
 }
#else
 assert( -1 != index );
 nrms[3 * id] = _local_coords[9 * index + 6];
 nrms[3 * id + 1] = _local_coords[9 * index + 7];
 nrms[3 * id + 2] = _local_coords[9 * index + 8];
#endif
 }
 }
 else
 {
 size_t id = 0;
 for( Range::iterator ivert = vertsh.begin(); ivert != vertsh.end(); ++ivert, ++id )
 {
 error = average_vertex_normal( *ivert, nrms + 3 * id );MB_CHK_ERR( error );
 }
 }
 return error;
}

References _hasderiv, _local_coords, _verts2rec, average_vertex_normal(), moab::Range::begin(), moab::Range::end(), moab::error(), ErrorCode, moab::Range::index(), MB_CHK_ERR, and MB_SUCCESS.

Referenced by polyfit3d_walf_surf_vertex().

get_tangents_curve()

ErrorCode moab::HiReconstruction::get_tangents_curve ( const Range &  vertsh, double *  tangs  )

protected

Return the tangent vectors of given vertices in a Range, writing to preallocated memory If tangent vectors have been computed and stored, just access them If not, compute on the fly

Parameters

vertsh	Range, EntityHandles of vertices
tangs	Pointer of array of doubles, size = 3*vertsh.size()

Definition at line 990 of file HiReconstruction.cpp.

{
 ErrorCode error = MB_SUCCESS;
 if( _hasderiv )
 {
 size_t id = 0;
 for( Range::iterator ivert = vertsh.begin(); ivert != vertsh.end(); ++ivert, ++id )
 {
 int index = _verts2rec.index( *ivert );
#ifdef MOAB_HAVE_MPI
 if( -1 != index )
 {
 tangs[3 * id] = _local_coords[3 * index];
 tangs[3 * id + 1] = _local_coords[3 * index + 1];
 tangs[3 * id + 2] = _local_coords[3 * index + 2];
 }
 else
 {
 error = average_vertex_tangent( *ivert, tangs + 3 * id );MB_CHK_ERR( error );
 }
#else
 assert( -1 != index );
 tangs[3 * id] = _local_coords[3 * index];
 tangs[3 * id + 1] = _local_coords[3 * index + 1];
 tangs[3 * id + 2] = _local_coords[3 * index + 2];
#endif
 }
 }
 else
 {
 size_t id = 0;
 for( Range::iterator ivert = vertsh.begin(); ivert != vertsh.end(); ++ivert, ++id )
 {
 error = average_vertex_tangent( *ivert, tangs + 3 * id );MB_CHK_ERR( error );
 }
 }
 return error;
}

References _hasderiv, _local_coords, _verts2rec, average_vertex_tangent(), moab::Range::begin(), moab::Range::end(), moab::error(), ErrorCode, moab::Range::index(), MB_CHK_ERR, and MB_SUCCESS.

Referenced by polyfit3d_walf_curve_vertex().

hiproj_walf_around_vertex()

ErrorCode moab::HiReconstruction::hiproj_walf_around_vertex	(	EntityHandle	vid,
		const int	npts2fit,
		const double *	coords2fit,
		double *	hiproj_new
	)

Perform high order projection of points around a vertex, using estimated geometry by reconstruction class.

Given an vertex on the input mesh, and new points around this vertex, estimate their position in surface. This is done by first projecting input points onto the local uv-plane around this vertex and use the precomputed local fitting to estimate the ideal position of input points. The result will be returned to the user preallocated memory

Parameters

vid	EntityHandle, the vertex around which to perform high order projection.
npts2fit	Integer, number of points lying around vid to be fitted.
coords2fit	Pointer to array of doubles, size=3*npts2fit, current coordinates of points to be projected.
newcoords	Pointer to array of doubles, preallocated by user, size=3*npts2fit, estimated positions of input points.

Definition at line 468 of file HiReconstruction.cpp.

{
 if( !_hasfittings )
 {
 MB_SET_ERR( MB_FAILURE, "There is no existing fitting results" );
 }
 else if( -1 == _verts2rec.index( vid ) )
 {
 std::ostringstream convert;
 convert << vid;
 std::string VID = convert.str();
 MB_SET_ERR( MB_FAILURE, "There is no existing fitting results for vertex " + VID );
 }
 ErrorCode error;
 // get center of local coordinates system
 double local_origin[3];
 error = mbImpl->get_coords( &vid, 1, local_origin );MB_CHK_ERR( error );
 // get local fitting parameters
 int index = _verts2rec.index( vid );
 bool interp = _interps[index];
 int local_deg = _degrees_out[index];
 double *uvw_coords, *local_coeffs;
 if( _geom == HISURFACE )
 {
 uvw_coords = &( _local_coords[9 * index] );
 // int ncoeffs = (local_deg+2)*(local_deg+1)>>1;
 size_t istr = _vertID2coeffID[index];
 local_coeffs = &( _local_fit_coeffs[istr] );
 walf3d_surf_vertex_eval( local_origin, uvw_coords, local_deg, local_coeffs, interp, npts2fit, coords2fit,
 hiproj_new );
 }
 else if( _geom == HI3DCURVE )
 {
 uvw_coords = &( _local_coords[3 * index] );
 size_t istr = _vertID2coeffID[index];
 local_coeffs = &( _local_fit_coeffs[istr] );
 walf3d_curve_vertex_eval( local_origin, uvw_coords, local_deg, local_coeffs, interp, npts2fit, coords2fit,
 hiproj_new );
 }
 return error;
}

References _degrees_out, _geom, _hasfittings, _interps, _local_coords, _local_fit_coeffs, _vertID2coeffID, _verts2rec, moab::error(), ErrorCode, moab::Core::get_coords(), moab::HI3DCURVE, moab::HISURFACE, moab::Range::index(), MB_CHK_ERR, MB_SET_ERR, mbImpl, walf3d_curve_vertex_eval(), and walf3d_surf_vertex_eval().

Referenced by hiproj_walf_in_element().

hiproj_walf_in_element()

ErrorCode moab::HiReconstruction::hiproj_walf_in_element	(	EntityHandle	elem,
		const int	nvpe,
		const int	npts2fit,
		const double *	naturalcoords2fit,
		double *	newcoords
	)

Perform high order projection of points in an element, using estimated geometry by reconstruction class.

Given an element on the input mesh, and new points in this element, represented as natural coordinates in element, estimate their position in surface. This is done by weighted averaging of local fittings: for each vertex of this elment, a fitting has been computed and the new points could be projected by this fitting. The final result of projection is the weighted average of these projections, weights are chosen as the barycentric coordinates of the point in this element. The result will be returned to the user preallocated memory

Parameters

elem	EntityHandle, the element on which to perform high order projection.
nvpe	Integer, number of nodes of this element, triangle is 3, quad is four.
npts2fit	Integer, number of points lying in elem to be projected.
naturalcoords2fit	Pointer to array of doubles, size=nvpe*npts2fit, natural coordinates in elem of points to be projected.
newcoords	Pointer to array of doubles, preallocated by user, size=3*npts2fit, estimated positions of input points.

Definition at line 389 of file HiReconstruction.cpp.

{
 assert( newcoords );
 ErrorCode error;
 // get connectivity table
 std::vector< EntityHandle > elemconn;
 error = mbImpl->get_connectivity( &elem, 1, elemconn );MB_CHK_ERR( error );
 if( nvpe != (int)elemconn.size() )
 {
 MB_SET_ERR( MB_FAILURE, "element connectivity table size doesn't match input size" );
 }
 
 if( !_hasfittings )
 {
 MB_SET_ERR( MB_FAILURE, "There is no existing fitting results" );
 }
 else
 {
 std::ostringstream convert;
 convert << elem;
 std::string ID = convert.str();
 for( int i = 0; i < nvpe; ++i )
 {
 if( -1 == _verts2rec.index( elemconn[i] ) )
 {
 MB_SET_ERR( MB_FAILURE, "There is no existing fitting results for element " + ID );
 }
 }
 }
 // check correctness of input
 for( int i = 0; i < npts2fit; ++i )
 {
 if( !check_barycentric_coords( nvpe, naturalcoords2fit + i * nvpe ) )
 {
 MB_SET_ERR( MB_FAILURE, "Wrong barycentric coordinates" );
 }
 }
 
 double* elemcoords = new double[nvpe * 3];
 error = mbImpl->get_coords( &( elemconn[0] ), nvpe, elemcoords );MB_CHK_ERR( error );
 
 double* coords2fit = new double[3 * npts2fit]();
 for( int i = 0; i < npts2fit; ++i )
 {
 for( int j = 0; j < nvpe; ++j )
 {
 coords2fit[3 * i] += naturalcoords2fit[i * nvpe + j] * elemcoords[3 * j];
 coords2fit[3 * i + 1] += naturalcoords2fit[i * nvpe + j] * elemcoords[3 * j + 1];
 coords2fit[3 * i + 2] += naturalcoords2fit[i * nvpe + j] * elemcoords[3 * j + 2];
 }
 }
 
 double* hiproj_new = new double[3 * npts2fit];
 // initialize output
 for( int i = 0; i < npts2fit; ++i )
 {
 newcoords[3 * i] = newcoords[3 * i + 1] = newcoords[3 * i + 2] = 0;
 }
 // for each input vertex, call nvpe fittings and take average
 for( int j = 0; j < nvpe; ++j )
 {
 error = hiproj_walf_around_vertex( elemconn[j], npts2fit, coords2fit, hiproj_new );MB_CHK_ERR( error );
 for( int i = 0; i < npts2fit; ++i )
 {
 newcoords[3 * i] += naturalcoords2fit[i * nvpe + j] * hiproj_new[3 * i];
 newcoords[3 * i + 1] += naturalcoords2fit[i * nvpe + j] * hiproj_new[3 * i + 1];
 newcoords[3 * i + 2] += naturalcoords2fit[i * nvpe + j] * hiproj_new[3 * i + 2];
 }
 }
 delete[] elemcoords;
 delete[] coords2fit;
 delete[] hiproj_new;
 return error;
}

References _hasfittings, _verts2rec, check_barycentric_coords(), moab::error(), ErrorCode, moab::Core::get_connectivity(), moab::Core::get_coords(), hiproj_walf_around_vertex(), moab::Range::index(), MB_CHK_ERR, MB_SET_ERR, and mbImpl.

initialize()

ErrorCode moab::HiReconstruction::initialize

(

bool

recwhole

)

Definition at line 51 of file HiReconstruction.cpp.

{
 ErrorCode error;
 
#ifdef HIREC_USE_AHF
 std::cout << "HIREC_USE_AHF: Initializing" << std::endl;
 ahf = new HalfFacetRep( mbImpl, pcomm, _mesh2rec, false );
 if( !ahf )
 {
 return MB_MEMORY_ALLOCATION_FAILED;
 }
 error = ahf->initialize();MB_CHK_ERR( error );
#else
 ahf = NULL;
#endif
 
 // error = ahf->get_entity_ranges(_inverts,_inedges,_infaces,_incells); MB_CHK_ERR(error);
 error = mbImpl->get_entities_by_dimension( _mesh2rec, 0, _inverts );MB_CHK_ERR( error );
 error = mbImpl->get_entities_by_dimension( _mesh2rec, 1, _inedges );MB_CHK_ERR( error );
 error = mbImpl->get_entities_by_dimension( _mesh2rec, 2, _infaces );MB_CHK_ERR( error );
 error = mbImpl->get_entities_by_dimension( _mesh2rec, 3, _incells );MB_CHK_ERR( error );
 if( _inedges.size() && _infaces.empty() && _incells.empty() )
 {
 _dim = 1;
 _MAXPNTS = 13;
 }
 else if( _infaces.size() && _incells.empty() )
 {
 _dim = 2;
 _MAXPNTS = 128;
 }
 else
 {
 MB_SET_ERR( MB_FAILURE, "Encountered a non-manifold mesh or a mesh with volume elements" );
 }
 
 // get locally hosted vertices by filtering pstatus
#ifdef MOAB_HAVE_MPI
 if( pcomm )
 {
 error = pcomm->filter_pstatus( _inverts, PSTATUS_GHOST, PSTATUS_NOT, -1, &_verts2rec );MB_CHK_ERR( error );
 }
 else
 {
 _verts2rec = _inverts;
 }
#else
 _verts2rec = _inverts;
#endif
 _nv2rec = _verts2rec.size();
 
 if( recwhole )
 {
 // compute normals(surface) or tangent vector(curve) for all locally hosted vertices
 if( 2 == _dim )
 {
 compute_average_vertex_normals_surf();
 }
 else if( 1 == _dim )
 {
 compute_average_vertex_tangents_curve();
 }
 else
 {
 MB_SET_ERR( MB_FAILURE, "Unknow space dimension" );
 }
 _hasderiv = true;
 }
 return error;
}

References _dim, _hasderiv, _incells, _inedges, _infaces, _inverts, _MAXPNTS, _mesh2rec, _nv2rec, _verts2rec, ahf, compute_average_vertex_normals_surf(), compute_average_vertex_tangents_curve(), moab::Range::empty(), moab::error(), ErrorCode, moab::ParallelComm::filter_pstatus(), moab::Core::get_entities_by_dimension(), moab::HalfFacetRep::initialize(), MB_CHK_ERR, MB_MEMORY_ALLOCATION_FAILED, MB_SET_ERR, mbImpl, pcomm, PSTATUS_GHOST, PSTATUS_NOT, and moab::Range::size().

Referenced by HiReconstruction().

initialize_3Dcurve_geom() [1/2]

void moab::HiReconstruction::initialize_3Dcurve_geom ( const int  degree )

protected

Definition at line 780 of file HiReconstruction.cpp.

{
 if( !_hasderiv )
 {
 compute_average_vertex_tangents_curve();
 _hasderiv = true;
 }
 if( !_initfittings )
 {
 int ncoeffspvpd = degree + 1;
 _degrees_out.assign( _nv2rec, 0 );
 _interps.assign( _nv2rec, false );
 _vertID2coeffID.resize( _nv2rec );
 _local_fit_coeffs.assign( _nv2rec * ncoeffspvpd * 3, 0 );
 for( size_t i = 0; i < _nv2rec; ++i )
 {
 _vertID2coeffID[i] = i * ncoeffspvpd * 3;
 }
 _initfittings = true;
 }
}

References _degrees_out, _hasderiv, _initfittings, _interps, _local_fit_coeffs, _nv2rec, _vertID2coeffID, and compute_average_vertex_tangents_curve().

Referenced by reconstruct3D_curve_geom().

initialize_3Dcurve_geom() [2/2]

void moab::HiReconstruction::initialize_3Dcurve_geom ( const size_t  npts, const int *  degrees  )

protected

Definition at line 802 of file HiReconstruction.cpp.

{
 if( !_hasderiv )
 {
 compute_average_vertex_tangents_curve();
 _hasderiv = true;
 }
 if( !_hasfittings )
 {
 assert( _nv2rec == npts );
 _degrees_out.assign( _nv2rec, 0 );
 _interps.assign( _nv2rec, false );
 _vertID2coeffID.reserve( _nv2rec );
 size_t index = 0;
 for( size_t i = 0; i < _nv2rec; ++i )
 {
 _vertID2coeffID[i] = index;
 index += 3 * ( degrees[i] + 1 );
 }
 _local_fit_coeffs.assign( index, 0 );
 _initfittings = true;
 }
}

References _degrees_out, _hasderiv, _hasfittings, _initfittings, _interps, _local_fit_coeffs, _nv2rec, _vertID2coeffID, and compute_average_vertex_tangents_curve().

initialize_surf_geom() [1/2]

void moab::HiReconstruction::initialize_surf_geom ( const int  degree )

protected

Initialize the storage for fitting results over _mesh2rec, curve/surface mesh Two options are provided: a) use uniform degree for all vertices b) use customized degrees for different vertices After calling of initializing functions, _initfitting is set to be true, the fitting result could be stored internally

Definition at line 734 of file HiReconstruction.cpp.

{
 if( !_hasderiv )
 {
 compute_average_vertex_normals_surf();
 _hasderiv = true;
 }
 if( !_initfittings )
 {
 int ncoeffspv = ( degree + 2 ) * ( degree + 1 ) / 2;
 _degrees_out.assign( _nv2rec, 0 );
 _interps.assign( _nv2rec, false );
 _vertID2coeffID.resize( _nv2rec );
 _local_fit_coeffs.assign( _nv2rec * ncoeffspv, 0 );
 for( size_t i = 0; i < _nv2rec; ++i )
 {
 _vertID2coeffID[i] = i * ncoeffspv;
 }
 _initfittings = true;
 }
}

References _degrees_out, _hasderiv, _initfittings, _interps, _local_fit_coeffs, _nv2rec, _vertID2coeffID, and compute_average_vertex_normals_surf().

Referenced by reconstruct3D_surf_geom().

initialize_surf_geom() [2/2]

void moab::HiReconstruction::initialize_surf_geom ( const size_t  npts, const int *  degrees  )

protected

Definition at line 756 of file HiReconstruction.cpp.

{
 if( !_hasderiv )
 {
 compute_average_vertex_normals_surf();
 _hasderiv = true;
 }
 if( !_initfittings )
 {
 assert( _nv2rec == npts );
 _degrees_out.assign( _nv2rec, 0 );
 _interps.assign( _nv2rec, false );
 _vertID2coeffID.resize( _nv2rec );
 size_t index = 0;
 for( size_t i = 0; i < _nv2rec; ++i )
 {
 _vertID2coeffID[i] = index;
 index += ( degrees[i] + 2 ) * ( degrees[i] + 1 ) / 2;
 }
 _local_fit_coeffs.assign( index, 0 );
 _initfittings = true;
 }
}

References _degrees_out, _hasderiv, _initfittings, _interps, _local_fit_coeffs, _nv2rec, _vertID2coeffID, and compute_average_vertex_normals_surf().

obtain_nring_ngbvs()

ErrorCode moab::HiReconstruction::obtain_nring_ngbvs ( const EntityHandle  vid, int  ring, const int  minpnts, Range &  ngbvs  )

protected

Get n-ring neighbor vertices, assuming curve/surface mesh, not volume mesh.

Given a vertex, find its n-ring neighbor vertices including itself in _mesrh2rec. 1-ring neighbor vertices of a vertex are the vertices connected with this vertex with an edge n-ring vertices are obtained first get the 1-ring vertices and then get the 1-ring of these vertices, and so on

Parameters

vid	EntityHandle, vertex around which to get n-ring vertices
ring	Integer, number of rings
minpnts	Integer, number of minimum vertices to obtain, if the input ring could not provide enough vertices, i.e. more than minpnts, then expand the number of rings
ngbvs	Range, the n-ring vertices of vid, including vid. If too many points found, i.e. more than _MAXPNTS, then terminate early.

Definition at line 672 of file HiReconstruction.cpp.

{
 ErrorCode error;
 std::deque< EntityHandle > todo;
 todo.push_back( vid );
 ngbvs.insert( vid );
 EntityHandle pre, nxt;
 for( int i = 1; i <= ring; ++i )
 {
 int count = todo.size();
 while( count )
 {
 EntityHandle center = todo.front();
 todo.pop_front();
 --count;
 std::vector< EntityHandle > adjents;
 error = vertex_get_incident_elements( center, _dim, adjents );MB_CHK_ERR( error );
 for( size_t j = 0; j < adjents.size(); ++j )
 {
 std::vector< EntityHandle > elemconn;
 error = mbImpl->get_connectivity( &adjents[j], 1, elemconn );MB_CHK_ERR( error );
 int nvpe = elemconn.size();
 for( int k = 0; k < nvpe; ++k )
 {
 if( elemconn[k] == center )
 {
 pre = k == 0 ? elemconn[nvpe - 1] : elemconn[k - 1];
 nxt = elemconn[( k + 1 ) % nvpe];
 if( ngbvs.find( pre ) == ngbvs.end() )
 {
 ngbvs.insert( pre );
 todo.push_back( pre );
 }
 if( ngbvs.find( nxt ) == ngbvs.end() )
 {
 ngbvs.insert( nxt );
 todo.push_back( nxt );
 }
 break;
 }
 }
 }
 }
 if( _MAXPNTS <= (int)ngbvs.size() )
 {
 // obtain enough points
 return error;
 }
 if( !todo.size() )
 {
 // current ring cannot introduce any points, return incase deadlock
 return error;
 }
 if( ( i == ring ) && ( minpnts > (int)ngbvs.size() ) )
 {
 // reach maximum ring but not enough points
 ++ring;
 }
 }
 return error;
}

References _dim, _MAXPNTS, center(), moab::Range::end(), moab::error(), ErrorCode, moab::Range::find(), moab::Core::get_connectivity(), moab::Range::insert(), MB_CHK_ERR, mbImpl, moab::Range::size(), and vertex_get_incident_elements().

Referenced by polyfit3d_walf_curve_vertex(), and polyfit3d_walf_surf_vertex().

operator=()

HiReconstruction& moab::HiReconstruction::operator= ( const HiReconstruction &  right )

protected

polyfit3d_curve_get_coeff()

void moab::HiReconstruction::polyfit3d_curve_get_coeff ( const int  nverts, const double *  ngbcors, const double *  ngbtangs, int  degree, const bool  interp, const bool  safeguard, const int  ncoords, double *  coords, const int  ncoeffs, double *  coeffs, int *  degree_out  )

protected

Compute local single variable coordinate system of a vertex, and perform vertex based polynomial fittings of three global coordinates axis.

This function take the first vertex of input as center, and build local u-line by estimating tangent vector Then other vertices form vectors originating from center and vectors then are projectd onto this u-plane to form three local height functions, one for each coordinates axis. Local fitting of these local height functions are performed in WLS sense, according if interpolation required or not.

Parameters

nverts	Integer, number of vertices of input
ngbcors	Pointer to array of doubles, size = 3*nverts, coordinates of input vertices, first will be center
ngbtangs	Pointer to array of doubles, size = 3*nverts, vertex tangent vectors of input vertices
degree	Integer, user specified fitting degree
interp	Boolean, user input, interpolation or not
safeguard	Boolean, true = use safeguarded numerical method in computing
ncoords	Integer, size of *coords, should be 3, used for check
coords	Pointer to array of doubles, preallocated memory for storing the glocal coordinates of local u axis direction
ncoeffs	Integer, size of coeffs, should be 3(degree+1), used for check
coeffs	Pointer to array of doubles, preallocated memory for storing coefficients of local fittings in monomial basis
degree_out	Pointer to integer, order of resulted polynomial of fitting, could be downgraded due to numerical issues

Definition at line 1318 of file HiReconstruction.cpp.

{
 if( !nverts )
 {
 return;
 }
 // step 1. compute local coordinates system
 double tang[3] = { ngbtangs[0], ngbtangs[1], ngbtangs[2] };
 if( coords && ncoords > 2 )
 {
 coords[0] = tang[0];
 coords[1] = tang[1];
 coords[2] = tang[2];
 }
 if( !coeffs || !ncoeffs )
 {
 return;
 }
 else
 {
 assert( ncoeffs >= 3 * ( degree + 1 ) );
 }
 // step 2. project onto local coordinate system
 int npts2fit = nverts - interp;
 if( !npts2fit )
 {
 *degree_out = 0;
 for( int i = 0; i < ncoeffs; ++i )
 {
 coeffs[0] = 0;
 }
 return;
 }
 std::vector< double > us( npts2fit ); // double *us = new double[npts2fit];
 std::vector< double > bs( npts2fit * 3 ); // double *bs = new double[npts2fit*3];
 double uu[3];
 for( int i = interp; i < nverts; ++i )
 {
 int k = i - interp;
 DGMSolver::vec_linear_operation( 3, 1, ngbcors + 3 * i, -1, ngbcors, uu );
 us[k] = DGMSolver::vec_innerprod( 3, uu, tang );
 bs[k] = uu[0];
 bs[npts2fit + k] = uu[1];
 bs[2 * npts2fit + k] = uu[2];
 }
 
 // step 3. copmute weights
 std::vector< double > ws( npts2fit );
 int nzeros = compute_weights( npts2fit, 1, &( us[0] ), nverts, ngbtangs, degree, _MINEPS, &( ws[0] ) );
 assert( nzeros <= npts2fit );
 
 // step 4. adjust according to number of zero-weights
 if( nzeros )
 {
 if( nzeros == npts2fit )
 {
 // singular case
 *degree_out = 0;
 for( int i = 0; i < ncoeffs; ++i )
 {
 coeffs[i] = 0;
 }
 return;
 }
 int npts_new = npts2fit - nzeros;
 std::vector< double > bs_temp( 3 * npts_new );
 int index = 0;
 for( int i = 0; i < npts2fit; ++i )
 {
 if( ws[i] )
 {
 if( i > index )
 {
 us[index] = us[i];
 ws[index] = ws[i];
 }
 bs_temp[index] = bs[i];
 bs_temp[index + npts_new] = bs[i + npts2fit];
 bs_temp[index + 2 * npts_new] = bs[i + 2 * npts2fit];
 ++index;
 }
 }
 assert( index == npts_new );
 npts2fit = npts_new;
 us.resize( index );
 ws.resize( index );
 bs = bs_temp;
 // destroy bs_temp;
 std::vector< double >().swap( bs_temp );
 }
 
 // step 5. fitting
 eval_vander_univar_cmf( npts2fit, &( us[0] ), 3, &( bs[0] ), degree, &( ws[0] ), interp, safeguard, degree_out );
 // step 6. write results to output pointers
 int ncoeffs_out_pvpd = *degree_out + 1;
 assert( ncoeffs >= 3 * ncoeffs_out_pvpd );
 // write to coeffs consecutively, bs's size is not changed by eval_vander_univar_cmf
 coeffs[0] = coeffs[ncoeffs_out_pvpd] = coeffs[2 * ncoeffs_out_pvpd] = 0;
 for( int i = 0; i < ncoeffs_out_pvpd - interp; ++i )
 {
 coeffs[i + interp] = bs[i];
 coeffs[i + interp + ncoeffs_out_pvpd] = bs[i + npts2fit];
 coeffs[i + interp + 2 * ncoeffs_out_pvpd] = bs[i + 2 * npts2fit];
 }
}

References _MINEPS, compute_weights(), eval_vander_univar_cmf(), moab::DGMSolver::vec_innerprod(), and moab::DGMSolver::vec_linear_operation().

Referenced by polyfit3d_walf_curve_vertex().

polyfit3d_surf_get_coeff()

void moab::HiReconstruction::polyfit3d_surf_get_coeff ( const int  nverts, const double *  ngbcors, const double *  ngbnrms, int  degree, const bool  interp, const bool  safeguard, const int  ncoords, double *  coords, const int  ncoeffs, double *  coeffs, int *  degree_out, int *  degree_pnt, int *  degree_qr  )

protected

Compute local coordinates system of a vertex, and perform vertex based polynomial fittings of local height function.

This function take the first vertex of input as center, and build local uv-plane by estimating vertex normals and tangent planes Then other vertices forms vectors starting from center and then are projectd onto this uv-plane to form a local height function. Local fitting of this local height function is performed in WLS sense, according if interpolation required or not.

Parameters

nverts	Integer, number of vertices of input
ngbcors	Pointer to array of doubles, size = 3*nverts, coordinates of input vertices, first will be center
ngbnrms	Pointer to array of doubles, size = 3*nverts, vertex normals of input vertices
degree	Integer, user specified fitting degree
interp	Boolean, user input, interpolation or not
safeguard	Boolean, true = use safeguarded numerical method in computing
ncoords	Integer, size of *coords, should be 9, used for check
coords	Pointer to array of doubles, preallocated memory for storing the glocal coordinates of local uvw axis directions
ncoeffs	Integer, size of *coeffs, should be (degree+2)(degree+1)/2, used for check
coeffs	Pointer to array of doubles, preallocated memory for storing coefficients of local fittings in monomial basis
degree_out	Pointer to integer, order of resulted polynomial of fitting, could be downgraded due to numerical issues
degree_pnt	Pointer to integer, polynomial fitting order determined by stencil size/number of points
degree_qr	Pointer to integer, polynomial fitting order determined by Vandermonde system condition number

Definition at line 1033 of file HiReconstruction.cpp.

{
 if( nverts <= 0 )
 {
 return;
 }
 
 // std::cout << "npnts in initial stencil = " << nverts << std::endl;
 // std::cout << "centered at (" << ngbcoords[0] << "," << ngbcoords[1] << "," << ngbcoords[2] <<
 // ")" << std::endl;
 
 // step 1. copmute local coordinate system
 double nrm[3] = { ngbnrms[0], ngbnrms[1], ngbnrms[2] }, tang1[3] = { 0, 0, 0 }, tang2[3] = { 0, 0, 0 };
 if( fabs( nrm[0] ) > fabs( nrm[1] ) && fabs( nrm[0] ) > fabs( nrm[2] ) )
 {
 tang1[1] = 1.0;
 }
 else
 {
 tang1[0] = 1.0;
 }
 
 DGMSolver::vec_projoff( 3, tang1, nrm, tang1 );
#ifndef NDEBUG
 assert( DGMSolver::vec_normalize( 3, tang1, tang1 ) );
#endif
 DGMSolver::vec_crossprod( nrm, tang1, tang2 );
 if( 9 <= ncoords && coords )
 {
 coords[0] = tang1[0];
 coords[1] = tang1[1];
 coords[2] = tang1[2];
 coords[3] = tang2[0];
 coords[4] = tang2[1];
 coords[5] = tang2[2];
 coords[6] = nrm[0];
 coords[7] = nrm[1];
 coords[8] = nrm[2];
 }
 if( !ncoeffs || !coeffs )
 {
 return;
 }
 else
 {
 assert( ncoeffs >= ( degree + 2 ) * ( degree + 1 ) / 2 );
 }
 
 // step 2. project onto local coordinates system
 int npts2fit = nverts - interp;
 if( 0 == npts2fit )
 {
 *degree_out = *degree_pnt = *degree_qr = 0;
 for( int i = 0; i < ncoeffs; ++i )
 {
 coeffs[i] = 0;
 }
 // coeffs[0] = 0;
 return;
 }
 std::vector< double > us( npts2fit * 2 ); // double *us = new double[npts2fit*2];
 std::vector< double > bs( npts2fit ); // double *bs = new double[npts2fit];
 for( int i = interp; i < nverts; ++i )
 {
 int k = i - interp;
 double uu[3];
 DGMSolver::vec_linear_operation( 3, 1, ngbcoords + 3 * i, -1, ngbcoords, uu );
 us[k * 2] = DGMSolver::vec_innerprod( 3, tang1, uu );
 us[k * 2 + 1] = DGMSolver::vec_innerprod( 3, tang2, uu );
 bs[k] = DGMSolver::vec_innerprod( 3, nrm, uu );
 }
 
 // step 3. compute weights
 std::vector< double > ws( npts2fit ); // double *ws = new double[npts2fit];
 int nzeros = compute_weights( npts2fit, 2, &( us[0] ), nverts, ngbnrms, degree, _MINEPS, &( ws[0] ) );
 
 // step 4. adjust according to zero-weights
 if( nzeros )
 {
 if( nzeros == npts2fit )
 {
 *degree_out = *degree_pnt = *degree_qr = 0;
 for( int i = 0; i < ncoeffs; ++i )
 {
 coeffs[i] = 0;
 }
 // coeffs[0] = 0;
 return;
 }
 int index = 0;
 for( int i = 0; i < npts2fit; ++i )
 {
 if( ws[i] )
 {
 if( i > index )
 {
 us[index * 2] = us[i * 2];
 us[index * 2 + 1] = us[i * 2 + 1];
 bs[index] = bs[i];
 ws[index] = ws[i];
 }
 ++index;
 }
 }
 npts2fit -= nzeros;
 assert( index == npts2fit );
 us.resize( npts2fit * 2 );
 bs.resize( npts2fit );
 ws.resize( npts2fit );
 /*us = realloc(us,npts2fit*2*sizeof(double));
 bs = realloc(bs,npts2fit*sizeof(double));
 ws = realloc(ws,npts2fit*sizeof(double));*/
 }
 
 // std::cout<<"npnts after weighting = "<<npts2fit<<std::endl;
 
 // step 5. fitting
 eval_vander_bivar_cmf( npts2fit, &( us[0] ), 1, &( bs[0] ), degree, &( ws[0] ), interp, safeguard, degree_out,
 degree_pnt, degree_qr );
 
 // step 6. organize output
 int ncoeffs_out = ( *degree_out + 2 ) * ( *degree_out + 1 ) / 2;
 assert( ncoeffs_out <= ncoeffs );
 coeffs[0] = 0;
 for( int j = 0; j < ncoeffs_out - interp; ++j )
 {
 coeffs[j + interp] = bs[j];
 }
 // delete [] us; delete [] bs; delete [] ws;
}

References _MINEPS, compute_weights(), eval_vander_bivar_cmf(), moab::DGMSolver::vec_crossprod(), moab::DGMSolver::vec_innerprod(), moab::DGMSolver::vec_linear_operation(), moab::DGMSolver::vec_normalize(), and moab::DGMSolver::vec_projoff().

Referenced by polyfit3d_walf_surf_vertex().

polyfit3d_walf_curve_vertex()

ErrorCode moab::HiReconstruction::polyfit3d_walf_curve_vertex	(	const EntityHandle	vid,
		const bool	interp,
		int	degree,
		int	minpnts,
		const bool	safeguard,
		const int	ncoords,
		double *	coords,
		int *	degree_out,
		const int	ncoeffs,
		double *	coeffs
	)

Construct vertex based polynomial fitting on a curve mesh.

Given a vertex on a curve mesh, construct three one-parameter local fittings for each coordinates axis around this vertex. Stencils around this vertex will be selected according to input degree and if data is noise. Local u-line, or the single parameter will be the estimated tangent line at this vertex. On each axis of xyz, a polynomial fitting will be performed according to user input. minpnts will be used to specify the minimum number allowed in the local stencil. The result will be returned to user by preallocated memory coords, degree_out, coeffs.

Parameters

vid	EntityHandle, the fittings will be performed around this vertex.
interp	Boolean, true=Interpolation, false=least square fitting.
degree	Integer, order of polynomials used for local fittings.
minpnts	Integer, the allowed minimum number of vertices in local stencil. If too small, the resulted fitting might be low order accurate. If too large, it may introduce overfitting.
safeguard	Boolean, true=using safe guarded method in numerical computing.
coords	Pointer to double, preallocated memory by user, should have at least 3 doubles; stores the global coordinates of local coordinate system u direction.
degree_out	Pointer to integer, used to store the degree of resulted fitting
coeffs,Pointer	to double, preallocated memory for coefficients of local fittings, should have at least 3*(degree+1) doubles.

Definition at line 344 of file HiReconstruction.cpp.

{
 ErrorCode error;
 int ring = estimate_num_rings( degree, interp );
 // get n-ring neighbors
 Range ngbvs;
 error = obtain_nring_ngbvs( vid, ring, minpnts, ngbvs );MB_CHK_ERR( error );
 // get coordinates
 size_t nverts = ngbvs.size();
 assert( nverts );
 double* ngbcoords = new double[nverts * 3];
 error = mbImpl->get_coords( ngbvs, ngbcoords );MB_CHK_ERR( error );
 // get tangent vectors
 double* ngbtangs = new double[nverts * 3];
 error = get_tangents_curve( ngbvs, ngbtangs );MB_CHK_ERR( error );
 // switch vid to first one
 int index = ngbvs.index( vid );
 assert( index != -1 );
 std::swap( ngbcoords[0], ngbcoords[3 * index] );
 std::swap( ngbcoords[1], ngbcoords[3 * index + 1] );
 std::swap( ngbcoords[2], ngbcoords[3 * index + 2] );
 std::swap( ngbtangs[0], ngbtangs[3 * index] );
 std::swap( ngbtangs[1], ngbtangs[3 * index + 1] );
 std::swap( ngbtangs[2], ngbtangs[3 * index + 2] );
 // local WLS fittings
 polyfit3d_curve_get_coeff( nverts, ngbcoords, ngbtangs, degree, interp, safeguard, ncoords, coords, ncoeffs, coeffs,
 degree_out );
 delete[] ngbcoords;
 delete[] ngbtangs;
 return error;
}

References moab::error(), ErrorCode, estimate_num_rings(), moab::Core::get_coords(), get_tangents_curve(), moab::Range::index(), MB_CHK_ERR, mbImpl, obtain_nring_ngbvs(), polyfit3d_curve_get_coeff(), and moab::Range::size().

Referenced by reconstruct3D_curve_geom().

polyfit3d_walf_surf_vertex()

ErrorCode moab::HiReconstruction::polyfit3d_walf_surf_vertex	(	const EntityHandle	vid,
		const bool	interp,
		int	degree,
		int	minpnts,
		const bool	safeguard,
		const int	ncoords,
		double *	coords,
		int *	degree_out,
		const int	ncoeffs,
		double *	coeffs
	)

Construct vertex based polynomial fitting on a surface mesh.

Given a vertex on a surface mesh, construct a local fitting around this vertex. Stencils around this vertex will be selected according to input degree and if data is noise. Local uv-plane will be the estimated tangent plane at this vertex. minpnts will be used to specify the minimum number allowed in the local stencil. The result will be returned to user by preallocated memory coords, degree_out, coeffs.

Parameters

vid	EntityHandle, the fitting will be performed around this vertex for the local height function over the uv-plane.
interp	Boolean, true=Interpolation, false=least square fitting.
degree	Integer, order of polynomials used for local fittings.
minpnts	Integer, the allowed minimum number of vertices in local stencil. If too small, the resulted fitting might be low order accurate. If too large, it may introduce overfitting.
safeguard	Boolean, true=using safe guarded method in numerical computing.
coords	Pointer to double, preallocated memory by user, should have at least 9 doubles; stores the global coordinates of local coordinates system uvw directions.
degree_out	Pointer to integer, used to store the degree of resulted fitting
coeffs,Pointer	to double, preallocated memory for coefficients of local fittings, should have at least (degree+2)(degree+1)/2 doubles.

Definition at line 294 of file HiReconstruction.cpp.

{
 assert( _dim == 2 );
 ErrorCode error;
 int ring = estimate_num_rings( degree, interp );
 // std::cout<<"ring = "<<ring<<std::endl;
 // get n-ring neighbors
 Range ngbvs;
 error = obtain_nring_ngbvs( vid, ring, minpnts, ngbvs );MB_CHK_ERR( error );
 // for debug
 /*if(_verts2rec.index(vid)==70){
 for(Range::iterator ingb=ngbvs.begin();ingb!=ngbvs.end();++ingb) std::cerr <<
 _verts2rec.index(*ingb) << " "; std::cout << std::endl;
 }*/
 
 // get coordinates;
 size_t nverts = ngbvs.size();
 assert( nverts );
 double* ngbcoords = new double[nverts * 3];
 error = mbImpl->get_coords( ngbvs, ngbcoords );MB_CHK_ERR( error );
 // get normals
 double* ngbnrms = new double[nverts * 3];
 error = get_normals_surf( ngbvs, ngbnrms );MB_CHK_ERR( error );
 // switch vid to first one
 int index = ngbvs.index( vid );
 assert( index != -1 );
 std::swap( ngbcoords[0], ngbcoords[3 * index] );
 std::swap( ngbcoords[1], ngbcoords[3 * index + 1] );
 std::swap( ngbcoords[2], ngbcoords[3 * index + 2] );
 std::swap( ngbnrms[0], ngbnrms[3 * index] );
 std::swap( ngbnrms[1], ngbnrms[3 * index + 1] );
 std::swap( ngbnrms[2], ngbnrms[3 * index + 2] );
 // local WLS fitting
 int degree_pnt, degree_qr;
 polyfit3d_surf_get_coeff( nverts, ngbcoords, ngbnrms, degree, interp, safeguard, ncoords, coords, ncoeffs, coeffs,
 degree_out, &degree_pnt, &degree_qr );
 delete[] ngbcoords;
 delete[] ngbnrms;
 return error;
}

References _dim, moab::error(), ErrorCode, estimate_num_rings(), moab::Core::get_coords(), get_normals_surf(), moab::Range::index(), MB_CHK_ERR, mbImpl, obtain_nring_ngbvs(), polyfit3d_surf_get_coeff(), and moab::Range::size().

Referenced by reconstruct3D_surf_geom().

reconstruct3D_curve_geom() [1/2]

ErrorCode moab::HiReconstruction::reconstruct3D_curve_geom	(	int	degree,
		bool	interp,
		bool	safeguard,
		bool	reset = false
	)

Reconstruct a high order curve on given curve mesh.

Given a curve mesh, compute vertex based polynomail fittings for all vertices hosted by current processor. The vertex based fitting is done by perfoming three one-parameter fittings along each axis, i.e. x,y,z. The result will be stored interally for later usage of evalution.

Parameters

degree	Integer, order of polynomials used for local fittings.
interp	Boolean, true=Interpolation, false=least square fitting.
safeguard	Boolean, true=using safe guarded method in numerical computing.
reset	Boolean, reset=true will recompute the fittings based on user input and replace the existing one.

Definition at line 221 of file HiReconstruction.cpp.

{
 assert( _dim == 1 );
 if( _hasfittings && !reset )
 {
 return MB_SUCCESS;
 }
 else
 {
 _initfittings = _hasfittings = false;
 }
 initialize_3Dcurve_geom( degree );
 ErrorCode error;
 double *coords = 0, *coeffs;
 int* degree_out;
 int ncoeffs = 3 * ( degree + 1 );
 for( Range::iterator ivert = _verts2rec.begin(); ivert != _verts2rec.end(); ++ivert )
 {
 int index = _verts2rec.index( *ivert );
 assert( index != -1 );
 size_t istr = _vertID2coeffID[index];
 coeffs = &( _local_fit_coeffs[istr] );
 degree_out = &( _degrees_out[index] );
 _interps[index] = interp;
 error = polyfit3d_walf_curve_vertex( *ivert, interp, degree, _MINPNTS, safeguard, 0, coords, degree_out,
 ncoeffs, coeffs );MB_CHK_ERR( error );
 }
 _geom = HI3DCURVE;
 _hasfittings = true;
 return error;
}

References _degrees_out, _dim, _geom, _hasfittings, _initfittings, _interps, _local_fit_coeffs, _MINPNTS, _vertID2coeffID, _verts2rec, moab::Range::begin(), moab::Range::end(), moab::error(), ErrorCode, moab::HI3DCURVE, moab::Range::index(), initialize_3Dcurve_geom(), MB_CHK_ERR, MB_SUCCESS, and polyfit3d_walf_curve_vertex().

reconstruct3D_curve_geom() [2/2]

ErrorCode moab::HiReconstruction::reconstruct3D_curve_geom	(	size_t	npts,
		int *	degrees,
		bool *	interps,
		bool	safeguard,
		bool	reset = false
	)

Reconstruct a high order curve on given curve mesh.

Given a curve mesh, compute vertex based polynomail fittings for all vertices hosted by current processor. The vertex based fitting is done by perfoming three one-parameter fittings along each axis, i.e. x,y,z. User could specify various degrees for different vertices. It assumes that the input degrees for vertices stored in the same order as that this class stores vertices: 1) reconstruction will be only performed at vertices hosted by current processor, thus input npts should match the number of hosted vertices. 2) all hosted vertices will be stored in a MOAB::Range object, degrees for all these vertices should be stored in degrees as the same order in the MOAB::Range object The result will be stored interally for later usage of evalution.

Parameters

npts	Integer size of array pointed by degrees, used for check
degrees	Integer arrray, order of polynomials for local fitting at all hosted vertices.
interp	Boolean, true=Interpolation, false=least square fitting.
safeguard	Boolean, true=using safe guarded method in numerical computing.
reset	Boolean, reset=true will recompute the fittings based on user input and replace the existing one.

Definition at line 253 of file HiReconstruction.cpp.

{
 assert( _dim == 1 );
 ErrorCode error;
 if( npts != _nv2rec )
 {
 MB_SET_ERR( MB_FAILURE, "Input number of degrees doesn't match the number of vertices" );
 }
 if( _hasfittings && !reset )
 {
 return MB_SUCCESS;
 }
 else
 {
 _initfittings = _hasfittings = false;
 }
 // initialize
 initialize_3Dcurve_geom( npts, degrees );
 double *coords = 0, *coeffs;
 int* degree_out;
 size_t i = 0;
 for( Range::iterator ivert = _verts2rec.begin(); ivert != _verts2rec.end(); ++ivert, ++i )
 {
 int index = _verts2rec.index( *ivert );
 size_t istr = _vertID2coeffID[index];
 coeffs = &( _local_fit_coeffs[istr] );
 degree_out = &( _degrees_out[index] );
 _interps[index] = interps[i];
 int ncoeffs = 3 * ( degrees[i] + 1 );
 error = polyfit3d_walf_curve_vertex( *ivert, interps[i], degrees[i], _MINPNTS, safeguard, 0, coords, degree_out,
 ncoeffs, coeffs );MB_CHK_ERR( error );
 }
 _geom = HI3DCURVE;
 _hasfittings = true;
 return error;
}

References _degrees_out, _dim, _geom, _hasfittings, _initfittings, _interps, _local_fit_coeffs, _MINPNTS, _nv2rec, _vertID2coeffID, _verts2rec, moab::Range::begin(), moab::Range::end(), moab::error(), ErrorCode, moab::HI3DCURVE, moab::Range::index(), initialize_3Dcurve_geom(), MB_CHK_ERR, MB_SET_ERR, MB_SUCCESS, and polyfit3d_walf_curve_vertex().

reconstruct3D_surf_geom() [1/2]

ErrorCode moab::HiReconstruction::reconstruct3D_surf_geom	(	int	degree,
		bool	interp,
		bool	safeguard,
		bool	reset = false
	)

Reconstruct a high order surface on given surface mesh.

Given a mesh, compute vertex based polynomial fittings for all vertices hosted by current processor. The result will be stored interally for later usage of evalution. The inputs are: a) degree, which is the order of polynomial used for vertex based fitting. b) interp, if it's true, then interpolation will be applied for local fitting, otherwise it's least square fitting. c) safeguard, specifies whether to use safeguarded numeric method. d) reset, if fittings have been computed and stored in current object, then reset=true will recompute the fittings based on user input and replace the existing one.

Parameters

degree	Integer, order of polynomials used for local fittings.
interp	Boolean, true=Interpolation, false=least square fitting.
safeguard	Boolean, true=using safe guarded method in numerical computing.
reset	Boolean, reset=true will recompute the fittings based on user input and replace the existing one.

Definition at line 126 of file HiReconstruction.cpp.

{
 assert( 2 == _dim );
 if( _hasfittings && !reset )
 {
 // This object has precomputed fitting results and user don't want to reset
 return MB_SUCCESS;
 }
 else
 {
 _initfittings = _hasfittings = false;
 }
 // initialize for geometric information
 initialize_surf_geom( degree );
 ErrorCode error;
 double *coeffs, *coords;
 int* degree_out;
 int ncoeffs = ( degree + 2 ) * ( degree + 1 ) / 2;
 
 // DBG
 int dcount = 0;
 
 for( Range::iterator ivert = _verts2rec.begin(); ivert != _verts2rec.end(); ++ivert )
 {
 int index = _verts2rec.index( *ivert );
 // for debug
 /*if(index==70){
 EntityHandle vid = *ivert;
 double vertcoords[3];
 error = mbImpl->get_coords(&vid,1,vertcoords);
 }*/
 
 size_t istr = _vertID2coeffID[index];
 coords = &( _local_coords[9 * index] );
 coeffs = &( _local_fit_coeffs[istr] );
 degree_out = &( _degrees_out[index] );
 _interps[index] = interp;
 error = polyfit3d_walf_surf_vertex( *ivert, interp, degree, _MINPNTS, safeguard, 9, coords, degree_out, ncoeffs,
 coeffs );MB_CHK_ERR( error );
 
 // DBG
 if( degree_out[0] < degree ) dcount += 1;
 }
 
 // DBG
 // std::cout<<"Total #points ="<<_verts2rec.size()<<", #degraded points = "<<dcount<<std::endl;
 
 _geom = HISURFACE;
 _hasfittings = true;
 return error;
}

References _degrees_out, _dim, _geom, _hasfittings, _initfittings, _interps, _local_coords, _local_fit_coeffs, _MINPNTS, _vertID2coeffID, _verts2rec, moab::Range::begin(), moab::Range::end(), moab::error(), ErrorCode, moab::HISURFACE, moab::Range::index(), initialize_surf_geom(), MB_CHK_ERR, MB_SUCCESS, and polyfit3d_walf_surf_vertex().

reconstruct3D_surf_geom() [2/2]

ErrorCode moab::HiReconstruction::reconstruct3D_surf_geom	(	size_t	npts,
		int *	degrees,
		bool *	interps,
		bool	safeguard,
		bool	reset = false
	)

Reconstruct a high order surface on given surface mesh.

Given a mesh, compute vertex based polynomial fittings for all vertices hosted by current processor. User could specify various degrees for different vertices. It assumes that the input degrees for vertices stored in the same order as that this class stores vertices: 1) reconstruction will be only performed at vertices hosted by current processor, thus input npts should match the number of hosted vertices. 2) all hosted vertices will be stored in a MOAB::Range object, degrees for all these vertices should be stored in degrees as the same order in the MOAB::Range object The result will be stored interally for later usage of evalution.

Parameters

npts	Integer size of array pointed by degrees, used for check
degrees	Integer arrray, order of polynomials for local fitting at all hosted vertices
interp	Boolean, true=Interpolation, false=least square fitting.
safeguard	Boolean, true=using safe guarded method in numerical computing.
reset	Boolean, reset=true will recompute the fittings based on user input and replace the existing one.

Definition at line 178 of file HiReconstruction.cpp.

{
 assert( _dim == 2 );
 if( npts != _nv2rec )
 {
 MB_SET_ERR( MB_FAILURE, "Input number of degrees doesn't match number of vertices" );
 }
 if( _hasfittings && !reset )
 {
 return MB_SUCCESS;
 }
 else
 {
 _initfittings = _hasfittings = false;
 }
 ErrorCode error;
 // initialization for fitting
 initialize_surf_geom( npts, degrees );
 double *coeffs, *coords;
 int* degree_out;
 size_t i = 0;
 for( Range::iterator ivert = _verts2rec.begin(); ivert != _verts2rec.end(); ++ivert, ++i )
 {
 int index = _verts2rec.index( *ivert );
 assert( -1 != index );
 size_t istr = _vertID2coeffID[index];
 coords = &( _local_coords[9 * index] );
 coeffs = &( _local_fit_coeffs[istr] );
 degree_out = &( _degrees_out[index] );
 _interps[index] = interps[i];
 int ncoeffs = ( degrees[i] + 2 ) * ( degrees[i] + 1 ) / 2;
 error = polyfit3d_walf_surf_vertex( *ivert, interps[i], degrees[i], _MINPNTS, safeguard, 9, coords, degree_out,
 ncoeffs, coeffs );MB_CHK_ERR( error );
 }
 _geom = HISURFACE;
 _hasfittings = true;
 return error;
}

References _degrees_out, _dim, _geom, _hasfittings, _initfittings, _interps, _local_coords, _local_fit_coeffs, _MINPNTS, _nv2rec, _vertID2coeffID, _verts2rec, moab::Range::begin(), moab::Range::end(), moab::error(), ErrorCode, moab::HISURFACE, moab::Range::index(), initialize_surf_geom(), MB_CHK_ERR, MB_SET_ERR, MB_SUCCESS, and polyfit3d_walf_surf_vertex().

vertex_get_incident_elements()

ErrorCode moab::HiReconstruction::vertex_get_incident_elements ( const EntityHandle &  vid, const int  elemdim, std::vector< EntityHandle > &  adjents  )

protected

Given a vertex, return the incident elements with dimension elemdim.

Wrapper of MOAB Core->get_adjacencies and HalfRep->get_up_adjacencies, depends on if USE_AHF is defined

Parameters

vid	EntityHandle of vertex
elemdim	Integer, dimension of elements incidented in vid
adjents	vector<EntityHandle>, container which push incident elements in

Definition at line 658 of file HiReconstruction.cpp.

{
 ErrorCode error;
 assert( elemdim == _dim );
#ifdef HIREC_USE_AHF
 error = ahf->get_up_adjacencies( vid, elemdim, adjents );MB_CHK_ERR( error );
#else
 error = mbImpl->get_adjacencies( &vid, 1, elemdim, false, adjents );MB_CHK_ERR( error );
#endif
 return error;
}

References _dim, ahf, moab::error(), ErrorCode, moab::Core::get_adjacencies(), moab::HalfFacetRep::get_up_adjacencies(), MB_CHK_ERR, and mbImpl.

Referenced by average_vertex_normal(), average_vertex_tangent(), and obtain_nring_ngbvs().

walf3d_curve_vertex_eval()

void moab::HiReconstruction::walf3d_curve_vertex_eval	(	const double *	local_origin,
		const double *	local_coords,
		const int	local_deg,
		const double *	local_coeffs,
		const bool	interp,
		const int	npts2fit,
		const double *	coords2fit,
		double *	hiproj_new
	)

Perform high order projection of points around a center vertex, assume geometry is curve.

Given a vertex position and the local one-parameter fittings parameter around this vertex, estimate the ideal position of input position according to the local fittings. This is done by first projecting input points onto the local u-direction at this vertex and then use the value u as parameter for the three fittings, one for each coordinates axis of xyz. The result will be returned to user preallocated memory

Parameters

local_origin	Pointer to 3 doubles, coordinates of the center vertex
local_coords	Pointer to 3 doubles, global coordinates of direction of local u coordinate axis at center vertex
local_deg	Integer, order of local polynomial fitting
local_coeffs	Pointer to array of doubles, size=3*(local_deg+1), coefficients of three local polynomial fittings, in monomial basis. For each fitting, local_deg+1 parameters.
interp	Boolean, true=local fitting is interpolation, false=local fitting is least square fitting
npts2fit	Integer, number of points to be estimated, around the center vertices
coords2fit	Pointer to array of doubles, size=3*npts2fit, current coordinates of points to be estimated
hiproj_new	Pointer to array of doubles, size=3*npts2fit, memory preallocated by user to store the fitting/estimated positions of input points.

Definition at line 569 of file HiReconstruction.cpp.

{
 assert( local_origin && local_coords && local_coeffs );
 int ncoeffspvpd = local_deg + 1;
 for( int i = 0; i < npts2fit; ++i )
 {
 // get the vector from center to current point, project to tangent line
 double vec[3], ans[3] = { 0, 0, 0 };
 DGMSolver::vec_linear_operation( 3, 1, coords2fit + 3 * i, -1, local_origin, vec );
 double u = DGMSolver::vec_innerprod( 3, local_coords, vec );
 // evaluate polynomials
 if( !interp )
 {
 ans[0] = local_coeffs[0];
 ans[1] = local_coeffs[ncoeffspvpd];
 ans[2] = local_coeffs[2 * ncoeffspvpd];
 }
 double uk = 1; // degree_out and degree different, stored in columnwise contiguously
 for( int j = 1; j < ncoeffspvpd; ++j )
 {
 uk *= u;
 ans[0] += uk * local_coeffs[j];
 ans[1] += uk * local_coeffs[j + ncoeffspvpd];
 ans[2] += uk * local_coeffs[j + 2 * ncoeffspvpd];
 }
 hiproj_new[3 * i] = ans[0] + local_origin[0];
 hiproj_new[3 * i + 1] = ans[1] + local_origin[1];
 hiproj_new[3 * i + 2] = ans[2] + local_origin[2];
 }
}

References moab::DGMSolver::vec_innerprod(), and moab::DGMSolver::vec_linear_operation().

Referenced by hiproj_walf_around_vertex().

walf3d_surf_vertex_eval()

void moab::HiReconstruction::walf3d_surf_vertex_eval	(	const double *	local_origin,
		const double *	local_coords,
		const int	local_deg,
		const double *	local_coeffs,
		const bool	interp,
		const int	npts2fit,
		const double *	coords2fit,
		double *	hiproj_new
	)

Perform high order projection of points around a center vertex, assume geometry is surface.

Given a vertex position and the local fitting parameter around this vertex, estimate the ideal position of input position according to the local fitting. This is done by first projecting input points onto the local uv-plane around this vertex and use the given fitting to estimate the ideal position of input points. The result will be returned to user preallocated memory

Parameters

local_origin	Pointer to 3 doubles, coordinates of the center vertex
local_coords	Pointer to 9 doubles, global coordinates of directions of local uvw coordinates axis at center vertex
local_deg	Integer, order of local polynomial fitting
local_coeffs	Pointer to array of doubles, size=(local_deg+2)(local_deg+1)/2, coefficients of local polynomial fittings, in monomial basis
interp	Boolean, true=local fitting is interpolation, false=local fitting is least square fitting
npts2fit	Integer, number of points to be estimated, around the center vertices
coords2fit	Pointer to array of doubles, size=3*npts2fit, current coordinates of points to be estimated
hiproj_new	Pointer to array of doubles, size=3*npts2fit, memory preallocated by user to store the fitting/estimated positions of input points.

Definition at line 513 of file HiReconstruction.cpp.

{
 double xaxis[3], yaxis[3], zaxis[3];
 for( int i = 0; i < 3; ++i )
 {
 xaxis[i] = local_coords[i];
 yaxis[i] = local_coords[3 + i];
 zaxis[i] = local_coords[6 + i];
 }
 // double *basis = new double[(local_deg+2)*(local_deg+1)/2-1];
 std::vector< double > basis( ( local_deg + 2 ) * ( local_deg + 1 ) / 2 - 1 );
 for( int i = 0; i < npts2fit; ++i )
 {
 double local_pos[3];
 for( int j = 0; j < 3; ++j )
 {
 local_pos[j] = coords2fit[3 * i + j] - local_origin[j];
 }
 double u, v, height = 0;
 u = DGMSolver::vec_innerprod( 3, local_pos, xaxis );
 v = DGMSolver::vec_innerprod( 3, local_pos, yaxis );
 basis[0] = u;
 basis[1] = v;
 int l = 1;
 for( int k = 2; k <= local_deg; ++k )
 {
 ++l;
 basis[l] = u * basis[l - k];
 for( int id = 0; id < k; ++id )
 {
 ++l;
 basis[l] = basis[l - k - 1] * v;
 }
 }
 if( !interp )
 {
 height = local_coeffs[0];
 }
 for( int p = 0; p <= l; ++p )
 {
 height += local_coeffs[p + 1] * basis[p];
 }
 hiproj_new[3 * i] = local_origin[0] + u * xaxis[0] + v * yaxis[0] + height * zaxis[0];
 hiproj_new[3 * i + 1] = local_origin[1] + u * xaxis[1] + v * yaxis[1] + height * zaxis[1];
 hiproj_new[3 * i + 2] = local_origin[2] + u * xaxis[2] + v * yaxis[2] + height * zaxis[2];
 }
 // delete [] basis;
}

References moab::DGMSolver::vec_innerprod().

Referenced by hiproj_walf_around_vertex().

Member Data Documentation

_degrees_out

std::vector< int > moab::HiReconstruction::_degrees_out

protected

Definition at line 312 of file HiReconstruction.hpp.

Referenced by get_fittings_data(), hiproj_walf_around_vertex(), initialize_3Dcurve_geom(), initialize_surf_geom(), reconstruct3D_curve_geom(), and reconstruct3D_surf_geom().

_dim

int moab::HiReconstruction::_dim

protected

Definition at line 306 of file HiReconstruction.hpp.

Referenced by HiReconstruction(), initialize(), obtain_nring_ngbvs(), polyfit3d_walf_surf_vertex(), reconstruct3D_curve_geom(), reconstruct3D_surf_geom(), and vertex_get_incident_elements().

_geom

GEOMTYPE moab::HiReconstruction::_geom

protected

Definition at line 305 of file HiReconstruction.hpp.

Referenced by get_fittings_data(), hiproj_walf_around_vertex(), reconstruct3D_curve_geom(), and reconstruct3D_surf_geom().

_hasderiv

bool moab::HiReconstruction::_hasderiv

protected

Definition at line 303 of file HiReconstruction.hpp.

Referenced by compute_average_vertex_normals_surf(), compute_average_vertex_tangents_curve(), get_normals_surf(), get_tangents_curve(), HiReconstruction(), initialize(), initialize_3Dcurve_geom(), and initialize_surf_geom().

_hasfittings

bool moab::HiReconstruction::_hasfittings

protected

Definition at line 307 of file HiReconstruction.hpp.

Referenced by get_fittings_data(), hiproj_walf_around_vertex(), hiproj_walf_in_element(), HiReconstruction(), initialize_3Dcurve_geom(), reconstruct3D_curve_geom(), and reconstruct3D_surf_geom().

_incells

Range moab::HiReconstruction::_incells

protected

Definition at line 295 of file HiReconstruction.hpp.

Referenced by initialize().

_inedges

Range moab::HiReconstruction::_inedges

protected

Definition at line 295 of file HiReconstruction.hpp.

Referenced by initialize().

_infaces

Range moab::HiReconstruction::_infaces

protected

Definition at line 295 of file HiReconstruction.hpp.

Referenced by initialize().

_initfittings

bool moab::HiReconstruction::_initfittings

protected

Definition at line 308 of file HiReconstruction.hpp.

Referenced by initialize_3Dcurve_geom(), initialize_surf_geom(), reconstruct3D_curve_geom(), and reconstruct3D_surf_geom().

_interps

std::vector< bool > moab::HiReconstruction::_interps

protected

Definition at line 313 of file HiReconstruction.hpp.

Referenced by get_fittings_data(), hiproj_walf_around_vertex(), initialize_3Dcurve_geom(), initialize_surf_geom(), reconstruct3D_curve_geom(), and reconstruct3D_surf_geom().

_inverts

Range moab::HiReconstruction::_inverts

protected

Definition at line 295 of file HiReconstruction.hpp.

Referenced by initialize().

_local_coords

std::vector< double > moab::HiReconstruction::_local_coords

protected

Definition at line 309 of file HiReconstruction.hpp.

Referenced by compute_average_vertex_normals_surf(), compute_average_vertex_tangents_curve(), get_fittings_data(), get_normals_surf(), get_tangents_curve(), hiproj_walf_around_vertex(), and reconstruct3D_surf_geom().

_local_fit_coeffs

std::vector< double > moab::HiReconstruction::_local_fit_coeffs

protected

Definition at line 310 of file HiReconstruction.hpp.

Referenced by get_fittings_data(), hiproj_walf_around_vertex(), initialize_3Dcurve_geom(), initialize_surf_geom(), reconstruct3D_curve_geom(), and reconstruct3D_surf_geom().

_MAXPNTS

int moab::HiReconstruction::_MAXPNTS

protected

Definition at line 298 of file HiReconstruction.hpp.

Referenced by compute_weights(), initialize(), and obtain_nring_ngbvs().

_mesh2rec

EntityHandle moab::HiReconstruction::_mesh2rec

protected

Definition at line 292 of file HiReconstruction.hpp.

Referenced by initialize().

_MINEPS

double moab::HiReconstruction::_MINEPS

protected

Definition at line 299 of file HiReconstruction.hpp.

Referenced by check_barycentric_coords(), HiReconstruction(), polyfit3d_curve_get_coeff(), and polyfit3d_surf_get_coeff().

_MINPNTS

int moab::HiReconstruction::_MINPNTS

protected

Definition at line 298 of file HiReconstruction.hpp.

Referenced by reconstruct3D_curve_geom(), and reconstruct3D_surf_geom().

_nv2rec

size_t moab::HiReconstruction::_nv2rec

protected

Definition at line 296 of file HiReconstruction.hpp.

Referenced by compute_average_vertex_normals_surf(), compute_average_vertex_tangents_curve(), initialize(), initialize_3Dcurve_geom(), initialize_surf_geom(), reconstruct3D_curve_geom(), and reconstruct3D_surf_geom().

_vertID2coeffID

std::vector< size_t > moab::HiReconstruction::_vertID2coeffID

protected

Definition at line 311 of file HiReconstruction.hpp.

Referenced by get_fittings_data(), hiproj_walf_around_vertex(), initialize_3Dcurve_geom(), initialize_surf_geom(), reconstruct3D_curve_geom(), and reconstruct3D_surf_geom().

_verts2rec

Range moab::HiReconstruction::_verts2rec

protected

Definition at line 295 of file HiReconstruction.hpp.

Referenced by compute_average_vertex_normals_surf(), compute_average_vertex_tangents_curve(), get_fittings_data(), get_normals_surf(), get_tangents_curve(), hiproj_walf_around_vertex(), hiproj_walf_in_element(), initialize(), reconstruct3D_curve_geom(), and reconstruct3D_surf_geom().

ahf

HalfFacetRep* moab::HiReconstruction::ahf

protected

Definition at line 286 of file HiReconstruction.hpp.

Referenced by initialize(), vertex_get_incident_elements(), and ~HiReconstruction().

mbImpl

Core* moab::HiReconstruction::mbImpl

protected

Definition at line 284 of file HiReconstruction.hpp.

Referenced by average_vertex_normal(), average_vertex_tangent(), hiproj_walf_around_vertex(), hiproj_walf_in_element(), HiReconstruction(), initialize(), obtain_nring_ngbvs(), polyfit3d_walf_curve_vertex(), polyfit3d_walf_surf_vertex(), and vertex_get_incident_elements().

pcomm

ParallelComm* moab::HiReconstruction::pcomm

protected

Definition at line 285 of file HiReconstruction.hpp.

Referenced by HiReconstruction(), and initialize().

---

The documentation for this class was generated from the following files:

• HiReconstruction.hpp
• HiReconstruction.cpp