#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.
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.
ErrorCode	reconstruct3D_curve_geom (int degree, bool interp, bool safeguard, bool reset=false)
	Reconstruct a high order curve on given curve mesh.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
bool	check_barycentric_coords (const int nws, const double *naturalcoords)
	Check the correctness of barycentric coordination, wi>=0 and sum(wi)=1.
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

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.

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

    : 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 );
    }
}

moab::HiReconstruction::~HiReconstruction ( )

Definition at line 42 of file HiReconstruction.cpp.

References ahf.

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

moab::HiReconstruction::HiReconstruction ( const HiReconstruction & source ) [protected]

Member Function Documentation

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.

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().

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

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.

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().

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

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.

References _MINEPS, and moab::sum().

Referenced by hiproj_walf_in_element().

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

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.

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().

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

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.

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().

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

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.

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

Referenced by polyfit3d_curve_get_coeff(), and polyfit3d_surf_get_coeff().

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

static int moab::HiReconstruction::estimate_num_ghost_layers	(	int	degree,
		bool	interp = `false`
	)		`[inline, static]`

Definition at line 277 of file HiReconstruction.hpp.

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

int moab::HiReconstruction::estimate_num_rings	(	int	degree,
		bool	interp
	)		`[protected]`

Definition at line 653 of file HiReconstruction.cpp.

Referenced by polyfit3d_walf_curve_vertex(), and polyfit3d_walf_surf_vertex().

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

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.

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().

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

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.

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().

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

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.

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

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

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.

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().

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

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.

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().

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

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.

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().

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

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.

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.

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

ErrorCode moab::HiReconstruction::initialize ( bool recwhole )

Definition at line 51 of file HiReconstruction.cpp.

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().

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

void moab::HiReconstruction::initialize_3Dcurve_geom ( const int degree ) [protected]

Definition at line 780 of file HiReconstruction.cpp.

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

Referenced by reconstruct3D_curve_geom().

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

void moab::HiReconstruction::initialize_3Dcurve_geom	(	const size_t	npts,
		const int *	degrees
	)		`[protected]`

Definition at line 802 of file HiReconstruction.cpp.

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

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

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.

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

Referenced by reconstruct3D_surf_geom().

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

void moab::HiReconstruction::initialize_surf_geom	(	const size_t	npts,
		const int *	degrees
	)		`[protected]`

Definition at line 756 of file HiReconstruction.cpp.

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

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

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.

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().

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

HiReconstruction& moab::HiReconstruction::operator= ( const HiReconstruction & right ) [protected]

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.

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

Referenced by polyfit3d_walf_curve_vertex().

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

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.

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().

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

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.

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().

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

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.

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().

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

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.

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().

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

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.

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().

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

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.

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().

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

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.

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().

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

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.

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().

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

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.

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

Referenced by hiproj_walf_around_vertex().

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

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.

References moab::DGMSolver::vec_innerprod().

Referenced by hiproj_walf_around_vertex().

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