Loading [MathJax]/extensions/tex2jax.js
Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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 ); 22  ErrorCode error; 23  _MINEPS = 1e-12; 24  _dim = 0; 25  _hasfittings = false; 26  _hasderiv = false; 27 #ifdef MOAB_HAVE_MPI 28  if( !pcomm ) 29  { 30  pcomm = moab::ParallelComm::get_pcomm( mbImpl, 0 ); 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 { 843  ErrorCode error; 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 { 945  ErrorCode error; 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  } 893  ErrorCode error; 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  } 980  ErrorCode error; 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 { 905  ErrorCode error = MB_SUCCESS; 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 { 992  ErrorCode error = MB_SUCCESS; 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  } 484  ErrorCode error; 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 ); 396  ErrorCode error; 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 { 53  ErrorCode error; 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  { 60  return MB_MEMORY_ALLOCATION_FAILED; 61  } 62  error = ahf->initialize();MB_CHK_ERR( error ); 63 #else 64  ahf = NULL; 65 #endif 66  67  // error = ahf->get_entity_ranges(_inverts,_inedges,_infaces,_incells); MB_CHK_ERR(error); 68  error = mbImpl->get_entities_by_dimension( _mesh2rec, 0, _inverts );MB_CHK_ERR( error ); 69  error = mbImpl->get_entities_by_dimension( _mesh2rec, 1, _inedges );MB_CHK_ERR( error ); 70  error = mbImpl->get_entities_by_dimension( _mesh2rec, 2, _infaces );MB_CHK_ERR( error ); 71  error = mbImpl->get_entities_by_dimension( _mesh2rec, 3, _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  { 91  error = pcomm->filter_pstatus( _inverts, PSTATUS_GHOST, PSTATUS_NOT, -1, &_verts2rec );MB_CHK_ERR( error ); 92  } 93  else 94  { 95  _verts2rec = _inverts; 96  } 97 #else 98  _verts2rec = _inverts; 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  { 107  compute_average_vertex_normals_surf(); 108  } 109  else if( 1 == _dim ) 110  { 111  compute_average_vertex_tangents_curve(); 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  { 784  compute_average_vertex_tangents_curve(); 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  { 806  compute_average_vertex_tangents_curve(); 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  { 738  compute_average_vertex_normals_surf(); 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  { 760  compute_average_vertex_normals_surf(); 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 { 674  ErrorCode error; 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; 688  error = vertex_get_incident_elements( center, _dim, adjents );MB_CHK_ERR( error ); 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 { 355  ErrorCode error; 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 ); 306  ErrorCode error; 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 ); 233  ErrorCode error; 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 ); 260  ErrorCode error; 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 ); 140  ErrorCode error; 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  } 197  ErrorCode error; 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 { 662  ErrorCode error; 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: