Mesh Oriented datABase  (version 5.6.0)
An array-based unstructured mesh library
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)
 
HiReconstructionoperator= (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

CorembImpl
 
ParallelCommpcomm
 
HalfFacetRepahf
 
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.

19  : mbImpl( impl ), pcomm( comm ), _mesh2rec( meshIn ), _MINPNTS( minpnts )
20 {
21  assert( NULL != impl );
23  _MINEPS = 1e-12;
24  _dim = 0;
25  _hasfittings = false;
26  _hasderiv = false;
27 #ifdef MOAB_HAVE_MPI
28  if( !pcomm )
29  {
31  }
32 #endif
33 
34  error = initialize( recwhole );
35  if( MB_SUCCESS != error )
36  {
37  std::cout << "Error initializing HiReconstruction\n" << std::endl;
38  exit( 1 );
39  }
40 }

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.

43 {
44 #ifdef MOAB_HAVE_AHF
45  ahf = NULL;
46 #else
47  delete ahf;
48 #endif
49 }

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
vidEntityHandle, a vertex in _mesh2rec, in _verts2rec
coordsPointer to double array, global coordinates of local uvw coordinate system axis directions
degree_outInteger, order of polynomial fittings for vid
coeffsPointer to double array, coefficients of local polynomial fittings in monomial basis
interpBoolean, 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
vidEntityHandle, vertex in _mesh2rec, might be ghost vertex
nrmPointer to 3-doubles array, preallocated by user

Definition at line 865 of file HiReconstruction.cpp.

866 {
868  std::vector< EntityHandle > adjfaces;
869  error = vertex_get_incident_elements( vid, 2, adjfaces );
870  MB_CHK_ERR( error );
871  int npolys = adjfaces.size();
872  if( !npolys )
873  {
874  MB_SET_ERR( MB_FAILURE, "Vertex has no incident 2D entities" );
875  }
876  else
877  {
878  double v1[3], v2[3], v3[3], a[3], b[3], c[3];
879  nrm[0] = nrm[1] = nrm[2] = 0;
880  for( int i = 0; i < npolys; ++i )
881  {
882  // get incident "triangles"
883  std::vector< EntityHandle > elemconn;
884  error = mbImpl->get_connectivity( &adjfaces[i], 1, elemconn );
885  MB_CHK_ERR( error );
886  EntityHandle pre, nxt;
887  int nvpe = elemconn.size();
888  for( int j = 0; j < nvpe; ++j )
889  {
890  if( vid == elemconn[j] )
891  {
892  pre = j == 0 ? elemconn[nvpe - 1] : elemconn[j - 1];
893  nxt = elemconn[( j + 1 ) % nvpe];
894  break;
895  }
896  }
897  // compute area weighted normals
898  error = mbImpl->get_coords( &pre, 1, a );
899  MB_CHK_ERR( error );
900  error = mbImpl->get_coords( &vid, 1, b );
901  MB_CHK_ERR( error );
902  error = mbImpl->get_coords( &nxt, 1, c );
903  MB_CHK_ERR( error );
904  DGMSolver::vec_linear_operation( 3, 1, c, -1, b, v1 );
905  DGMSolver::vec_linear_operation( 3, 1, a, -1, b, v2 );
906  DGMSolver::vec_crossprod( v1, v2, v3 );
907  DGMSolver::vec_linear_operation( 3, 1, nrm, 1, v3, nrm );
908  }
909 #ifndef NDEBUG
910  assert( DGMSolver::vec_normalize( 3, nrm, nrm ) );
911 #endif
912  }
913  return error;
914 }

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
vidEntityHandle, vertex in _mesh2rec, might be ghost vertex
tangPointer to 3-doubles array, preallocated by user

Definition at line 975 of file HiReconstruction.cpp.

976 {
978  std::vector< EntityHandle > adjedges;
979  error = vertex_get_incident_elements( vid, 1, adjedges );
980  MB_CHK_ERR( error );
981  int nedges = adjedges.size();
982  if( !nedges )
983  {
984  MB_SET_ERR( MB_FAILURE, "Vertex has no incident edges" );
985  }
986  else
987  {
988  assert( nedges <= 2 );
989  tang[0] = tang[1] = tang[2] = 0;
990  for( int i = 0; i < nedges; ++i )
991  {
992  std::vector< EntityHandle > edgeconn;
993  error = mbImpl->get_connectivity( &adjedges[i], 1, edgeconn );
994  double istr[3], iend[3], t[3];
995  error = mbImpl->get_coords( &( edgeconn[0] ), 1, istr );
996  error = mbImpl->get_coords( &( edgeconn[1] ), 1, iend );
997  DGMSolver::vec_linear_operation( 3, 1, iend, -1, istr, t );
998  DGMSolver::vec_linear_operation( 3, 1, tang, 1, t, tang );
999  }
1000 #ifndef NDEBUG
1001  assert( DGMSolver::vec_normalize( 3, tang, tang ) );
1002 #endif
1003  }
1004  return error;
1005 }

References moab::error(), ErrorCode, moab::Core::get_connectivity(), moab::Core::get_coords(), MB_CHK_ERR, MB_SET_ERR, mbImpl, 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 1629 of file HiReconstruction.cpp.

1630 {
1631  double sum = 0;
1632  for( int i = 0; i < nws; ++i )
1633  {
1634  if( naturalcoords[i] < -_MINEPS )
1635  {
1636  return false;
1637  }
1638  sum += naturalcoords[i];
1639  }
1640  if( fabs( 1 - sum ) > _MINEPS )
1641  {
1642  return false;
1643  }
1644  else
1645  {
1646  return true;
1647  }
1648 }

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 916 of file HiReconstruction.cpp.

917 {
918  if( _hasderiv )
919  {
920  return MB_SUCCESS;
921  }
923  _local_coords.assign( 9 * _nv2rec, 0 );
924  size_t index = 0;
925  for( Range::iterator ivert = _verts2rec.begin(); ivert != _verts2rec.end(); ++ivert, ++index )
926  {
927  error = average_vertex_normal( *ivert, &( _local_coords[9 * index + 6] ) );
928  MB_CHK_ERR( error );
929  }
930  return error;
931 }

References _hasderiv, _local_coords, _nv2rec, _verts2rec, average_vertex_normal(), moab::Range::begin(), moab::Range::end(), moab::error(), ErrorCode, moab::index, 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 1007 of file HiReconstruction.cpp.

1008 {
1009  if( _hasderiv )
1010  {
1011  return MB_SUCCESS;
1012  }
1013  ErrorCode error;
1014  _local_coords.assign( 3 * _nv2rec, 0 );
1015  size_t index = 0;
1016  for( Range::iterator ivert = _verts2rec.begin(); ivert != _verts2rec.end(); ++ivert, ++index )
1017  {
1018  error = average_vertex_tangent( *ivert, &( _local_coords[3 * index] ) );
1019  MB_CHK_ERR( error );
1020  }
1021  return error;
1022 }

References _hasderiv, _local_coords, _nv2rec, _verts2rec, average_vertex_tangent(), moab::Range::begin(), moab::Range::end(), moab::error(), ErrorCode, moab::index, 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 1580 of file HiReconstruction.cpp.

1588 {
1589  assert( nrows <= _MAXPNTS && ws );
1590  bool interp = false;
1591  if( nngbs != nrows )
1592  {
1593  assert( nngbs == nrows + 1 );
1594  interp = true;
1595  }
1596  double epsilon = 1e-2;
1597 
1598  // First, compute squared distance from each input piont to the center
1599  for( int i = 0; i < nrows; ++i )
1600  {
1601  ws[i] = DGMSolver::vec_innerprod( ncols, us + i * ncols, us + i * ncols );
1602  }
1603 
1604  // Second, compute a small correction termt o guard against zero
1605  double h = 0;
1606  for( int i = 0; i < nrows; ++i )
1607  {
1608  h += ws[i];
1609  }
1610  h /= (double)nrows;
1611 
1612  // Finally, compute the weights for each vertex
1613  int nzeros = 0;
1614  for( int i = 0; i < nrows; ++i )
1615  {
1616  double costheta = DGMSolver::vec_innerprod( 3, ngbnrms, ngbnrms + 3 * ( i + interp ) );
1617  if( costheta > toler )
1618  {
1619  ws[i] = costheta * pow( ws[i] / h + epsilon, -1 * (double)degree / 2.0 );
1620  }
1621  else
1622  {
1623  ws[i] = 0;
1624  ++nzeros;
1625  }
1626  }
1627  return nzeros;
1628 }

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.

278  {
279  return 1 + ( interp ? ( ( degree + 1 ) >> 1 ) + ( ( degree + 1 ) & 1 )
280  : ( ( degree + 2 ) >> 1 ) + ( ( degree + 2 ) & 1 ) );
281  };

◆ estimate_num_rings()

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

Definition at line 673 of file HiReconstruction.cpp.

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

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 1212 of file HiReconstruction.cpp.

1223 {
1224  // step 1. adjust the degree according to number of points to fit
1225  int ncols = ( ( ( degree + 2 ) * ( degree + 1 ) ) >> 1 ) - interp;
1226  while( 1 < degree && npts2fit < ncols )
1227  {
1228  --degree;
1229  ncols = ( ( ( degree + 2 ) * ( degree + 1 ) ) >> 1 ) - interp;
1230  }
1231  *degree_pnt = degree;
1232 
1233  // std::cout << "degree_pnt: " << *degree_pnt << std::endl;
1234 
1235  // step 2. construct Vandermonde matrix, stored in columnwise
1236  std::vector< double > V; // V(npts2fit*(ncols+interp)); //double *V_init = new double[npts2fit*(ncols+interp)];
1237  DGMSolver::gen_vander_multivar( npts2fit, 2, us, degree, V );
1238  // remove the first column of 1s if interpolation
1239  if( interp )
1240  {
1241  V.erase( V.begin(), V.begin() + npts2fit );
1242  }
1243  /*double* V;
1244  if(interp){
1245  V = new double[npts2fit*ncols];
1246  std::memcpy(V,V_init+npts2fit,ncols*npts2fit*sizeof(double));
1247  delete [] V_init; V_init = 0;
1248  }else{
1249  V = V_init;
1250  }*/
1251 
1252  // step 3. Scale rows to assign different weights to different points
1253  for( int i = 0; i < npts2fit; ++i )
1254  {
1255  for( int j = 0; j < ncols; ++j )
1256  {
1257  V[j * npts2fit + i] *= ws[i];
1258  }
1259  for( int k = 0; k < ndim; ++k )
1260  {
1261  bs[k * npts2fit + i] *= ws[i];
1262  }
1263  }
1264 
1265  // step 4. scale columns to reduce condition number
1266  std::vector< double > ts( ncols ); // double *ts = new double[ncols];
1267  DGMSolver::rescale_matrix( npts2fit, ncols, &( V[0] ), &( ts[0] ) );
1268 
1269  // step 5. Perform Householder QR factorization
1270  std::vector< double > D( ncols ); // double *D = new double[ncols];
1271  int rank;
1272  DGMSolver::qr_polyfit_safeguarded( npts2fit, ncols, &( V[0] ), &( D[0] ), &rank );
1273 
1274  // step 6. adjust degree of fitting according to rank of Vandermonde matrix
1275  int ncols_sub = ncols;
1276  while( rank < ncols_sub )
1277  {
1278  --degree;
1279  if( degree == 0 )
1280  {
1281  // surface is flat, return 0
1282  *degree_out = *degree_qr = degree;
1283  for( int i = 0; i < npts2fit; ++i )
1284  {
1285  for( int k = 0; k < ndim; ++k )
1286  {
1287  bs[k * npts2fit + i] = 0;
1288  }
1289  }
1290  return;
1291  }
1292  else
1293  {
1294  ncols_sub = ( ( ( degree + 2 ) * ( degree + 1 ) ) >> 1 ) - interp;
1295  }
1296  }
1297  *degree_qr = degree;
1298 
1299  // std::cout << "degree_qr: " << *degree_qr << std::endl;
1300 
1301  /* DBG
1302  * std::cout<<"before Qtb"<<std::endl;
1303  std::cout<<std::endl;
1304  std::cout<<"bs = "<<std::endl;
1305  std::cout<<std::endl;
1306  for (int k=0; k< ndim; k++){
1307  for (int j=0; j<npts2fit; ++j){
1308  std::cout<<" "<<bs[npts2fit*k+j]<<std::endl;
1309  }
1310  }
1311  std::cout<<std::endl;*/
1312 
1313  // step 7. compute Q'b
1314  DGMSolver::compute_qtransposeB( npts2fit, ncols_sub, &( V[0] ), ndim, bs );
1315 
1316  /* DBG
1317  * std::cout<<"after Qtb"<<std::endl;
1318  std::cout<<"bs = "<<std::endl;
1319  std::cout<<std::endl;
1320  for (int k=0; k< ndim; k++){
1321  for (int j=0; j<npts2fit; ++j){
1322  std::cout<<" "<<bs[npts2fit*k+j]<<std::endl;
1323  }
1324  }
1325  std::cout<<std::endl;*/
1326 
1327  // step 8. perform backward substitution and scale the solution
1328  // assign diagonals of V
1329  for( int i = 0; i < ncols_sub; ++i )
1330  {
1331  V[i * npts2fit + i] = D[i];
1332  }
1333 
1334  // backsolve
1335  if( safeguard )
1336  {
1337  // for debug
1338  // std::cout << "ts size " << ts.size() << std::endl;
1339  DGMSolver::backsolve_polyfit_safeguarded( 2, degree, interp, npts2fit, ncols_sub, &( V[0] ), ndim, bs,
1340  &( ts[0] ), degree_out );
1341  }
1342  else
1343  {
1344  DGMSolver::backsolve( npts2fit, ncols_sub, &( V[0] ), 1, bs, &( ts[0] ) );
1345  *degree_out = degree;
1346  }
1347  /*if(V_init){
1348  delete [] V_init;
1349  }else{
1350  delete [] V;
1351  }*/
1352 }

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 1470 of file HiReconstruction.cpp.

1479 {
1480  // step 1. determine degree of polynomials to fit according to number of points
1481  int ncols = degree + 1 - interp;
1482  while( npts2fit < ncols && degree >= 1 )
1483  {
1484  --degree;
1485  ncols = degree + 1 - interp;
1486  }
1487  if( !degree )
1488  {
1489  if( interp )
1490  {
1491  for( int icol = 0; icol < ndim; ++icol )
1492  {
1493  bs[icol * npts2fit] = 0;
1494  }
1495  }
1496  for( int irow = 1; irow < npts2fit; ++irow )
1497  {
1498  for( int icol = 0; icol < ndim; ++icol )
1499  {
1500  bs[icol * npts2fit + irow] = 0;
1501  }
1502  }
1503  *degree_out = 0;
1504  return;
1505  }
1506  // step 2. construct Vandermonde matrix
1507  std::vector< double > V; // V(npts2fit*(ncols+interp));
1508  DGMSolver::gen_vander_multivar( npts2fit, 1, us, degree, V );
1509 
1510  if( interp )
1511  {
1512  V.erase( V.begin(), V.begin() + npts2fit );
1513  }
1514 
1515  // step 3. scale rows with respect to weights
1516  for( int i = 0; i < npts2fit; ++i )
1517  {
1518  for( int j = 0; j < ncols; ++j )
1519  {
1520  V[j * npts2fit + i] *= ws[i];
1521  }
1522  for( int k = 0; k < ndim; ++k )
1523  {
1524  bs[k * npts2fit + i] *= ws[i];
1525  }
1526  }
1527 
1528  // step 4. scale columns to reduce condition number
1529  std::vector< double > ts( ncols );
1530  DGMSolver::rescale_matrix( npts2fit, ncols, &( V[0] ), &( ts[0] ) );
1531 
1532  // step 5. perform Householder QR factorization
1533  std::vector< double > D( ncols );
1534  int rank;
1535  DGMSolver::qr_polyfit_safeguarded( npts2fit, ncols, &( V[0] ), &( D[0] ), &rank );
1536 
1537  // step 6. adjust degree of fitting
1538  int ncols_sub = ncols;
1539  while( rank < ncols_sub )
1540  {
1541  --degree;
1542  if( !degree )
1543  {
1544  // matrix is singular, consider curve on flat plane passing center and orthogonal to
1545  // tangent line
1546  *degree_out = 0;
1547  for( int i = 0; i < npts2fit; ++i )
1548  {
1549  for( int k = 0; k < ndim; ++k )
1550  {
1551  bs[k * npts2fit + i] = 0;
1552  }
1553  }
1554  }
1555  ncols_sub = degree + 1 - interp;
1556  }
1557 
1558  // step 7. compute Q'*bs
1559  DGMSolver::compute_qtransposeB( npts2fit, ncols_sub, &( V[0] ), ndim, bs );
1560 
1561  // step 8. perform backward substitution and scale solutions
1562  // assign diagonals of V
1563  for( int i = 0; i < ncols_sub; ++i )
1564  {
1565  V[i * npts2fit + i] = D[i];
1566  }
1567  // backsolve
1568  if( safeguard )
1569  {
1570  DGMSolver::backsolve_polyfit_safeguarded( 1, degree, interp, npts2fit, ncols, &( V[0] ), ndim, bs, ws,
1571  degree_out );
1572  }
1573  else
1574  {
1575  DGMSolver::backsolve( npts2fit, ncols_sub, &( V[0] ), ndim, bs, &( ts[0] ) );
1576  *degree_out = degree;
1577  }
1578 }

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
vidEntityHandle, a vertex in _verts2rec
geomtypeGEOMTYPE, one of HISURFACE,HI3DCURVE,HI2DCURVE
coordsvector, global coordinates of local uvw coordinate system axis directions will be appended to the end of coords
degree_outReference to Integer, order of polynomial fittings for vid
coeffsvector, coefficients of local polynomial fittings in monomial basis will be appended to the end of coeffs
interpReference to Boolean, true = interpolation

Definition at line 627 of file HiReconstruction.cpp.

633 {
634  if( !_hasfittings )
635  {
636  return false;
637  }
638  else
639  {
640  int index = _verts2rec.index( vid );
641  if( -1 == index )
642  {
643  std::cout << "Input vertex is not locally hosted vertex in this mesh set" << std::endl;
644  return false;
645  }
646  geomtype = _geom;
647  if( HISURFACE == _geom )
648  {
649  coords.insert( coords.end(), _local_coords.begin() + 9 * index, _local_coords.begin() + 9 * index + 9 );
650  degree_out = _degrees_out[index];
651  interp = _interps[index];
652  int ncoeffs = ( degree_out + 2 ) * ( degree_out + 1 ) >> 1;
653  size_t istr = _vertID2coeffID[index];
654  coeffs.insert( coeffs.end(), _local_fit_coeffs.begin() + istr, _local_fit_coeffs.begin() + istr + ncoeffs );
655  }
656  else if( HI3DCURVE == _geom )
657  {
658  coords.insert( coords.end(), _local_coords.begin() + 3 * index, _local_coords.begin() + 3 * index + 3 );
659  degree_out = _degrees_out[index];
660  interp = _interps[index];
661  int ncoeffs = 3 * ( degree_out + 1 );
662  size_t istr = _vertID2coeffID[index];
663  coeffs.insert( coeffs.end(), _local_fit_coeffs.begin() + istr, _local_fit_coeffs.begin() + istr + ncoeffs );
664  }
665  return true;
666  }
667 }

References _degrees_out, _geom, _hasfittings, _interps, _local_coords, _local_fit_coeffs, _vertID2coeffID, _verts2rec, moab::HI3DCURVE, moab::HISURFACE, moab::index, 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
vertshRange, EntityHandles of vertices
nrmsPointer of array of doubles, size = 3*vertsh.size()

Definition at line 933 of file HiReconstruction.cpp.

934 {
936  if( _hasderiv )
937  {
938  size_t id = 0;
939  for( Range::iterator ivert = vertsh.begin(); ivert != vertsh.end(); ++ivert, ++id )
940  {
941  int index = _verts2rec.index( *ivert );
942 #ifdef MOAB_HAVE_MPI
943  if( -1 == index )
944  {
945  // ghost vertex
946  error = average_vertex_normal( *ivert, nrms + 3 * id );
947  MB_CHK_ERR( error );
948  }
949  else
950  {
951  nrms[3 * id] = _local_coords[9 * index + 6];
952  nrms[3 * id + 1] = _local_coords[9 * index + 7];
953  nrms[3 * id + 2] = _local_coords[9 * index + 8];
954  }
955 #else
956  assert( -1 != index );
957  nrms[3 * id] = _local_coords[9 * index + 6];
958  nrms[3 * id + 1] = _local_coords[9 * index + 7];
959  nrms[3 * id + 2] = _local_coords[9 * index + 8];
960 #endif
961  }
962  }
963  else
964  {
965  size_t id = 0;
966  for( Range::iterator ivert = vertsh.begin(); ivert != vertsh.end(); ++ivert, ++id )
967  {
968  error = average_vertex_normal( *ivert, nrms + 3 * id );
969  MB_CHK_ERR( error );
970  }
971  }
972  return error;
973 }

References _hasderiv, _local_coords, _verts2rec, average_vertex_normal(), moab::Range::begin(), moab::Range::end(), moab::error(), ErrorCode, moab::index, 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
vertshRange, EntityHandles of vertices
tangsPointer of array of doubles, size = 3*vertsh.size()

Definition at line 1024 of file HiReconstruction.cpp.

1025 {
1027  if( _hasderiv )
1028  {
1029  size_t id = 0;
1030  for( Range::iterator ivert = vertsh.begin(); ivert != vertsh.end(); ++ivert, ++id )
1031  {
1032  int index = _verts2rec.index( *ivert );
1033 #ifdef MOAB_HAVE_MPI
1034  if( -1 != index )
1035  {
1036  tangs[3 * id] = _local_coords[3 * index];
1037  tangs[3 * id + 1] = _local_coords[3 * index + 1];
1038  tangs[3 * id + 2] = _local_coords[3 * index + 2];
1039  }
1040  else
1041  {
1042  error = average_vertex_tangent( *ivert, tangs + 3 * id );
1043  MB_CHK_ERR( error );
1044  }
1045 #else
1046  assert( -1 != index );
1047  tangs[3 * id] = _local_coords[3 * index];
1048  tangs[3 * id + 1] = _local_coords[3 * index + 1];
1049  tangs[3 * id + 2] = _local_coords[3 * index + 2];
1050 #endif
1051  }
1052  }
1053  else
1054  {
1055  size_t id = 0;
1056  for( Range::iterator ivert = vertsh.begin(); ivert != vertsh.end(); ++ivert, ++id )
1057  {
1058  error = average_vertex_tangent( *ivert, tangs + 3 * id );
1059  MB_CHK_ERR( error );
1060  }
1061  }
1062  return error;
1063 }

References _hasderiv, _local_coords, _verts2rec, average_vertex_tangent(), moab::Range::begin(), moab::Range::end(), moab::error(), ErrorCode, moab::index, 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
vidEntityHandle, the vertex around which to perform high order projection.
npts2fitInteger, number of points lying around vid to be fitted.
coords2fitPointer to array of doubles, size=3*npts2fit, current coordinates of points to be projected.
newcoordsPointer to array of doubles, preallocated by user, size=3*npts2fit, estimated positions of input points.

Definition at line 487 of file HiReconstruction.cpp.

491 {
492  if( !_hasfittings )
493  {
494  MB_SET_ERR( MB_FAILURE, "There is no existing fitting results" );
495  }
496  else if( -1 == _verts2rec.index( vid ) )
497  {
498  std::ostringstream convert;
499  convert << vid;
500  std::string VID = convert.str();
501  MB_SET_ERR( MB_FAILURE, "There is no existing fitting results for vertex " + VID );
502  }
504  // get center of local coordinates system
505  double local_origin[3];
506  error = mbImpl->get_coords( &vid, 1, local_origin );
507  MB_CHK_ERR( error );
508  // get local fitting parameters
509  int index = _verts2rec.index( vid );
510  bool interp = _interps[index];
511  int local_deg = _degrees_out[index];
512  double *uvw_coords, *local_coeffs;
513  if( _geom == HISURFACE )
514  {
515  uvw_coords = &( _local_coords[9 * index] );
516  // int ncoeffs = (local_deg+2)*(local_deg+1)>>1;
517  size_t istr = _vertID2coeffID[index];
518  local_coeffs = &( _local_fit_coeffs[istr] );
519  walf3d_surf_vertex_eval( local_origin, uvw_coords, local_deg, local_coeffs, interp, npts2fit, coords2fit,
520  hiproj_new );
521  }
522  else if( _geom == HI3DCURVE )
523  {
524  uvw_coords = &( _local_coords[3 * index] );
525  size_t istr = _vertID2coeffID[index];
526  local_coeffs = &( _local_fit_coeffs[istr] );
527  walf3d_curve_vertex_eval( local_origin, uvw_coords, local_deg, local_coeffs, interp, npts2fit, coords2fit,
528  hiproj_new );
529  }
530  return error;
531 }

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::index, 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
elemEntityHandle, the element on which to perform high order projection.
nvpeInteger, number of nodes of this element, triangle is 3, quad is four.
npts2fitInteger, number of points lying in elem to be projected.
naturalcoords2fitPointer to array of doubles, size=nvpe*npts2fit, natural coordinates in elem of points to be projected.
newcoordsPointer to array of doubles, preallocated by user, size=3*npts2fit, estimated positions of input points.

Definition at line 405 of file HiReconstruction.cpp.

410 {
411  assert( newcoords );
413  // get connectivity table
414  std::vector< EntityHandle > elemconn;
415  error = mbImpl->get_connectivity( &elem, 1, elemconn );
416  MB_CHK_ERR( error );
417  if( nvpe != (int)elemconn.size() )
418  {
419  MB_SET_ERR( MB_FAILURE, "element connectivity table size doesn't match input size" );
420  }
421 
422  if( !_hasfittings )
423  {
424  MB_SET_ERR( MB_FAILURE, "There is no existing fitting results" );
425  }
426  else
427  {
428  std::ostringstream convert;
429  convert << elem;
430  std::string ID = convert.str();
431  for( int i = 0; i < nvpe; ++i )
432  {
433  if( -1 == _verts2rec.index( elemconn[i] ) )
434  {
435  MB_SET_ERR( MB_FAILURE, "There is no existing fitting results for element " + ID );
436  }
437  }
438  }
439  // check correctness of input
440  for( int i = 0; i < npts2fit; ++i )
441  {
442  if( !check_barycentric_coords( nvpe, naturalcoords2fit + i * nvpe ) )
443  {
444  MB_SET_ERR( MB_FAILURE, "Wrong barycentric coordinates" );
445  }
446  }
447 
448  double* elemcoords = new double[nvpe * 3];
449  error = mbImpl->get_coords( &( elemconn[0] ), nvpe, elemcoords );
450  MB_CHK_ERR( error );
451 
452  double* coords2fit = new double[3 * npts2fit]();
453  for( int i = 0; i < npts2fit; ++i )
454  {
455  for( int j = 0; j < nvpe; ++j )
456  {
457  coords2fit[3 * i] += naturalcoords2fit[i * nvpe + j] * elemcoords[3 * j];
458  coords2fit[3 * i + 1] += naturalcoords2fit[i * nvpe + j] * elemcoords[3 * j + 1];
459  coords2fit[3 * i + 2] += naturalcoords2fit[i * nvpe + j] * elemcoords[3 * j + 2];
460  }
461  }
462 
463  double* hiproj_new = new double[3 * npts2fit];
464  // initialize output
465  for( int i = 0; i < npts2fit; ++i )
466  {
467  newcoords[3 * i] = newcoords[3 * i + 1] = newcoords[3 * i + 2] = 0;
468  }
469  // for each input vertex, call nvpe fittings and take average
470  for( int j = 0; j < nvpe; ++j )
471  {
472  error = hiproj_walf_around_vertex( elemconn[j], npts2fit, coords2fit, hiproj_new );
473  MB_CHK_ERR( error );
474  for( int i = 0; i < npts2fit; ++i )
475  {
476  newcoords[3 * i] += naturalcoords2fit[i * nvpe + j] * hiproj_new[3 * i];
477  newcoords[3 * i + 1] += naturalcoords2fit[i * nvpe + j] * hiproj_new[3 * i + 1];
478  newcoords[3 * i + 2] += naturalcoords2fit[i * nvpe + j] * hiproj_new[3 * i + 2];
479  }
480  }
481  delete[] elemcoords;
482  delete[] coords2fit;
483  delete[] hiproj_new;
484  return error;
485 }

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.

52 {
54 
55 #ifdef HIREC_USE_AHF
56  std::cout << "HIREC_USE_AHF: Initializing" << std::endl;
57  ahf = new HalfFacetRep( mbImpl, pcomm, _mesh2rec, false );
58  if( !ahf )
59  {
61  }
62  error = ahf->initialize();
63  MB_CHK_ERR( error );
64 #else
65  ahf = NULL;
66 #endif
67 
68  // error = ahf->get_entity_ranges(_inverts,_inedges,_infaces,_incells); MB_CHK_ERR(error);
70  MB_CHK_ERR( error );
72  MB_CHK_ERR( error );
74  MB_CHK_ERR( error );
76  MB_CHK_ERR( error );
77  if( _inedges.size() && _infaces.empty() && _incells.empty() )
78  {
79  _dim = 1;
80  _MAXPNTS = 13;
81  }
82  else if( _infaces.size() && _incells.empty() )
83  {
84  _dim = 2;
85  _MAXPNTS = 128;
86  }
87  else
88  {
89  MB_SET_ERR( MB_FAILURE, "Encountered a non-manifold mesh or a mesh with volume elements" );
90  }
91 
92  // get locally hosted vertices by filtering pstatus
93 #ifdef MOAB_HAVE_MPI
94  if( pcomm )
95  {
97  MB_CHK_ERR( error );
98  }
99  else
100  {
102  }
103 #else
105 #endif
106  _nv2rec = _verts2rec.size();
107 
108  if( recwhole )
109  {
110  // compute normals(surface) or tangent vector(curve) for all locally hosted vertices
111  if( 2 == _dim )
112  {
114  }
115  else if( 1 == _dim )
116  {
118  }
119  else
120  {
121  MB_SET_ERR( MB_FAILURE, "Unknow space dimension" );
122  }
123  _hasderiv = true;
124  }
125  return error;
126 }

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 804 of file HiReconstruction.cpp.

805 {
806  if( !_hasderiv )
807  {
809  _hasderiv = true;
810  }
811  if( !_initfittings )
812  {
813  int ncoeffspvpd = degree + 1;
814  _degrees_out.assign( _nv2rec, 0 );
815  _interps.assign( _nv2rec, false );
816  _vertID2coeffID.resize( _nv2rec );
817  _local_fit_coeffs.assign( _nv2rec * ncoeffspvpd * 3, 0 );
818  for( size_t i = 0; i < _nv2rec; ++i )
819  {
820  _vertID2coeffID[i] = i * ncoeffspvpd * 3;
821  }
822  _initfittings = true;
823  }
824 }

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 826 of file HiReconstruction.cpp.

827 {
828  if( !_hasderiv )
829  {
831  _hasderiv = true;
832  }
833  if( !_hasfittings )
834  {
835  assert( _nv2rec == npts );
836  _degrees_out.assign( _nv2rec, 0 );
837  _interps.assign( _nv2rec, false );
838  _vertID2coeffID.reserve( _nv2rec );
839  size_t index = 0;
840  for( size_t i = 0; i < _nv2rec; ++i )
841  {
842  _vertID2coeffID[i] = index;
843  index += 3 * ( degrees[i] + 1 );
844  }
845  _local_fit_coeffs.assign( index, 0 );
846  _initfittings = true;
847  }
848 }

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

◆ 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 758 of file HiReconstruction.cpp.

759 {
760  if( !_hasderiv )
761  {
763  _hasderiv = true;
764  }
765  if( !_initfittings )
766  {
767  int ncoeffspv = ( degree + 2 ) * ( degree + 1 ) / 2;
768  _degrees_out.assign( _nv2rec, 0 );
769  _interps.assign( _nv2rec, false );
770  _vertID2coeffID.resize( _nv2rec );
771  _local_fit_coeffs.assign( _nv2rec * ncoeffspv, 0 );
772  for( size_t i = 0; i < _nv2rec; ++i )
773  {
774  _vertID2coeffID[i] = i * ncoeffspv;
775  }
776  _initfittings = true;
777  }
778 }

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 780 of file HiReconstruction.cpp.

781 {
782  if( !_hasderiv )
783  {
785  _hasderiv = true;
786  }
787  if( !_initfittings )
788  {
789  assert( _nv2rec == npts );
790  _degrees_out.assign( _nv2rec, 0 );
791  _interps.assign( _nv2rec, false );
792  _vertID2coeffID.resize( _nv2rec );
793  size_t index = 0;
794  for( size_t i = 0; i < _nv2rec; ++i )
795  {
796  _vertID2coeffID[i] = index;
797  index += ( degrees[i] + 2 ) * ( degrees[i] + 1 ) / 2;
798  }
799  _local_fit_coeffs.assign( index, 0 );
800  _initfittings = true;
801  }
802 }

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

◆ 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
vidEntityHandle, vertex around which to get n-ring vertices
ringInteger, number of rings
minpntsInteger, 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
ngbvsRange, the n-ring vertices of vid, including vid. If too many points found, i.e. more than _MAXPNTS, then terminate early.

Definition at line 694 of file HiReconstruction.cpp.

695 {
697  std::deque< EntityHandle > todo;
698  todo.push_back( vid );
699  ngbvs.insert( vid );
700  EntityHandle pre, nxt;
701  for( int i = 1; i <= ring; ++i )
702  {
703  int count = todo.size();
704  while( count )
705  {
706  EntityHandle center = todo.front();
707  todo.pop_front();
708  --count;
709  std::vector< EntityHandle > adjents;
711  MB_CHK_ERR( error );
712  for( size_t j = 0; j < adjents.size(); ++j )
713  {
714  std::vector< EntityHandle > elemconn;
715  error = mbImpl->get_connectivity( &adjents[j], 1, elemconn );
716  MB_CHK_ERR( error );
717  int nvpe = elemconn.size();
718  for( int k = 0; k < nvpe; ++k )
719  {
720  if( elemconn[k] == center )
721  {
722  pre = k == 0 ? elemconn[nvpe - 1] : elemconn[k - 1];
723  nxt = elemconn[( k + 1 ) % nvpe];
724  if( ngbvs.find( pre ) == ngbvs.end() )
725  {
726  ngbvs.insert( pre );
727  todo.push_back( pre );
728  }
729  if( ngbvs.find( nxt ) == ngbvs.end() )
730  {
731  ngbvs.insert( nxt );
732  todo.push_back( nxt );
733  }
734  break;
735  }
736  }
737  }
738  }
739  if( _MAXPNTS <= (int)ngbvs.size() )
740  {
741  // obtain enough points
742  return error;
743  }
744  if( !todo.size() )
745  {
746  // current ring cannot introduce any points, return incase deadlock
747  return error;
748  }
749  if( ( i == ring ) && ( minpnts > (int)ngbvs.size() ) )
750  {
751  // reach maximum ring but not enough points
752  ++ring;
753  }
754  }
755  return error;
756 }

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
nvertsInteger, number of vertices of input
ngbcorsPointer to array of doubles, size = 3*nverts, coordinates of input vertices, first will be center
ngbtangsPointer to array of doubles, size = 3*nverts, vertex tangent vectors of input vertices
degreeInteger, user specified fitting degree
interpBoolean, user input, interpolation or not
safeguardBoolean, true = use safeguarded numerical method in computing
ncoordsInteger, size of *coords, should be 3, used for check
coordsPointer to array of doubles, preallocated memory for storing the glocal coordinates of local u axis direction
ncoeffsInteger, size of coeffs, should be 3(degree+1), used for check
coeffsPointer to array of doubles, preallocated memory for storing coefficients of local fittings in monomial basis
degree_outPointer to integer, order of resulted polynomial of fitting, could be downgraded due to numerical issues

Definition at line 1354 of file HiReconstruction.cpp.

1365 {
1366  if( !nverts )
1367  {
1368  return;
1369  }
1370  // step 1. compute local coordinates system
1371  double tang[3] = { ngbtangs[0], ngbtangs[1], ngbtangs[2] };
1372  if( coords && ncoords > 2 )
1373  {
1374  coords[0] = tang[0];
1375  coords[1] = tang[1];
1376  coords[2] = tang[2];
1377  }
1378  if( !coeffs || !ncoeffs )
1379  {
1380  return;
1381  }
1382  else
1383  {
1384  assert( ncoeffs >= 3 * ( degree + 1 ) );
1385  }
1386  // step 2. project onto local coordinate system
1387  int npts2fit = nverts - interp;
1388  if( !npts2fit )
1389  {
1390  *degree_out = 0;
1391  for( int i = 0; i < ncoeffs; ++i )
1392  {
1393  coeffs[0] = 0;
1394  }
1395  return;
1396  }
1397  std::vector< double > us( npts2fit ); // double *us = new double[npts2fit];
1398  std::vector< double > bs( npts2fit * 3 ); // double *bs = new double[npts2fit*3];
1399  double uu[3];
1400  for( int i = interp; i < nverts; ++i )
1401  {
1402  int k = i - interp;
1403  DGMSolver::vec_linear_operation( 3, 1, ngbcors + 3 * i, -1, ngbcors, uu );
1404  us[k] = DGMSolver::vec_innerprod( 3, uu, tang );
1405  bs[k] = uu[0];
1406  bs[npts2fit + k] = uu[1];
1407  bs[2 * npts2fit + k] = uu[2];
1408  }
1409 
1410  // step 3. copmute weights
1411  std::vector< double > ws( npts2fit );
1412  int nzeros = compute_weights( npts2fit, 1, &( us[0] ), nverts, ngbtangs, degree, _MINEPS, &( ws[0] ) );
1413  assert( nzeros <= npts2fit );
1414 
1415  // step 4. adjust according to number of zero-weights
1416  if( nzeros )
1417  {
1418  if( nzeros == npts2fit )
1419  {
1420  // singular case
1421  *degree_out = 0;
1422  for( int i = 0; i < ncoeffs; ++i )
1423  {
1424  coeffs[i] = 0;
1425  }
1426  return;
1427  }
1428  int npts_new = npts2fit - nzeros;
1429  std::vector< double > bs_temp( 3 * npts_new );
1430  int index = 0;
1431  for( int i = 0; i < npts2fit; ++i )
1432  {
1433  if( ws[i] )
1434  {
1435  if( i > index )
1436  {
1437  us[index] = us[i];
1438  ws[index] = ws[i];
1439  }
1440  bs_temp[index] = bs[i];
1441  bs_temp[index + npts_new] = bs[i + npts2fit];
1442  bs_temp[index + 2 * npts_new] = bs[i + 2 * npts2fit];
1443  ++index;
1444  }
1445  }
1446  assert( index == npts_new );
1447  npts2fit = npts_new;
1448  us.resize( index );
1449  ws.resize( index );
1450  bs = bs_temp;
1451  // destroy bs_temp;
1452  std::vector< double >().swap( bs_temp );
1453  }
1454 
1455  // step 5. fitting
1456  eval_vander_univar_cmf( npts2fit, &( us[0] ), 3, &( bs[0] ), degree, &( ws[0] ), interp, safeguard, degree_out );
1457  // step 6. write results to output pointers
1458  int ncoeffs_out_pvpd = *degree_out + 1;
1459  assert( ncoeffs >= 3 * ncoeffs_out_pvpd );
1460  // write to coeffs consecutively, bs's size is not changed by eval_vander_univar_cmf
1461  coeffs[0] = coeffs[ncoeffs_out_pvpd] = coeffs[2 * ncoeffs_out_pvpd] = 0;
1462  for( int i = 0; i < ncoeffs_out_pvpd - interp; ++i )
1463  {
1464  coeffs[i + interp] = bs[i];
1465  coeffs[i + interp + ncoeffs_out_pvpd] = bs[i + npts2fit];
1466  coeffs[i + interp + 2 * ncoeffs_out_pvpd] = bs[i + 2 * npts2fit];
1467  }
1468 }

References _MINEPS, compute_weights(), eval_vander_univar_cmf(), moab::index, 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
nvertsInteger, number of vertices of input
ngbcorsPointer to array of doubles, size = 3*nverts, coordinates of input vertices, first will be center
ngbnrmsPointer to array of doubles, size = 3*nverts, vertex normals of input vertices
degreeInteger, user specified fitting degree
interpBoolean, user input, interpolation or not
safeguardBoolean, true = use safeguarded numerical method in computing
ncoordsInteger, size of *coords, should be 9, used for check
coordsPointer to array of doubles, preallocated memory for storing the glocal coordinates of local uvw axis directions
ncoeffsInteger, size of *coeffs, should be (degree+2)(degree+1)/2, used for check
coeffsPointer to array of doubles, preallocated memory for storing coefficients of local fittings in monomial basis
degree_outPointer to integer, order of resulted polynomial of fitting, could be downgraded due to numerical issues
degree_pntPointer to integer, polynomial fitting order determined by stencil size/number of points
degree_qrPointer to integer, polynomial fitting order determined by Vandermonde system condition number

Definition at line 1069 of file HiReconstruction.cpp.

1082 {
1083  if( nverts <= 0 )
1084  {
1085  return;
1086  }
1087 
1088  // std::cout << "npnts in initial stencil = " << nverts << std::endl;
1089  // std::cout << "centered at (" << ngbcoords[0] << "," << ngbcoords[1] << "," << ngbcoords[2] <<
1090  // ")" << std::endl;
1091 
1092  // step 1. copmute local coordinate system
1093  double nrm[3] = { ngbnrms[0], ngbnrms[1], ngbnrms[2] }, tang1[3] = { 0, 0, 0 }, tang2[3] = { 0, 0, 0 };
1094  if( fabs( nrm[0] ) > fabs( nrm[1] ) && fabs( nrm[0] ) > fabs( nrm[2] ) )
1095  {
1096  tang1[1] = 1.0;
1097  }
1098  else
1099  {
1100  tang1[0] = 1.0;
1101  }
1102 
1103  DGMSolver::vec_projoff( 3, tang1, nrm, tang1 );
1104 #ifndef NDEBUG
1105  assert( DGMSolver::vec_normalize( 3, tang1, tang1 ) );
1106 #endif
1107  DGMSolver::vec_crossprod( nrm, tang1, tang2 );
1108  if( 9 <= ncoords && coords )
1109  {
1110  coords[0] = tang1[0];
1111  coords[1] = tang1[1];
1112  coords[2] = tang1[2];
1113  coords[3] = tang2[0];
1114  coords[4] = tang2[1];
1115  coords[5] = tang2[2];
1116  coords[6] = nrm[0];
1117  coords[7] = nrm[1];
1118  coords[8] = nrm[2];
1119  }
1120  if( !ncoeffs || !coeffs )
1121  {
1122  return;
1123  }
1124  else
1125  {
1126  assert( ncoeffs >= ( degree + 2 ) * ( degree + 1 ) / 2 );
1127  }
1128 
1129  // step 2. project onto local coordinates system
1130  int npts2fit = nverts - interp;
1131  if( 0 == npts2fit )
1132  {
1133  *degree_out = *degree_pnt = *degree_qr = 0;
1134  for( int i = 0; i < ncoeffs; ++i )
1135  {
1136  coeffs[i] = 0;
1137  }
1138  // coeffs[0] = 0;
1139  return;
1140  }
1141  std::vector< double > us( npts2fit * 2 ); // double *us = new double[npts2fit*2];
1142  std::vector< double > bs( npts2fit ); // double *bs = new double[npts2fit];
1143  for( int i = interp; i < nverts; ++i )
1144  {
1145  int k = i - interp;
1146  double uu[3];
1147  DGMSolver::vec_linear_operation( 3, 1, ngbcoords + 3 * i, -1, ngbcoords, uu );
1148  us[k * 2] = DGMSolver::vec_innerprod( 3, tang1, uu );
1149  us[k * 2 + 1] = DGMSolver::vec_innerprod( 3, tang2, uu );
1150  bs[k] = DGMSolver::vec_innerprod( 3, nrm, uu );
1151  }
1152 
1153  // step 3. compute weights
1154  std::vector< double > ws( npts2fit ); // double *ws = new double[npts2fit];
1155  int nzeros = compute_weights( npts2fit, 2, &( us[0] ), nverts, ngbnrms, degree, _MINEPS, &( ws[0] ) );
1156 
1157  // step 4. adjust according to zero-weights
1158  if( nzeros )
1159  {
1160  if( nzeros == npts2fit )
1161  {
1162  *degree_out = *degree_pnt = *degree_qr = 0;
1163  for( int i = 0; i < ncoeffs; ++i )
1164  {
1165  coeffs[i] = 0;
1166  }
1167  // coeffs[0] = 0;
1168  return;
1169  }
1170  int index = 0;
1171  for( int i = 0; i < npts2fit; ++i )
1172  {
1173  if( ws[i] )
1174  {
1175  if( i > index )
1176  {
1177  us[index * 2] = us[i * 2];
1178  us[index * 2 + 1] = us[i * 2 + 1];
1179  bs[index] = bs[i];
1180  ws[index] = ws[i];
1181  }
1182  ++index;
1183  }
1184  }
1185  npts2fit -= nzeros;
1186  assert( index == npts2fit );
1187  us.resize( npts2fit * 2 );
1188  bs.resize( npts2fit );
1189  ws.resize( npts2fit );
1190  /*us = realloc(us,npts2fit*2*sizeof(double));
1191  bs = realloc(bs,npts2fit*sizeof(double));
1192  ws = realloc(ws,npts2fit*sizeof(double));*/
1193  }
1194 
1195  // std::cout<<"npnts after weighting = "<<npts2fit<<std::endl;
1196 
1197  // step 5. fitting
1198  eval_vander_bivar_cmf( npts2fit, &( us[0] ), 1, &( bs[0] ), degree, &( ws[0] ), interp, safeguard, degree_out,
1199  degree_pnt, degree_qr );
1200 
1201  // step 6. organize output
1202  int ncoeffs_out = ( *degree_out + 2 ) * ( *degree_out + 1 ) / 2;
1203  assert( ncoeffs_out <= ncoeffs );
1204  coeffs[0] = 0;
1205  for( int j = 0; j < ncoeffs_out - interp; ++j )
1206  {
1207  coeffs[j + interp] = bs[j];
1208  }
1209  // delete [] us; delete [] bs; delete [] ws;
1210 }

References _MINEPS, compute_weights(), eval_vander_bivar_cmf(), moab::index, 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
vidEntityHandle, the fittings will be performed around this vertex.
interpBoolean, true=Interpolation, false=least square fitting.
degreeInteger, order of polynomials used for local fittings.
minpntsInteger, 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.
safeguardBoolean, true=using safe guarded method in numerical computing.
coordsPointer to double, preallocated memory by user, should have at least 3 doubles; stores the global coordinates of local coordinate system u direction.
degree_outPointer to integer, used to store the degree of resulted fitting
coeffs,Pointerto double, preallocated memory for coefficients of local fittings, should have at least 3*(degree+1) doubles.

Definition at line 357 of file HiReconstruction.cpp.

367 {
369  int ring = estimate_num_rings( degree, interp );
370  // get n-ring neighbors
371  Range ngbvs;
372  error = obtain_nring_ngbvs( vid, ring, minpnts, ngbvs );
373  MB_CHK_ERR( error );
374  // get coordinates
375  size_t nverts = ngbvs.size();
376  assert( nverts );
377  double* ngbcoords = new double[nverts * 3];
378  error = mbImpl->get_coords( ngbvs, ngbcoords );
379  MB_CHK_ERR( error );
380  // get tangent vectors
381  double* ngbtangs = new double[nverts * 3];
382  error = get_tangents_curve( ngbvs, ngbtangs );
383  MB_CHK_ERR( error );
384  // switch vid to first one
385  int index = ngbvs.index( vid );
386  assert( index != -1 );
387  std::swap( ngbcoords[0], ngbcoords[3 * index] );
388  std::swap( ngbcoords[1], ngbcoords[3 * index + 1] );
389  std::swap( ngbcoords[2], ngbcoords[3 * index + 2] );
390  std::swap( ngbtangs[0], ngbtangs[3 * index] );
391  std::swap( ngbtangs[1], ngbtangs[3 * index + 1] );
392  std::swap( ngbtangs[2], ngbtangs[3 * index + 2] );
393  // local WLS fittings
394  polyfit3d_curve_get_coeff( nverts, ngbcoords, ngbtangs, degree, interp, safeguard, ncoords, coords, ncoeffs, coeffs,
395  degree_out );
396  delete[] ngbcoords;
397  delete[] ngbtangs;
398  return error;
399 }

References moab::error(), ErrorCode, estimate_num_rings(), moab::Core::get_coords(), get_tangents_curve(), moab::index, 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
vidEntityHandle, the fitting will be performed around this vertex for the local height function over the uv-plane.
interpBoolean, true=Interpolation, false=least square fitting.
degreeInteger, order of polynomials used for local fittings.
minpntsInteger, 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.
safeguardBoolean, true=using safe guarded method in numerical computing.
coordsPointer to double, preallocated memory by user, should have at least 9 doubles; stores the global coordinates of local coordinates system uvw directions.
degree_outPointer to integer, used to store the degree of resulted fitting
coeffs,Pointerto double, preallocated memory for coefficients of local fittings, should have at least (degree+2)(degree+1)/2 doubles.

Definition at line 304 of file HiReconstruction.cpp.

314 {
315  assert( _dim == 2 );
317  int ring = estimate_num_rings( degree, interp );
318  // std::cout<<"ring = "<<ring<<std::endl;
319  // get n-ring neighbors
320  Range ngbvs;
321  error = obtain_nring_ngbvs( vid, ring, minpnts, ngbvs );
322  MB_CHK_ERR( error );
323  // for debug
324  /*if(_verts2rec.index(vid)==70){
325  for(Range::iterator ingb=ngbvs.begin();ingb!=ngbvs.end();++ingb) std::cerr <<
326  _verts2rec.index(*ingb) << " "; std::cout << std::endl;
327  }*/
328 
329  // get coordinates;
330  size_t nverts = ngbvs.size();
331  assert( nverts );
332  double* ngbcoords = new double[nverts * 3];
333  error = mbImpl->get_coords( ngbvs, ngbcoords );
334  MB_CHK_ERR( error );
335  // get normals
336  double* ngbnrms = new double[nverts * 3];
337  error = get_normals_surf( ngbvs, ngbnrms );
338  MB_CHK_ERR( error );
339  // switch vid to first one
340  int index = ngbvs.index( vid );
341  assert( index != -1 );
342  std::swap( ngbcoords[0], ngbcoords[3 * index] );
343  std::swap( ngbcoords[1], ngbcoords[3 * index + 1] );
344  std::swap( ngbcoords[2], ngbcoords[3 * index + 2] );
345  std::swap( ngbnrms[0], ngbnrms[3 * index] );
346  std::swap( ngbnrms[1], ngbnrms[3 * index + 1] );
347  std::swap( ngbnrms[2], ngbnrms[3 * index + 2] );
348  // local WLS fitting
349  int degree_pnt, degree_qr;
350  polyfit3d_surf_get_coeff( nverts, ngbcoords, ngbnrms, degree, interp, safeguard, ncoords, coords, ncoeffs, coeffs,
351  degree_out, &degree_pnt, &degree_qr );
352  delete[] ngbcoords;
353  delete[] ngbnrms;
354  return error;
355 }

References _dim, moab::error(), ErrorCode, estimate_num_rings(), moab::Core::get_coords(), get_normals_surf(), moab::index, 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
degreeInteger, order of polynomials used for local fittings.
interpBoolean, true=Interpolation, false=least square fitting.
safeguardBoolean, true=using safe guarded method in numerical computing.
resetBoolean, reset=true will recompute the fittings based on user input and replace the existing one.

Definition at line 229 of file HiReconstruction.cpp.

230 {
231  assert( _dim == 1 );
232  if( _hasfittings && !reset )
233  {
234  return MB_SUCCESS;
235  }
236  else
237  {
238  _initfittings = _hasfittings = false;
239  }
240  initialize_3Dcurve_geom( degree );
242  double *coords = 0, *coeffs;
243  int* degree_out;
244  int ncoeffs = 3 * ( degree + 1 );
245  for( Range::iterator ivert = _verts2rec.begin(); ivert != _verts2rec.end(); ++ivert )
246  {
247  int index = _verts2rec.index( *ivert );
248  assert( index != -1 );
249  size_t istr = _vertID2coeffID[index];
250  coeffs = &( _local_fit_coeffs[istr] );
251  degree_out = &( _degrees_out[index] );
252  _interps[index] = interp;
253  error = polyfit3d_walf_curve_vertex( *ivert, interp, degree, _MINPNTS, safeguard, 0, coords, degree_out,
254  ncoeffs, coeffs );
255  MB_CHK_ERR( error );
256  }
257  _geom = HI3DCURVE;
258  _hasfittings = true;
259  return error;
260 }

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::index, 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
nptsInteger size of array pointed by degrees, used for check
degreesInteger arrray, order of polynomials for local fitting at all hosted vertices.
interpBoolean, true=Interpolation, false=least square fitting.
safeguardBoolean, true=using safe guarded method in numerical computing.
resetBoolean, reset=true will recompute the fittings based on user input and replace the existing one.

Definition at line 262 of file HiReconstruction.cpp.

267 {
268  assert( _dim == 1 );
270  if( npts != _nv2rec )
271  {
272  MB_SET_ERR( MB_FAILURE, "Input number of degrees doesn't match the number of vertices" );
273  }
274  if( _hasfittings && !reset )
275  {
276  return MB_SUCCESS;
277  }
278  else
279  {
280  _initfittings = _hasfittings = false;
281  }
282  // initialize
283  initialize_3Dcurve_geom( npts, degrees );
284  double *coords = 0, *coeffs;
285  int* degree_out;
286  size_t i = 0;
287  for( Range::iterator ivert = _verts2rec.begin(); ivert != _verts2rec.end(); ++ivert, ++i )
288  {
289  int index = _verts2rec.index( *ivert );
290  size_t istr = _vertID2coeffID[index];
291  coeffs = &( _local_fit_coeffs[istr] );
292  degree_out = &( _degrees_out[index] );
293  _interps[index] = interps[i];
294  int ncoeffs = 3 * ( degrees[i] + 1 );
295  error = polyfit3d_walf_curve_vertex( *ivert, interps[i], degrees[i], _MINPNTS, safeguard, 0, coords, degree_out,
296  ncoeffs, coeffs );
297  MB_CHK_ERR( error );
298  }
299  _geom = HI3DCURVE;
300  _hasfittings = true;
301  return error;
302 }

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::index, 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
degreeInteger, order of polynomials used for local fittings.
interpBoolean, true=Interpolation, false=least square fitting.
safeguardBoolean, true=using safe guarded method in numerical computing.
resetBoolean, reset=true will recompute the fittings based on user input and replace the existing one.

Definition at line 132 of file HiReconstruction.cpp.

133 {
134  assert( 2 == _dim );
135  if( _hasfittings && !reset )
136  {
137  // This object has precomputed fitting results and user don't want to reset
138  return MB_SUCCESS;
139  }
140  else
141  {
142  _initfittings = _hasfittings = false;
143  }
144  // initialize for geometric information
145  initialize_surf_geom( degree );
147  double *coeffs, *coords;
148  int* degree_out;
149  int ncoeffs = ( degree + 2 ) * ( degree + 1 ) / 2;
150 
151  // DBG
152  // int dcount = 0;
153 
154  for( Range::iterator ivert = _verts2rec.begin(); ivert != _verts2rec.end(); ++ivert )
155  {
156  int index = _verts2rec.index( *ivert );
157  // for debug
158  /*if(index==70){
159  EntityHandle vid = *ivert;
160  double vertcoords[3];
161  error = mbImpl->get_coords(&vid,1,vertcoords);
162  }*/
163 
164  size_t istr = _vertID2coeffID[index];
165  coords = &( _local_coords[9 * index] );
166  coeffs = &( _local_fit_coeffs[istr] );
167  degree_out = &( _degrees_out[index] );
168  _interps[index] = interp;
169  error = polyfit3d_walf_surf_vertex( *ivert, interp, degree, _MINPNTS, safeguard, 9, coords, degree_out, ncoeffs,
170  coeffs );
171  MB_CHK_ERR( error );
172 
173  // DBG
174  // if( degree_out[0] < degree ) dcount += 1;
175  }
176 
177  // DBG
178  // std::cout<<"Total #points ="<<_verts2rec.size()<<", #degraded points = "<<dcount<<std::endl;
179 
180  _geom = HISURFACE;
181  _hasfittings = true;
182  return error;
183 }

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::index, 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
nptsInteger size of array pointed by degrees, used for check
degreesInteger arrray, order of polynomials for local fitting at all hosted vertices
interpBoolean, true=Interpolation, false=least square fitting.
safeguardBoolean, true=using safe guarded method in numerical computing.
resetBoolean, reset=true will recompute the fittings based on user input and replace the existing one.

Definition at line 185 of file HiReconstruction.cpp.

190 {
191  assert( _dim == 2 );
192  if( npts != _nv2rec )
193  {
194  MB_SET_ERR( MB_FAILURE, "Input number of degrees doesn't match number of vertices" );
195  }
196  if( _hasfittings && !reset )
197  {
198  return MB_SUCCESS;
199  }
200  else
201  {
202  _initfittings = _hasfittings = false;
203  }
205  // initialization for fitting
206  initialize_surf_geom( npts, degrees );
207  double *coeffs, *coords;
208  int* degree_out;
209  size_t i = 0;
210  for( Range::iterator ivert = _verts2rec.begin(); ivert != _verts2rec.end(); ++ivert, ++i )
211  {
212  int index = _verts2rec.index( *ivert );
213  assert( -1 != index );
214  size_t istr = _vertID2coeffID[index];
215  coords = &( _local_coords[9 * index] );
216  coeffs = &( _local_fit_coeffs[istr] );
217  degree_out = &( _degrees_out[index] );
218  _interps[index] = interps[i];
219  int ncoeffs = ( degrees[i] + 2 ) * ( degrees[i] + 1 ) / 2;
220  error = polyfit3d_walf_surf_vertex( *ivert, interps[i], degrees[i], _MINPNTS, safeguard, 9, coords, degree_out,
221  ncoeffs, coeffs );
222  MB_CHK_ERR( error );
223  }
224  _geom = HISURFACE;
225  _hasfittings = true;
226  return error;
227 }

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::index, 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
vidEntityHandle of vertex
elemdimInteger, dimension of elements incidented in vid
adjentsvector<EntityHandle>, container which push incident elements in

Definition at line 678 of file HiReconstruction.cpp.

681 {
683  assert( elemdim == _dim );
684 #ifdef HIREC_USE_AHF
685  error = ahf->get_up_adjacencies( vid, elemdim, adjents );
686  MB_CHK_ERR( error );
687 #else
688  error = mbImpl->get_adjacencies( &vid, 1, elemdim, false, adjents );
689  MB_CHK_ERR( error );
690 #endif
691  return error;
692 }

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_originPointer to 3 doubles, coordinates of the center vertex
local_coordsPointer to 3 doubles, global coordinates of direction of local u coordinate axis at center vertex
local_degInteger, order of local polynomial fitting
local_coeffsPointer 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.
interpBoolean, true=local fitting is interpolation, false=local fitting is least square fitting
npts2fitInteger, number of points to be estimated, around the center vertices
coords2fitPointer to array of doubles, size=3*npts2fit, current coordinates of points to be estimated
hiproj_newPointer to array of doubles, size=3*npts2fit, memory preallocated by user to store the fitting/estimated positions of input points.

Definition at line 589 of file HiReconstruction.cpp.

597 {
598  assert( local_origin && local_coords && local_coeffs );
599  int ncoeffspvpd = local_deg + 1;
600  for( int i = 0; i < npts2fit; ++i )
601  {
602  // get the vector from center to current point, project to tangent line
603  double vec[3], ans[3] = { 0, 0, 0 };
604  DGMSolver::vec_linear_operation( 3, 1, coords2fit + 3 * i, -1, local_origin, vec );
605  double u = DGMSolver::vec_innerprod( 3, local_coords, vec );
606  // evaluate polynomials
607  if( !interp )
608  {
609  ans[0] = local_coeffs[0];
610  ans[1] = local_coeffs[ncoeffspvpd];
611  ans[2] = local_coeffs[2 * ncoeffspvpd];
612  }
613  double uk = 1; // degree_out and degree different, stored in columnwise contiguously
614  for( int j = 1; j < ncoeffspvpd; ++j )
615  {
616  uk *= u;
617  ans[0] += uk * local_coeffs[j];
618  ans[1] += uk * local_coeffs[j + ncoeffspvpd];
619  ans[2] += uk * local_coeffs[j + 2 * ncoeffspvpd];
620  }
621  hiproj_new[3 * i] = ans[0] + local_origin[0];
622  hiproj_new[3 * i + 1] = ans[1] + local_origin[1];
623  hiproj_new[3 * i + 2] = ans[2] + local_origin[2];
624  }
625 }

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_originPointer to 3 doubles, coordinates of the center vertex
local_coordsPointer to 9 doubles, global coordinates of directions of local uvw coordinates axis at center vertex
local_degInteger, order of local polynomial fitting
local_coeffsPointer to array of doubles, size=(local_deg+2)(local_deg+1)/2, coefficients of local polynomial fittings, in monomial basis
interpBoolean, true=local fitting is interpolation, false=local fitting is least square fitting
npts2fitInteger, number of points to be estimated, around the center vertices
coords2fitPointer to array of doubles, size=3*npts2fit, current coordinates of points to be estimated
hiproj_newPointer to array of doubles, size=3*npts2fit, memory preallocated by user to store the fitting/estimated positions of input points.

Definition at line 533 of file HiReconstruction.cpp.

541 {
542  double xaxis[3], yaxis[3], zaxis[3];
543  for( int i = 0; i < 3; ++i )
544  {
545  xaxis[i] = local_coords[i];
546  yaxis[i] = local_coords[3 + i];
547  zaxis[i] = local_coords[6 + i];
548  }
549  // double *basis = new double[(local_deg+2)*(local_deg+1)/2-1];
550  std::vector< double > basis( ( local_deg + 2 ) * ( local_deg + 1 ) / 2 - 1 );
551  for( int i = 0; i < npts2fit; ++i )
552  {
553  double local_pos[3];
554  for( int j = 0; j < 3; ++j )
555  {
556  local_pos[j] = coords2fit[3 * i + j] - local_origin[j];
557  }
558  double u, v, height = 0;
559  u = DGMSolver::vec_innerprod( 3, local_pos, xaxis );
560  v = DGMSolver::vec_innerprod( 3, local_pos, yaxis );
561  basis[0] = u;
562  basis[1] = v;
563  int l = 1;
564  for( int k = 2; k <= local_deg; ++k )
565  {
566  ++l;
567  basis[l] = u * basis[l - k];
568  for( int id = 0; id < k; ++id )
569  {
570  ++l;
571  basis[l] = basis[l - k - 1] * v;
572  }
573  }
574  if( !interp )
575  {
576  height = local_coeffs[0];
577  }
578  for( int p = 0; p <= l; ++p )
579  {
580  height += local_coeffs[p + 1] * basis[p];
581  }
582  hiproj_new[3 * i] = local_origin[0] + u * xaxis[0] + v * yaxis[0] + height * zaxis[0];
583  hiproj_new[3 * i + 1] = local_origin[1] + u * xaxis[1] + v * yaxis[1] + height * zaxis[1];
584  hiproj_new[3 * i + 2] = local_origin[2] + u * xaxis[2] + v * yaxis[2] + height * zaxis[2];
585  }
586  // delete [] basis;
587 }

References moab::DGMSolver::vec_innerprod().

Referenced by hiproj_walf_around_vertex().

Member Data Documentation

◆ _degrees_out

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

◆ _dim

◆ _geom

GEOMTYPE moab::HiReconstruction::_geom
protected

◆ _hasderiv

◆ _hasfittings

◆ _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

◆ _interps

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

◆ _inverts

Range moab::HiReconstruction::_inverts
protected

Definition at line 295 of file HiReconstruction.hpp.

Referenced by initialize().

◆ _local_coords

◆ _local_fit_coeffs

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

◆ _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

◆ _MINPNTS

int moab::HiReconstruction::_MINPNTS
protected

Definition at line 298 of file HiReconstruction.hpp.

Referenced by reconstruct3D_curve_geom(), and reconstruct3D_surf_geom().

◆ _nv2rec

◆ _vertID2coeffID

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

◆ _verts2rec

◆ ahf

HalfFacetRep* moab::HiReconstruction::ahf
protected

◆ mbImpl

◆ 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: