Mesh Oriented datABase  (version 5.5.1)
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 841 of file HiReconstruction.cpp.

842 {
844  std::vector< EntityHandle > adjfaces;
845  error = vertex_get_incident_elements( vid, 2, adjfaces );MB_CHK_ERR( error );
846  int npolys = adjfaces.size();
847  if( !npolys )
848  {
849  MB_SET_ERR( MB_FAILURE, "Vertex has no incident 2D entities" );
850  }
851  else
852  {
853  double v1[3], v2[3], v3[3], a[3], b[3], c[3];
854  nrm[0] = nrm[1] = nrm[2] = 0;
855  for( int i = 0; i < npolys; ++i )
856  {
857  // get incident "triangles"
858  std::vector< EntityHandle > elemconn;
859  error = mbImpl->get_connectivity( &adjfaces[i], 1, elemconn );MB_CHK_ERR( error );
860  EntityHandle pre, nxt;
861  int nvpe = elemconn.size();
862  for( int j = 0; j < nvpe; ++j )
863  {
864  if( vid == elemconn[j] )
865  {
866  pre = j == 0 ? elemconn[nvpe - 1] : elemconn[j - 1];
867  nxt = elemconn[( j + 1 ) % nvpe];
868  break;
869  }
870  }
871  // compute area weighted normals
872  error = mbImpl->get_coords( &pre, 1, a );MB_CHK_ERR( error );
873  error = mbImpl->get_coords( &vid, 1, b );MB_CHK_ERR( error );
874  error = mbImpl->get_coords( &nxt, 1, c );MB_CHK_ERR( error );
875  DGMSolver::vec_linear_operation( 3, 1, c, -1, b, v1 );
876  DGMSolver::vec_linear_operation( 3, 1, a, -1, b, v2 );
877  DGMSolver::vec_crossprod( v1, v2, v3 );
878  DGMSolver::vec_linear_operation( 3, 1, nrm, 1, v3, nrm );
879  }
880 #ifndef NDEBUG
881  assert( DGMSolver::vec_normalize( 3, nrm, nrm ) );
882 #endif
883  }
884  return error;
885 }

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

944 {
946  std::vector< EntityHandle > adjedges;
947  error = vertex_get_incident_elements( vid, 1, adjedges );MB_CHK_ERR( error );
948  int nedges = adjedges.size();
949  if( !nedges )
950  {
951  MB_SET_ERR( MB_FAILURE, "Vertex has no incident edges" );
952  }
953  else
954  {
955  assert( nedges <= 2 );
956  tang[0] = tang[1] = tang[2] = 0;
957  for( int i = 0; i < nedges; ++i )
958  {
959  std::vector< EntityHandle > edgeconn;
960  error = mbImpl->get_connectivity( &adjedges[i], 1, edgeconn );
961  double istr[3], iend[3], t[3];
962  error = mbImpl->get_coords( &( edgeconn[0] ), 1, istr );
963  error = mbImpl->get_coords( &( edgeconn[1] ), 1, iend );
964  DGMSolver::vec_linear_operation( 3, 1, iend, -1, istr, t );
965  DGMSolver::vec_linear_operation( 3, 1, tang, 1, t, tang );
966  }
967 #ifndef NDEBUG
968  assert( DGMSolver::vec_normalize( 3, tang, tang ) );
969 #endif
970  }
971  return error;
972 }

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

1594 {
1595  double sum = 0;
1596  for( int i = 0; i < nws; ++i )
1597  {
1598  if( naturalcoords[i] < -_MINEPS )
1599  {
1600  return false;
1601  }
1602  sum += naturalcoords[i];
1603  }
1604  if( fabs( 1 - sum ) > _MINEPS )
1605  {
1606  return false;
1607  }
1608  else
1609  {
1610  return true;
1611  }
1612 }

References _MINEPS, and moab::sum().

Referenced by hiproj_walf_in_element().

◆ compute_average_vertex_normals_surf()

ErrorCode moab::HiReconstruction::compute_average_vertex_normals_surf ( )
protected

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

Definition at line 887 of file HiReconstruction.cpp.

888 {
889  if( _hasderiv )
890  {
891  return MB_SUCCESS;
892  }
894  _local_coords.assign( 9 * _nv2rec, 0 );
895  size_t index = 0;
896  for( Range::iterator ivert = _verts2rec.begin(); ivert != _verts2rec.end(); ++ivert, ++index )
897  {
898  error = average_vertex_normal( *ivert, &( _local_coords[9 * index + 6] ) );MB_CHK_ERR( error );
899  }
900  return error;
901 }

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

Referenced by initialize(), and initialize_surf_geom().

◆ compute_average_vertex_tangents_curve()

ErrorCode moab::HiReconstruction::compute_average_vertex_tangents_curve ( )
protected

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

Definition at line 974 of file HiReconstruction.cpp.

975 {
976  if( _hasderiv )
977  {
978  return MB_SUCCESS;
979  }
981  _local_coords.assign( 3 * _nv2rec, 0 );
982  size_t index = 0;
983  for( Range::iterator ivert = _verts2rec.begin(); ivert != _verts2rec.end(); ++ivert, ++index )
984  {
985  error = average_vertex_tangent( *ivert, &( _local_coords[3 * index] ) );MB_CHK_ERR( error );
986  }
987  return error;
988 }

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

Referenced by initialize(), and initialize_3Dcurve_geom().

◆ compute_weights()

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

Compute weights for points selected in weighted least square fittigns.

Definition at line 1544 of file HiReconstruction.cpp.

1552 {
1553  assert( nrows <= _MAXPNTS && ws );
1554  bool interp = false;
1555  if( nngbs != nrows )
1556  {
1557  assert( nngbs == nrows + 1 );
1558  interp = true;
1559  }
1560  double epsilon = 1e-2;
1561 
1562  // First, compute squared distance from each input piont to the center
1563  for( int i = 0; i < nrows; ++i )
1564  {
1565  ws[i] = DGMSolver::vec_innerprod( ncols, us + i * ncols, us + i * ncols );
1566  }
1567 
1568  // Second, compute a small correction termt o guard against zero
1569  double h = 0;
1570  for( int i = 0; i < nrows; ++i )
1571  {
1572  h += ws[i];
1573  }
1574  h /= (double)nrows;
1575 
1576  // Finally, compute the weights for each vertex
1577  int nzeros = 0;
1578  for( int i = 0; i < nrows; ++i )
1579  {
1580  double costheta = DGMSolver::vec_innerprod( 3, ngbnrms, ngbnrms + 3 * ( i + interp ) );
1581  if( costheta > toler )
1582  {
1583  ws[i] = costheta * pow( ws[i] / h + epsilon, -1 * (double)degree / 2.0 );
1584  }
1585  else
1586  {
1587  ws[i] = 0;
1588  ++nzeros;
1589  }
1590  }
1591  return nzeros;
1592 }

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

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

Referenced by polyfit3d_walf_curve_vertex(), and polyfit3d_walf_surf_vertex().

◆ eval_vander_bivar_cmf()

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

Form and solve Vandermonde system of bi-variables.

Definition at line 1176 of file HiReconstruction.cpp.

1187 {
1188  // step 1. adjust the degree according to number of points to fit
1189  int ncols = ( ( ( degree + 2 ) * ( degree + 1 ) ) >> 1 ) - interp;
1190  while( 1 < degree && npts2fit < ncols )
1191  {
1192  --degree;
1193  ncols = ( ( ( degree + 2 ) * ( degree + 1 ) ) >> 1 ) - interp;
1194  }
1195  *degree_pnt = degree;
1196 
1197  // std::cout << "degree_pnt: " << *degree_pnt << std::endl;
1198 
1199  // step 2. construct Vandermonde matrix, stored in columnwise
1200  std::vector< double > V; // V(npts2fit*(ncols+interp)); //double *V_init = new double[npts2fit*(ncols+interp)];
1201  DGMSolver::gen_vander_multivar( npts2fit, 2, us, degree, V );
1202  // remove the first column of 1s if interpolation
1203  if( interp )
1204  {
1205  V.erase( V.begin(), V.begin() + npts2fit );
1206  }
1207  /*double* V;
1208  if(interp){
1209  V = new double[npts2fit*ncols];
1210  std::memcpy(V,V_init+npts2fit,ncols*npts2fit*sizeof(double));
1211  delete [] V_init; V_init = 0;
1212  }else{
1213  V = V_init;
1214  }*/
1215 
1216  // step 3. Scale rows to assign different weights to different points
1217  for( int i = 0; i < npts2fit; ++i )
1218  {
1219  for( int j = 0; j < ncols; ++j )
1220  {
1221  V[j * npts2fit + i] *= ws[i];
1222  }
1223  for( int k = 0; k < ndim; ++k )
1224  {
1225  bs[k * npts2fit + i] *= ws[i];
1226  }
1227  }
1228 
1229  // step 4. scale columns to reduce condition number
1230  std::vector< double > ts( ncols ); // double *ts = new double[ncols];
1231  DGMSolver::rescale_matrix( npts2fit, ncols, &( V[0] ), &( ts[0] ) );
1232 
1233  // step 5. Perform Householder QR factorization
1234  std::vector< double > D( ncols ); // double *D = new double[ncols];
1235  int rank;
1236  DGMSolver::qr_polyfit_safeguarded( npts2fit, ncols, &( V[0] ), &( D[0] ), &rank );
1237 
1238  // step 6. adjust degree of fitting according to rank of Vandermonde matrix
1239  int ncols_sub = ncols;
1240  while( rank < ncols_sub )
1241  {
1242  --degree;
1243  if( degree == 0 )
1244  {
1245  // surface is flat, return 0
1246  *degree_out = *degree_qr = degree;
1247  for( int i = 0; i < npts2fit; ++i )
1248  {
1249  for( int k = 0; k < ndim; ++k )
1250  {
1251  bs[k * npts2fit + i] = 0;
1252  }
1253  }
1254  return;
1255  }
1256  else
1257  {
1258  ncols_sub = ( ( ( degree + 2 ) * ( degree + 1 ) ) >> 1 ) - interp;
1259  }
1260  }
1261  *degree_qr = degree;
1262 
1263  // std::cout << "degree_qr: " << *degree_qr << std::endl;
1264 
1265  /* DBG
1266  * std::cout<<"before Qtb"<<std::endl;
1267  std::cout<<std::endl;
1268  std::cout<<"bs = "<<std::endl;
1269  std::cout<<std::endl;
1270  for (int k=0; k< ndim; k++){
1271  for (int j=0; j<npts2fit; ++j){
1272  std::cout<<" "<<bs[npts2fit*k+j]<<std::endl;
1273  }
1274  }
1275  std::cout<<std::endl;*/
1276 
1277  // step 7. compute Q'b
1278  DGMSolver::compute_qtransposeB( npts2fit, ncols_sub, &( V[0] ), ndim, bs );
1279 
1280  /* DBG
1281  * std::cout<<"after Qtb"<<std::endl;
1282  std::cout<<"bs = "<<std::endl;
1283  std::cout<<std::endl;
1284  for (int k=0; k< ndim; k++){
1285  for (int j=0; j<npts2fit; ++j){
1286  std::cout<<" "<<bs[npts2fit*k+j]<<std::endl;
1287  }
1288  }
1289  std::cout<<std::endl;*/
1290 
1291  // step 8. perform backward substitution and scale the solution
1292  // assign diagonals of V
1293  for( int i = 0; i < ncols_sub; ++i )
1294  {
1295  V[i * npts2fit + i] = D[i];
1296  }
1297 
1298  // backsolve
1299  if( safeguard )
1300  {
1301  // for debug
1302  // std::cout << "ts size " << ts.size() << std::endl;
1303  DGMSolver::backsolve_polyfit_safeguarded( 2, degree, interp, npts2fit, ncols_sub, &( V[0] ), ndim, bs,
1304  &( ts[0] ), degree_out );
1305  }
1306  else
1307  {
1308  DGMSolver::backsolve( npts2fit, ncols_sub, &( V[0] ), 1, bs, &( ts[0] ) );
1309  *degree_out = degree;
1310  }
1311  /*if(V_init){
1312  delete [] V_init;
1313  }else{
1314  delete [] V;
1315  }*/
1316 }

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

Referenced by polyfit3d_surf_get_coeff().

◆ eval_vander_univar_cmf()

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

Form and solve Vandermonde system of single-variables.

Definition at line 1434 of file HiReconstruction.cpp.

1443 {
1444  // step 1. determine degree of polynomials to fit according to number of points
1445  int ncols = degree + 1 - interp;
1446  while( npts2fit < ncols && degree >= 1 )
1447  {
1448  --degree;
1449  ncols = degree + 1 - interp;
1450  }
1451  if( !degree )
1452  {
1453  if( interp )
1454  {
1455  for( int icol = 0; icol < ndim; ++icol )
1456  {
1457  bs[icol * npts2fit] = 0;
1458  }
1459  }
1460  for( int irow = 1; irow < npts2fit; ++irow )
1461  {
1462  for( int icol = 0; icol < ndim; ++icol )
1463  {
1464  bs[icol * npts2fit + irow] = 0;
1465  }
1466  }
1467  *degree_out = 0;
1468  return;
1469  }
1470  // step 2. construct Vandermonde matrix
1471  std::vector< double > V; // V(npts2fit*(ncols+interp));
1472  DGMSolver::gen_vander_multivar( npts2fit, 1, us, degree, V );
1473 
1474  if( interp )
1475  {
1476  V.erase( V.begin(), V.begin() + npts2fit );
1477  }
1478 
1479  // step 3. scale rows with respect to weights
1480  for( int i = 0; i < npts2fit; ++i )
1481  {
1482  for( int j = 0; j < ncols; ++j )
1483  {
1484  V[j * npts2fit + i] *= ws[i];
1485  }
1486  for( int k = 0; k < ndim; ++k )
1487  {
1488  bs[k * npts2fit + i] *= ws[i];
1489  }
1490  }
1491 
1492  // step 4. scale columns to reduce condition number
1493  std::vector< double > ts( ncols );
1494  DGMSolver::rescale_matrix( npts2fit, ncols, &( V[0] ), &( ts[0] ) );
1495 
1496  // step 5. perform Householder QR factorization
1497  std::vector< double > D( ncols );
1498  int rank;
1499  DGMSolver::qr_polyfit_safeguarded( npts2fit, ncols, &( V[0] ), &( D[0] ), &rank );
1500 
1501  // step 6. adjust degree of fitting
1502  int ncols_sub = ncols;
1503  while( rank < ncols_sub )
1504  {
1505  --degree;
1506  if( !degree )
1507  {
1508  // matrix is singular, consider curve on flat plane passing center and orthogonal to
1509  // tangent line
1510  *degree_out = 0;
1511  for( int i = 0; i < npts2fit; ++i )
1512  {
1513  for( int k = 0; k < ndim; ++k )
1514  {
1515  bs[k * npts2fit + i] = 0;
1516  }
1517  }
1518  }
1519  ncols_sub = degree + 1 - interp;
1520  }
1521 
1522  // step 7. compute Q'*bs
1523  DGMSolver::compute_qtransposeB( npts2fit, ncols_sub, &( V[0] ), ndim, bs );
1524 
1525  // step 8. perform backward substitution and scale solutions
1526  // assign diagonals of V
1527  for( int i = 0; i < ncols_sub; ++i )
1528  {
1529  V[i * npts2fit + i] = D[i];
1530  }
1531  // backsolve
1532  if( safeguard )
1533  {
1534  DGMSolver::backsolve_polyfit_safeguarded( 1, degree, interp, npts2fit, ncols, &( V[0] ), ndim, bs, ws,
1535  degree_out );
1536  }
1537  else
1538  {
1539  DGMSolver::backsolve( npts2fit, ncols_sub, &( V[0] ), ndim, bs, &( ts[0] ) );
1540  *degree_out = degree;
1541  }
1542 }

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

613 {
614  if( !_hasfittings )
615  {
616  return false;
617  }
618  else
619  {
620  int index = _verts2rec.index( vid );
621  if( -1 == index )
622  {
623  std::cout << "Input vertex is not locally hosted vertex in this mesh set" << std::endl;
624  return false;
625  }
626  geomtype = _geom;
627  if( HISURFACE == _geom )
628  {
629  coords.insert( coords.end(), _local_coords.begin() + 9 * index, _local_coords.begin() + 9 * index + 9 );
630  degree_out = _degrees_out[index];
631  interp = _interps[index];
632  int ncoeffs = ( degree_out + 2 ) * ( degree_out + 1 ) >> 1;
633  size_t istr = _vertID2coeffID[index];
634  coeffs.insert( coeffs.end(), _local_fit_coeffs.begin() + istr, _local_fit_coeffs.begin() + istr + ncoeffs );
635  }
636  else if( HI3DCURVE == _geom )
637  {
638  coords.insert( coords.end(), _local_coords.begin() + 3 * index, _local_coords.begin() + 3 * index + 3 );
639  degree_out = _degrees_out[index];
640  interp = _interps[index];
641  int ncoeffs = 3 * ( degree_out + 1 );
642  size_t istr = _vertID2coeffID[index];
643  coeffs.insert( coeffs.end(), _local_fit_coeffs.begin() + istr, _local_fit_coeffs.begin() + istr + ncoeffs );
644  }
645  return true;
646  }
647 }

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

◆ get_normals_surf()

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

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

Parameters
vertshRange, EntityHandles of vertices
nrmsPointer of array of doubles, size = 3*vertsh.size()

Definition at line 903 of file HiReconstruction.cpp.

904 {
906  if( _hasderiv )
907  {
908  size_t id = 0;
909  for( Range::iterator ivert = vertsh.begin(); ivert != vertsh.end(); ++ivert, ++id )
910  {
911  int index = _verts2rec.index( *ivert );
912 #ifdef MOAB_HAVE_MPI
913  if( -1 == index )
914  {
915  // ghost vertex
916  error = average_vertex_normal( *ivert, nrms + 3 * id );MB_CHK_ERR( error );
917  }
918  else
919  {
920  nrms[3 * id] = _local_coords[9 * index + 6];
921  nrms[3 * id + 1] = _local_coords[9 * index + 7];
922  nrms[3 * id + 2] = _local_coords[9 * index + 8];
923  }
924 #else
925  assert( -1 != index );
926  nrms[3 * id] = _local_coords[9 * index + 6];
927  nrms[3 * id + 1] = _local_coords[9 * index + 7];
928  nrms[3 * id + 2] = _local_coords[9 * index + 8];
929 #endif
930  }
931  }
932  else
933  {
934  size_t id = 0;
935  for( Range::iterator ivert = vertsh.begin(); ivert != vertsh.end(); ++ivert, ++id )
936  {
937  error = average_vertex_normal( *ivert, nrms + 3 * id );MB_CHK_ERR( error );
938  }
939  }
940  return error;
941 }

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

Referenced by polyfit3d_walf_surf_vertex().

◆ get_tangents_curve()

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

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

Parameters
vertshRange, EntityHandles of vertices
tangsPointer of array of doubles, size = 3*vertsh.size()

Definition at line 990 of file HiReconstruction.cpp.

991 {
993  if( _hasderiv )
994  {
995  size_t id = 0;
996  for( Range::iterator ivert = vertsh.begin(); ivert != vertsh.end(); ++ivert, ++id )
997  {
998  int index = _verts2rec.index( *ivert );
999 #ifdef MOAB_HAVE_MPI
1000  if( -1 != index )
1001  {
1002  tangs[3 * id] = _local_coords[3 * index];
1003  tangs[3 * id + 1] = _local_coords[3 * index + 1];
1004  tangs[3 * id + 2] = _local_coords[3 * index + 2];
1005  }
1006  else
1007  {
1008  error = average_vertex_tangent( *ivert, tangs + 3 * id );MB_CHK_ERR( error );
1009  }
1010 #else
1011  assert( -1 != index );
1012  tangs[3 * id] = _local_coords[3 * index];
1013  tangs[3 * id + 1] = _local_coords[3 * index + 1];
1014  tangs[3 * id + 2] = _local_coords[3 * index + 2];
1015 #endif
1016  }
1017  }
1018  else
1019  {
1020  size_t id = 0;
1021  for( Range::iterator ivert = vertsh.begin(); ivert != vertsh.end(); ++ivert, ++id )
1022  {
1023  error = average_vertex_tangent( *ivert, tangs + 3 * id );MB_CHK_ERR( error );
1024  }
1025  }
1026  return error;
1027 }

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

Referenced by polyfit3d_walf_curve_vertex().

◆ hiproj_walf_around_vertex()

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

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

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

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

472 {
473  if( !_hasfittings )
474  {
475  MB_SET_ERR( MB_FAILURE, "There is no existing fitting results" );
476  }
477  else if( -1 == _verts2rec.index( vid ) )
478  {
479  std::ostringstream convert;
480  convert << vid;
481  std::string VID = convert.str();
482  MB_SET_ERR( MB_FAILURE, "There is no existing fitting results for vertex " + VID );
483  }
485  // get center of local coordinates system
486  double local_origin[3];
487  error = mbImpl->get_coords( &vid, 1, local_origin );MB_CHK_ERR( error );
488  // get local fitting parameters
489  int index = _verts2rec.index( vid );
490  bool interp = _interps[index];
491  int local_deg = _degrees_out[index];
492  double *uvw_coords, *local_coeffs;
493  if( _geom == HISURFACE )
494  {
495  uvw_coords = &( _local_coords[9 * index] );
496  // int ncoeffs = (local_deg+2)*(local_deg+1)>>1;
497  size_t istr = _vertID2coeffID[index];
498  local_coeffs = &( _local_fit_coeffs[istr] );
499  walf3d_surf_vertex_eval( local_origin, uvw_coords, local_deg, local_coeffs, interp, npts2fit, coords2fit,
500  hiproj_new );
501  }
502  else if( _geom == HI3DCURVE )
503  {
504  uvw_coords = &( _local_coords[3 * index] );
505  size_t istr = _vertID2coeffID[index];
506  local_coeffs = &( _local_fit_coeffs[istr] );
507  walf3d_curve_vertex_eval( local_origin, uvw_coords, local_deg, local_coeffs, interp, npts2fit, coords2fit,
508  hiproj_new );
509  }
510  return error;
511 }

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

Referenced by hiproj_walf_in_element().

◆ hiproj_walf_in_element()

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

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

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

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

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

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  }
63 #else
64  ahf = NULL;
65 #endif
66 
67  // error = ahf->get_entity_ranges(_inverts,_inedges,_infaces,_incells); MB_CHK_ERR(error);
72  if( _inedges.size() && _infaces.empty() && _incells.empty() )
73  {
74  _dim = 1;
75  _MAXPNTS = 13;
76  }
77  else if( _infaces.size() && _incells.empty() )
78  {
79  _dim = 2;
80  _MAXPNTS = 128;
81  }
82  else
83  {
84  MB_SET_ERR( MB_FAILURE, "Encountered a non-manifold mesh or a mesh with volume elements" );
85  }
86 
87  // get locally hosted vertices by filtering pstatus
88 #ifdef MOAB_HAVE_MPI
89  if( pcomm )
90  {
92  }
93  else
94  {
96  }
97 #else
99 #endif
100  _nv2rec = _verts2rec.size();
101 
102  if( recwhole )
103  {
104  // compute normals(surface) or tangent vector(curve) for all locally hosted vertices
105  if( 2 == _dim )
106  {
108  }
109  else if( 1 == _dim )
110  {
112  }
113  else
114  {
115  MB_SET_ERR( MB_FAILURE, "Unknow space dimension" );
116  }
117  _hasderiv = true;
118  }
119  return error;
120 }

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

Referenced by HiReconstruction().

◆ initialize_3Dcurve_geom() [1/2]

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

Definition at line 780 of file HiReconstruction.cpp.

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

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

Referenced by reconstruct3D_curve_geom().

◆ initialize_3Dcurve_geom() [2/2]

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

Definition at line 802 of file HiReconstruction.cpp.

803 {
804  if( !_hasderiv )
805  {
807  _hasderiv = true;
808  }
809  if( !_hasfittings )
810  {
811  assert( _nv2rec == npts );
812  _degrees_out.assign( _nv2rec, 0 );
813  _interps.assign( _nv2rec, false );
814  _vertID2coeffID.reserve( _nv2rec );
815  size_t index = 0;
816  for( size_t i = 0; i < _nv2rec; ++i )
817  {
818  _vertID2coeffID[i] = index;
819  index += 3 * ( degrees[i] + 1 );
820  }
821  _local_fit_coeffs.assign( index, 0 );
822  _initfittings = true;
823  }
824 }

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

◆ initialize_surf_geom() [1/2]

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

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

Definition at line 734 of file HiReconstruction.cpp.

735 {
736  if( !_hasderiv )
737  {
739  _hasderiv = true;
740  }
741  if( !_initfittings )
742  {
743  int ncoeffspv = ( degree + 2 ) * ( degree + 1 ) / 2;
744  _degrees_out.assign( _nv2rec, 0 );
745  _interps.assign( _nv2rec, false );
746  _vertID2coeffID.resize( _nv2rec );
747  _local_fit_coeffs.assign( _nv2rec * ncoeffspv, 0 );
748  for( size_t i = 0; i < _nv2rec; ++i )
749  {
750  _vertID2coeffID[i] = i * ncoeffspv;
751  }
752  _initfittings = true;
753  }
754 }

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

Referenced by reconstruct3D_surf_geom().

◆ initialize_surf_geom() [2/2]

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

Definition at line 756 of file HiReconstruction.cpp.

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

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

◆ obtain_nring_ngbvs()

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

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

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

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

673 {
675  std::deque< EntityHandle > todo;
676  todo.push_back( vid );
677  ngbvs.insert( vid );
678  EntityHandle pre, nxt;
679  for( int i = 1; i <= ring; ++i )
680  {
681  int count = todo.size();
682  while( count )
683  {
684  EntityHandle center = todo.front();
685  todo.pop_front();
686  --count;
687  std::vector< EntityHandle > adjents;
689  for( size_t j = 0; j < adjents.size(); ++j )
690  {
691  std::vector< EntityHandle > elemconn;
692  error = mbImpl->get_connectivity( &adjents[j], 1, elemconn );MB_CHK_ERR( error );
693  int nvpe = elemconn.size();
694  for( int k = 0; k < nvpe; ++k )
695  {
696  if( elemconn[k] == center )
697  {
698  pre = k == 0 ? elemconn[nvpe - 1] : elemconn[k - 1];
699  nxt = elemconn[( k + 1 ) % nvpe];
700  if( ngbvs.find( pre ) == ngbvs.end() )
701  {
702  ngbvs.insert( pre );
703  todo.push_back( pre );
704  }
705  if( ngbvs.find( nxt ) == ngbvs.end() )
706  {
707  ngbvs.insert( nxt );
708  todo.push_back( nxt );
709  }
710  break;
711  }
712  }
713  }
714  }
715  if( _MAXPNTS <= (int)ngbvs.size() )
716  {
717  // obtain enough points
718  return error;
719  }
720  if( !todo.size() )
721  {
722  // current ring cannot introduce any points, return incase deadlock
723  return error;
724  }
725  if( ( i == ring ) && ( minpnts > (int)ngbvs.size() ) )
726  {
727  // reach maximum ring but not enough points
728  ++ring;
729  }
730  }
731  return error;
732 }

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

1329 {
1330  if( !nverts )
1331  {
1332  return;
1333  }
1334  // step 1. compute local coordinates system
1335  double tang[3] = { ngbtangs[0], ngbtangs[1], ngbtangs[2] };
1336  if( coords && ncoords > 2 )
1337  {
1338  coords[0] = tang[0];
1339  coords[1] = tang[1];
1340  coords[2] = tang[2];
1341  }
1342  if( !coeffs || !ncoeffs )
1343  {
1344  return;
1345  }
1346  else
1347  {
1348  assert( ncoeffs >= 3 * ( degree + 1 ) );
1349  }
1350  // step 2. project onto local coordinate system
1351  int npts2fit = nverts - interp;
1352  if( !npts2fit )
1353  {
1354  *degree_out = 0;
1355  for( int i = 0; i < ncoeffs; ++i )
1356  {
1357  coeffs[0] = 0;
1358  }
1359  return;
1360  }
1361  std::vector< double > us( npts2fit ); // double *us = new double[npts2fit];
1362  std::vector< double > bs( npts2fit * 3 ); // double *bs = new double[npts2fit*3];
1363  double uu[3];
1364  for( int i = interp; i < nverts; ++i )
1365  {
1366  int k = i - interp;
1367  DGMSolver::vec_linear_operation( 3, 1, ngbcors + 3 * i, -1, ngbcors, uu );
1368  us[k] = DGMSolver::vec_innerprod( 3, uu, tang );
1369  bs[k] = uu[0];
1370  bs[npts2fit + k] = uu[1];
1371  bs[2 * npts2fit + k] = uu[2];
1372  }
1373 
1374  // step 3. copmute weights
1375  std::vector< double > ws( npts2fit );
1376  int nzeros = compute_weights( npts2fit, 1, &( us[0] ), nverts, ngbtangs, degree, _MINEPS, &( ws[0] ) );
1377  assert( nzeros <= npts2fit );
1378 
1379  // step 4. adjust according to number of zero-weights
1380  if( nzeros )
1381  {
1382  if( nzeros == npts2fit )
1383  {
1384  // singular case
1385  *degree_out = 0;
1386  for( int i = 0; i < ncoeffs; ++i )
1387  {
1388  coeffs[i] = 0;
1389  }
1390  return;
1391  }
1392  int npts_new = npts2fit - nzeros;
1393  std::vector< double > bs_temp( 3 * npts_new );
1394  int index = 0;
1395  for( int i = 0; i < npts2fit; ++i )
1396  {
1397  if( ws[i] )
1398  {
1399  if( i > index )
1400  {
1401  us[index] = us[i];
1402  ws[index] = ws[i];
1403  }
1404  bs_temp[index] = bs[i];
1405  bs_temp[index + npts_new] = bs[i + npts2fit];
1406  bs_temp[index + 2 * npts_new] = bs[i + 2 * npts2fit];
1407  ++index;
1408  }
1409  }
1410  assert( index == npts_new );
1411  npts2fit = npts_new;
1412  us.resize( index );
1413  ws.resize( index );
1414  bs = bs_temp;
1415  // destroy bs_temp;
1416  std::vector< double >().swap( bs_temp );
1417  }
1418 
1419  // step 5. fitting
1420  eval_vander_univar_cmf( npts2fit, &( us[0] ), 3, &( bs[0] ), degree, &( ws[0] ), interp, safeguard, degree_out );
1421  // step 6. write results to output pointers
1422  int ncoeffs_out_pvpd = *degree_out + 1;
1423  assert( ncoeffs >= 3 * ncoeffs_out_pvpd );
1424  // write to coeffs consecutively, bs's size is not changed by eval_vander_univar_cmf
1425  coeffs[0] = coeffs[ncoeffs_out_pvpd] = coeffs[2 * ncoeffs_out_pvpd] = 0;
1426  for( int i = 0; i < ncoeffs_out_pvpd - interp; ++i )
1427  {
1428  coeffs[i + interp] = bs[i];
1429  coeffs[i + interp + ncoeffs_out_pvpd] = bs[i + npts2fit];
1430  coeffs[i + interp + 2 * ncoeffs_out_pvpd] = bs[i + 2 * npts2fit];
1431  }
1432 }

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

Referenced by polyfit3d_walf_curve_vertex().

◆ polyfit3d_surf_get_coeff()

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

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

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

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

1046 {
1047  if( nverts <= 0 )
1048  {
1049  return;
1050  }
1051 
1052  // std::cout << "npnts in initial stencil = " << nverts << std::endl;
1053  // std::cout << "centered at (" << ngbcoords[0] << "," << ngbcoords[1] << "," << ngbcoords[2] <<
1054  // ")" << std::endl;
1055 
1056  // step 1. copmute local coordinate system
1057  double nrm[3] = { ngbnrms[0], ngbnrms[1], ngbnrms[2] }, tang1[3] = { 0, 0, 0 }, tang2[3] = { 0, 0, 0 };
1058  if( fabs( nrm[0] ) > fabs( nrm[1] ) && fabs( nrm[0] ) > fabs( nrm[2] ) )
1059  {
1060  tang1[1] = 1.0;
1061  }
1062  else
1063  {
1064  tang1[0] = 1.0;
1065  }
1066 
1067  DGMSolver::vec_projoff( 3, tang1, nrm, tang1 );
1068 #ifndef NDEBUG
1069  assert( DGMSolver::vec_normalize( 3, tang1, tang1 ) );
1070 #endif
1071  DGMSolver::vec_crossprod( nrm, tang1, tang2 );
1072  if( 9 <= ncoords && coords )
1073  {
1074  coords[0] = tang1[0];
1075  coords[1] = tang1[1];
1076  coords[2] = tang1[2];
1077  coords[3] = tang2[0];
1078  coords[4] = tang2[1];
1079  coords[5] = tang2[2];
1080  coords[6] = nrm[0];
1081  coords[7] = nrm[1];
1082  coords[8] = nrm[2];
1083  }
1084  if( !ncoeffs || !coeffs )
1085  {
1086  return;
1087  }
1088  else
1089  {
1090  assert( ncoeffs >= ( degree + 2 ) * ( degree + 1 ) / 2 );
1091  }
1092 
1093  // step 2. project onto local coordinates system
1094  int npts2fit = nverts - interp;
1095  if( 0 == npts2fit )
1096  {
1097  *degree_out = *degree_pnt = *degree_qr = 0;
1098  for( int i = 0; i < ncoeffs; ++i )
1099  {
1100  coeffs[i] = 0;
1101  }
1102  // coeffs[0] = 0;
1103  return;
1104  }
1105  std::vector< double > us( npts2fit * 2 ); // double *us = new double[npts2fit*2];
1106  std::vector< double > bs( npts2fit ); // double *bs = new double[npts2fit];
1107  for( int i = interp; i < nverts; ++i )
1108  {
1109  int k = i - interp;
1110  double uu[3];
1111  DGMSolver::vec_linear_operation( 3, 1, ngbcoords + 3 * i, -1, ngbcoords, uu );
1112  us[k * 2] = DGMSolver::vec_innerprod( 3, tang1, uu );
1113  us[k * 2 + 1] = DGMSolver::vec_innerprod( 3, tang2, uu );
1114  bs[k] = DGMSolver::vec_innerprod( 3, nrm, uu );
1115  }
1116 
1117  // step 3. compute weights
1118  std::vector< double > ws( npts2fit ); // double *ws = new double[npts2fit];
1119  int nzeros = compute_weights( npts2fit, 2, &( us[0] ), nverts, ngbnrms, degree, _MINEPS, &( ws[0] ) );
1120 
1121  // step 4. adjust according to zero-weights
1122  if( nzeros )
1123  {
1124  if( nzeros == npts2fit )
1125  {
1126  *degree_out = *degree_pnt = *degree_qr = 0;
1127  for( int i = 0; i < ncoeffs; ++i )
1128  {
1129  coeffs[i] = 0;
1130  }
1131  // coeffs[0] = 0;
1132  return;
1133  }
1134  int index = 0;
1135  for( int i = 0; i < npts2fit; ++i )
1136  {
1137  if( ws[i] )
1138  {
1139  if( i > index )
1140  {
1141  us[index * 2] = us[i * 2];
1142  us[index * 2 + 1] = us[i * 2 + 1];
1143  bs[index] = bs[i];
1144  ws[index] = ws[i];
1145  }
1146  ++index;
1147  }
1148  }
1149  npts2fit -= nzeros;
1150  assert( index == npts2fit );
1151  us.resize( npts2fit * 2 );
1152  bs.resize( npts2fit );
1153  ws.resize( npts2fit );
1154  /*us = realloc(us,npts2fit*2*sizeof(double));
1155  bs = realloc(bs,npts2fit*sizeof(double));
1156  ws = realloc(ws,npts2fit*sizeof(double));*/
1157  }
1158 
1159  // std::cout<<"npnts after weighting = "<<npts2fit<<std::endl;
1160 
1161  // step 5. fitting
1162  eval_vander_bivar_cmf( npts2fit, &( us[0] ), 1, &( bs[0] ), degree, &( ws[0] ), interp, safeguard, degree_out,
1163  degree_pnt, degree_qr );
1164 
1165  // step 6. organize output
1166  int ncoeffs_out = ( *degree_out + 2 ) * ( *degree_out + 1 ) / 2;
1167  assert( ncoeffs_out <= ncoeffs );
1168  coeffs[0] = 0;
1169  for( int j = 0; j < ncoeffs_out - interp; ++j )
1170  {
1171  coeffs[j + interp] = bs[j];
1172  }
1173  // delete [] us; delete [] bs; delete [] ws;
1174 }

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

Referenced by polyfit3d_walf_surf_vertex().

◆ polyfit3d_walf_curve_vertex()

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

Construct vertex based polynomial fitting on a curve mesh.

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

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

354 {
356  int ring = estimate_num_rings( degree, interp );
357  // get n-ring neighbors
358  Range ngbvs;
359  error = obtain_nring_ngbvs( vid, ring, minpnts, ngbvs );MB_CHK_ERR( error );
360  // get coordinates
361  size_t nverts = ngbvs.size();
362  assert( nverts );
363  double* ngbcoords = new double[nverts * 3];
364  error = mbImpl->get_coords( ngbvs, ngbcoords );MB_CHK_ERR( error );
365  // get tangent vectors
366  double* ngbtangs = new double[nverts * 3];
367  error = get_tangents_curve( ngbvs, ngbtangs );MB_CHK_ERR( error );
368  // switch vid to first one
369  int index = ngbvs.index( vid );
370  assert( index != -1 );
371  std::swap( ngbcoords[0], ngbcoords[3 * index] );
372  std::swap( ngbcoords[1], ngbcoords[3 * index + 1] );
373  std::swap( ngbcoords[2], ngbcoords[3 * index + 2] );
374  std::swap( ngbtangs[0], ngbtangs[3 * index] );
375  std::swap( ngbtangs[1], ngbtangs[3 * index + 1] );
376  std::swap( ngbtangs[2], ngbtangs[3 * index + 2] );
377  // local WLS fittings
378  polyfit3d_curve_get_coeff( nverts, ngbcoords, ngbtangs, degree, interp, safeguard, ncoords, coords, ncoeffs, coeffs,
379  degree_out );
380  delete[] ngbcoords;
381  delete[] ngbtangs;
382  return error;
383 }

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

Referenced by reconstruct3D_curve_geom().

◆ polyfit3d_walf_surf_vertex()

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

Construct vertex based polynomial fitting on a surface mesh.

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

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

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

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

Referenced by reconstruct3D_surf_geom().

◆ reconstruct3D_curve_geom() [1/2]

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

Reconstruct a high order curve on given curve mesh.

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

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

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

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

◆ reconstruct3D_curve_geom() [2/2]

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

Reconstruct a high order curve on given curve mesh.

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

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

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

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

◆ reconstruct3D_surf_geom() [1/2]

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

Reconstruct a high order surface on given surface mesh.

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

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

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

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

◆ reconstruct3D_surf_geom() [2/2]

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

Reconstruct a high order surface on given surface mesh.

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

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

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

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

◆ vertex_get_incident_elements()

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

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

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

Parameters
vidEntityHandle of vertex
elemdimInteger, dimension of elements incidented in vid
adjentsvector<EntityHandle>, container which push incident elements in

Definition at line 658 of file HiReconstruction.cpp.

661 {
663  assert( elemdim == _dim );
664 #ifdef HIREC_USE_AHF
665  error = ahf->get_up_adjacencies( vid, elemdim, adjents );MB_CHK_ERR( error );
666 #else
667  error = mbImpl->get_adjacencies( &vid, 1, elemdim, false, adjents );MB_CHK_ERR( error );
668 #endif
669  return error;
670 }

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

577 {
578  assert( local_origin && local_coords && local_coeffs );
579  int ncoeffspvpd = local_deg + 1;
580  for( int i = 0; i < npts2fit; ++i )
581  {
582  // get the vector from center to current point, project to tangent line
583  double vec[3], ans[3] = { 0, 0, 0 };
584  DGMSolver::vec_linear_operation( 3, 1, coords2fit + 3 * i, -1, local_origin, vec );
585  double u = DGMSolver::vec_innerprod( 3, local_coords, vec );
586  // evaluate polynomials
587  if( !interp )
588  {
589  ans[0] = local_coeffs[0];
590  ans[1] = local_coeffs[ncoeffspvpd];
591  ans[2] = local_coeffs[2 * ncoeffspvpd];
592  }
593  double uk = 1; // degree_out and degree different, stored in columnwise contiguously
594  for( int j = 1; j < ncoeffspvpd; ++j )
595  {
596  uk *= u;
597  ans[0] += uk * local_coeffs[j];
598  ans[1] += uk * local_coeffs[j + ncoeffspvpd];
599  ans[2] += uk * local_coeffs[j + 2 * ncoeffspvpd];
600  }
601  hiproj_new[3 * i] = ans[0] + local_origin[0];
602  hiproj_new[3 * i + 1] = ans[1] + local_origin[1];
603  hiproj_new[3 * i + 2] = ans[2] + local_origin[2];
604  }
605 }

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

521 {
522  double xaxis[3], yaxis[3], zaxis[3];
523  for( int i = 0; i < 3; ++i )
524  {
525  xaxis[i] = local_coords[i];
526  yaxis[i] = local_coords[3 + i];
527  zaxis[i] = local_coords[6 + i];
528  }
529  // double *basis = new double[(local_deg+2)*(local_deg+1)/2-1];
530  std::vector< double > basis( ( local_deg + 2 ) * ( local_deg + 1 ) / 2 - 1 );
531  for( int i = 0; i < npts2fit; ++i )
532  {
533  double local_pos[3];
534  for( int j = 0; j < 3; ++j )
535  {
536  local_pos[j] = coords2fit[3 * i + j] - local_origin[j];
537  }
538  double u, v, height = 0;
539  u = DGMSolver::vec_innerprod( 3, local_pos, xaxis );
540  v = DGMSolver::vec_innerprod( 3, local_pos, yaxis );
541  basis[0] = u;
542  basis[1] = v;
543  int l = 1;
544  for( int k = 2; k <= local_deg; ++k )
545  {
546  ++l;
547  basis[l] = u * basis[l - k];
548  for( int id = 0; id < k; ++id )
549  {
550  ++l;
551  basis[l] = basis[l - k - 1] * v;
552  }
553  }
554  if( !interp )
555  {
556  height = local_coeffs[0];
557  }
558  for( int p = 0; p <= l; ++p )
559  {
560  height += local_coeffs[p + 1] * basis[p];
561  }
562  hiproj_new[3 * i] = local_origin[0] + u * xaxis[0] + v * yaxis[0] + height * zaxis[0];
563  hiproj_new[3 * i + 1] = local_origin[1] + u * xaxis[1] + v * yaxis[1] + height * zaxis[1];
564  hiproj_new[3 * i + 2] = local_origin[2] + u * xaxis[2] + v * yaxis[2] + height * zaxis[2];
565  }
566  // delete [] basis;
567 }

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: