#include <FaceOffset_3.h>

Inheritance diagram for FaceOffset_3:

Collaboration diagram for FaceOffset_3:

Public Member Functions
	FaceOffset_3 ()
	Default constructor. More...

	FaceOffset_3 (Manifold wm, COM::Window buf)
	Construct an object from a window manifold and a buffer window. More...

virtual	~FaceOffset_3 ()

virtual double	time_stepping (const COM::Attribute spds, double dt, COM::Attribute disps, int *smoothed=NULL)
	Main entry of the algorithm. More...

void	set_wavefrontal_expansion (bool b)
	Set the wavefrontal switch. More...

bool	get_wavefrontal_expansion () const
	Obtain the wavefrontal switch. More...

void	set_normal_diffusion (char i)
	Set number of iterations for normal diffusion. More...

void	set_feature_layer (bool b)
	Set whether or not to use feature layer. More...

int	get_normal_diffuion () const
	Get number of iterations for normal diffusion. More...

void	set_eigen_threshold (double eps)
	Set the threshold of eigenvalues for splitting primary and null spaces. More...

double	get_eigen_threshold () const
	Get the threshold of eigenvalues for splitting primary and null spaces. More...

void	set_courant_constant (double c)
	Set the constant in front of CFL condition (0<=c<1) More...

double	get_courant_constant () const
	Get the constant in front of CFL condition (0<=c<=1) More...

void	set_weighting_scheme (int w)
	Weighting scheme. More...

int	get_weighting_scheme () const

void	set_fangle_weak (double psi)
	Set the threshold for weak-feature angle. More...

void	set_fangle_strong (double psi)
	Set the threshold for weak-feature angle. More...

void	set_fangle_turn (double psi)
	Set the threshold for weak-feature angle. More...

double	get_fangle_weak_radians () const
	Get the threshold for weak-feature angle. More...

double	get_fangle_strong_radians () const
	Get the threshold for strong-feature angle. More...

double	get_fangle_turn_radians () const
	Set the threshold for weak-feature angle. More...

void	set_smoother (int smoother)
	Set the choice of mesh smoother. More...

void	set_conserve (int b)
	Set mass conservation. More...

void	reset_default_parameters (bool b=false)

Public Member Functions inherited from Propagation_3
	Propagation_3 ()

virtual	~Propagation_3 ()

	Propagation_3 (Manifold wm, COM::Window buf)
	Construct an object from a window manifold. More...

virtual void	set_constraints (const COM::Attribute *cnstr_types)
	Set the types and directions of nodal constraints. More...

void	set_bounds (const COM::Attribute *bnd)
	Set the bounds. More...

void	bound_facial_speed (COM::Attribute *fa)

virtual void	enforce_nodal_constraints (COM::Attribute *du)
	Enforces the nodal constraints by projecting motion onto given direction. More...

virtual void	bound_nodal_motion (COM::Attribute *disps)

void	set_verbose (bool b)
	Set the verbose level. More...

Protected Member Functions
void	numerical_dissipation (COM::Attribute *nodal_buffer)
	Introduce numerical dissipation into equation. More...

void	compute_quadrics (double dt, const COM::Attribute spds, COM::Attribute rhs, COM::Attribute *predicted)
	Compute quadrics for every vertex. More...

void	compute_directions (COM::Attribute *b_offset, bool in_principle)
	Compute mean normals. More...

void	compute_volume_vector (const COM::Attribute disps, COM::Attribute bs)
	Compute volumn-vector. More...

bool	obtain_constrained_directions (COM::Attribute disps, COM::Attribute vcenters)
	Decompose propagation directions based on constraints. More...

double	rescale_displacements (COM::Attribute disps, COM::Attribute vcenters, int depth=0)
	Reduce time step based on stability limit, and rescale displacements by the reduced time step and diffusing expansion. More...

void	filter_and_identify_ridge_edges (bool filter_curves)
	Filter out isolated ridge vertices and identify ridge edges. More...

void	reclassify_ridge_vertices (const bool upgrade_corners, const bool upgrade_ridge, COM::Attribute neighbor_feas, COM::Attribute tr_attr, bool to_filter)

int	build_ridge_neighbor_list (const COM::Attribute maxtrnangv, const COM::Attribute mintrnangv, const std::vector< std::map< Edge_ID, double > > &edge_maps, const COM::Attribute maxranglev, const COM::Attribute minranglev, const std::vector< std::map< Edge_ID, double > > &rangle_maps, std::vector< ObscendSet > &obscends)
	Build the list of ridge nighbors and obtain a list of weak-ended vertices Return the number of obscended vertices. More...

int	append_obscure_ends (const COM::Attribute maxtrnangv, const COM::Attribute mintrnangv, const std::vector< std::map< Edge_ID, double > > &edge_maps, const COM::Attribute maxranglev, const COM::Attribute minranglev, const std::vector< std::map< Edge_ID, double > > &rangle_maps)

bool	filter_obscended_ridge (const std::vector< std::map< Edge_ID, double > > &edge_maps, const std::vector< std::map< Edge_ID, double > > &diangl_maps, const std::vector< ObscendSet > &obscends)
	Filter out weak-ended curves. More...

int	remove_obscure_curves (const std::vector< ObscendSet > &obscends)

void	mark_weak_vertices ()
	Mark vertices that are weak. More...

void	nulaplacian_smooth (const COM::Attribute vert_normals, const COM::Attribute tangranks, const COM::Attribute disps, COM::Attribute vcenters, COM::Attribute buf, COM::Attribute vert_weights_buf, const COM::Attribute *edge_weights=NULL)
	Redistribute smooth vertices by NuLaplacian smoothing. More...

void	balance_mass ()

void	distribute_volume_e2n (const COM::Attribute fvol, const COM::Attribute tranks, COM::Attribute *vvol)

void	update_face_volumes (const COM::Attribute fnormal, const COM::Attribute vdsps, COM::Attribute *fvol)

void	adjust_wavefrontal_motion (COM::Attribute *disps)

std::pair< double, double >	comp_wavefrontal_motion (double det, double w, double delta, const Vector_3 &disp_nz, const Vector_3 &nrm_nz)

int	eigen_analyze_vertex (Vector_3 A_io[3], Vector_3 &b_io, double angle_defect=0)

void	obtain_directions (Vector_3 es[3], const Vector_3 &lambdas, int nrank, Vector_3 &b_medial, Vector_3 *b_offset=NULL) const

void	obtain_face_offset (const Point_3 pnts, const double spd_ptr, const COM::Attribute spd, double dt, const Vector_3 dirs, COM::Element_node_enumerator &ene, Vector_3 &ns_nz, Point_3 &cnt, Vector_3 disps, Vector_3 ns)

void	get_tangents (const Vector_3 es, const Vector_3 &lambdas, const Vector_3 &mn, int nrank, Vector_3 vs[2], double scales[2], int is_strong=NULL) const

void	compute_anisotropic_vertex_centers (const COM::Attribute *disps_buf)

void	denoise_surface (const COM::Attribute disps_buf, COM::Attribute normal_motion)

void	get_primary_component (const Vector_3 &nrm, const Vector_3 es[], int trank, Vector_3 &prim) const

int	insert_boundary_edges (COM::Attribute *tr_attr)

void	update_vertex_centers ()

std::pair< double, double >	solve_quadratic_func (double a, double b, double c)

double	sign (double x)

Vector_3	compute_weighted_normal (Vector_3 es[3], const Vector_3 &lambdas, int trank, const Vector_3 &mean_nrm, const Vector_3 &face_nrm)

bool	is_acute_ridge (const Vector_3 *es, const Vector_3 &lambdas, const Vector_3 &mn, int nrank) const

bool	is_acute_corner (const Vector_3 *es, const Vector_3 &lambdas, const Vector_3 &mn, int nrank) const

bool	is_strong (const Vector_3 *es, const Vector_3 &lambdas, const Vector_3 &mn, int nrank) const

double	get_nonuniform_weight (const Vector_3 &lambdas, int trank)

Protected Member Functions inherited from Propagation_3
void	convert_constraints (const COM::Attribute ctypes_faces, COM::Attribute ctypes_nodes)
	Convert facial constraints to nodal constraints. More...

void	determine_constraint_boundary (const COM::Attribute ctypes_faces, COM::Attribute ctypes_bndry_edges, COM::Attribute *ctypes_bndry_nodes)
	Convert facial or panel constraints to nodal constraints. More...

Static Protected Member Functions
template<class FT >
static void	solve (const FT &a1, const FT &a2, const FT &b1, const FT &b2, const FT &c1, const FT &c2, FT &x, FT &y)

template<class FT >
static void	solve (const FT &a1, const FT &a2, const FT &a3, const FT &b1, const FT &b2, const FT &b3, const FT &c1, const FT &c2, const FT &c3, const FT &d1, const FT &d2, const FT &d3, FT &x, FT &y, FT &z)

static void	solve (const Vector_3 A[3], const Vector_3 &q, Vector_3 &x)

static Vector_2	proj (const Vector_3 &v, const Vector_3 &d1, const Vector_3 &d2)

static double	eval_angle (const Vector_3 &v1, const Vector_3 &v2)

static double	eval_angle (const Vector_2 &v1, const Vector_2 &v2)

static void	compute_eigenvectors (Vector_3 A[3], Vector_3 &lambdas)

Static Protected Member Functions inherited from Propagation_3
static void	get_constraint_directions (int type, int &ndirs, Vector_3 dirs[2])
	Get orthonormals of the constraints. More...

static void	enforce_nodal_constraint (int type, Vector_3 &du)
	Enforce constraint for a specific vector. More...

static void	bound_nodal_motion (const Point_3 &pnt, const double *bnd, Vector_3 &du, double eps=0)

static bool	check_spherical_bound (const Point_3 &pnt, const Point_3 &org, const double rad_min, const double rad_max, double eps=0)

static void	bound_spherical_disp (const Point_3 &pnt, const Point_3 &org, const double rad_min, const double rad_max, Vector_3 &disp, double eps=0)

static bool	check_radial_bound (const double x, const double y, const double bnd_min, const double bnd_max, double eps=0)

static void	bound_radial_disp (const double x, const double y, const double bnd_min, const double bnd_max, double &dx, double &dy, double eps=0)

static bool	check_axial_bound (const double x, const double bnd_min, const double bnd_max, double eps=0)

static void	bound_axial_disp (const double x, const double bnd_min, const double bnd_max, double &dx, double eps=0)

static bool	reached_nodal_bound (const Point_3 &pnt, const double *bnd, double eps=0)

static bool	in_bounding_box (const Point_3 &pnt, const Point_3 &lb, const Point_3 &ub)

static double	square (double x)

Protected Attributes
bool	_wf_expn

char	_wght_scheme

char	_nrm_dfsn

char	_smoother

double	_courant

double	_inv_sqrt_c

double	_dir_thres_weak

double	_dir_thres_strong

double	_tol_angle_weak

double	_tol_angle_strong

int	_tol_kstrong

double	_tol_kangle

double	_tol_eangle

double	_tol_cos_osta

double	_tol_turn_angle

double	_tol_ad

double	_eig_thres

double	_eig_weak_lb

double	_eig_weak_ub

bool	_conserv

bool	_is_strtd

bool	_feature_layer

COM::Attribute *	_As

COM::Attribute *	_boffset

COM::Attribute *	_bmedial

COM::Attribute *	_eigvalues

COM::Attribute *	_eigvecs

COM::Attribute *	_vnormals

COM::Attribute *	_facenormals

COM::Attribute *	_facecenters

COM::Attribute *	_faceareas

COM::Attribute *	_vcenters

COM::Attribute *	_tangranks

COM::Attribute *	_ctangranks

COM::Attribute *	_weak

COM::Attribute *	_strong

COM::Attribute *	_ridges

COM::Attribute *	_ridgeneighbors

COM::Attribute *	_scales

COM::Attribute *	_weights

std::vector< std::set< Edge_ID > >	_edges

Protected Attributes inherited from Propagation_3
Manifold *	_surf

SURF::Access_Mode	_mode

COM::Window *	_buf

int	_rank

bool	_verb

bool	_cnstr_set

bool	_bnd_set

COM::Attribute *	_cnstr_nodes

COM::Attribute *	_cnstr_faces

COM::Attribute *	_cnstr_bndry_edges

COM::Attribute *	_cnstr_bndry_nodes

COM::Attribute *	_cnstr_bound

std::vector< COM::Pane * >	_panes

Static Protected Attributes
static const double	pi = 3.14159265358979

Private Types
typedef std::map< Edge_ID, std::pair< int, int > >	ObscendSet

Detailed Description

Definition at line 39 of file FaceOffset_3.h.

Member Typedef Documentation

typedef std::map<Edge_ID,std::pair<int, int> > ObscendSet

private

Definition at line 40 of file FaceOffset_3.h.

Constructor & Destructor Documentation

FaceOffset_3 ( )

inline

Default constructor.

Definition at line 43 of file FaceOffset_3.h.

References reset_default_parameters().

                  : _is_strtd(0), _As(NULL), _boffset(NULL), _bmedial(NULL), 
                    _eigvalues(NULL), _eigvecs(NULL), _vnormals(NULL),
                    _facenormals(NULL), _facecenters(NULL), 
                    _faceareas(NULL), _vcenters(NULL), _tangranks(NULL), 
                    _ctangranks(NULL), _scales(NULL), _weights(NULL)
   { reset_default_parameters(); }

Here is the call graph for this function:

FaceOffset_3	(	Manifold *	wm,
		COM::Window *	buf
	)

Construct an object from a window manifold and a buffer window.

Definition at line 81 of file FaceOffset_3.C.

References _As, _bmedial, _boffset, Propagation_3::_buf, _ctangranks, _eigvalues, _eigvecs, _faceareas, _facecenters, _facenormals, _is_strtd, Propagation_3::_panes, _ridgeneighbors, _ridges, _scales, _strong, _tangranks, _vcenters, _vnormals, _weak, _weights, COM_CHAR, COM_DOUBLE, COM_INT, COMMPI_Initialized(), i, MPI_MIN, reset_default_parameters(), and s.

   : Propagation_3( wm, buf)
 {
   // Determine whether the mesh is structured
   _is_strtd = true;
   std::vector< COM::Pane*>::iterator it = _panes.begin();
   for (int i=0, local_npanes = _panes.size(); i<local_npanes; ++i, ++it) { 
     COM::Pane *pane = *it;
     if ( !pane->is_structured()) _is_strtd = false;
   }
   if ( COMMPI_Initialized()) {
     int s = _is_strtd;
     MPI_Allreduce( &s, &_is_strtd, 1, MPI_INT, MPI_MIN, 
                    _buf->get_communicator());
   }
 
   // Set parameters to default values. 
   reset_default_parameters( _is_strtd);
 
   // Extend buffer window
   _facenormals = _buf->new_attribute( "facenormals", 'e', COM_DOUBLE, 3, "");
   _buf->resize_array( _facenormals, 0);
 
   _facecenters = _buf->new_attribute( "facecenters", 'e', COM_DOUBLE, 3, "");
   _buf->resize_array( _facecenters, 0);
 
   _faceareas = _buf->new_attribute( "faceareas", 'e', COM_DOUBLE, 1, "");
   _buf->resize_array( _faceareas, 0);
 
   // This is used to store v-center for smoothing offset.
   _vcenters = _buf->new_attribute( "vcenters", 'n', COM_DOUBLE, 3, "");
   _buf->resize_array( _vcenters, 0);
 
   _eigvalues = _buf->new_attribute( "lambda", 'n', COM_DOUBLE, 3, "");
   _buf->resize_array( _eigvalues, 0);
 
   _vnormals = _buf->new_attribute( "vnormals", 'n', COM_DOUBLE, 3, "");
   _buf->resize_array( _vnormals, 0);
 
   _ridges = _buf->new_attribute( "ridges", 'p', COM_INT, 2, "");
   
   _ridgeneighbors = _buf->new_attribute( "ridgeneighbors", 'n', COM_INT, 4, "");
   _buf->resize_array( _ridgeneighbors, 0);
 
   _weak = _buf->new_attribute( "weakfeature", 'n', COM_CHAR, 1, "");
   _buf->resize_array( _weak, 0);
 
   _strong = _buf->new_attribute( "strongfeature", 'n', COM_CHAR, 1, "");
   _buf->resize_array( _strong, 0);
 
   _As = _buf->new_attribute( "As", 'n', COM_DOUBLE, 9, "");
   _buf->resize_array( _As, 0);
 
   _boffset = _buf->new_attribute( "boffset", 'n', COM_DOUBLE, 3, "");
   _buf->resize_array( _boffset, 0);
 
   _bmedial = _buf->new_attribute( "bmedial", 'n', COM_DOUBLE, 3, "");
   _buf->resize_array( _bmedial, 0);
 
   _eigvecs = _buf->new_attribute( "eigvecs", 'n', COM_DOUBLE, 9, "");
   _buf->resize_array( _eigvecs, 0);
 
   _tangranks = _buf->new_attribute( "tangranks", 'n', COM_CHAR, 1, "");
   _buf->resize_array( _tangranks, 0);
 
   _ctangranks = _buf->new_attribute( "ctangranks", 'n', COM_CHAR, 1, "");
   _buf->resize_array( _ctangranks, 0);
 
   _scales = _buf->new_attribute( "scales", 'n', COM_DOUBLE, 1, "");
   _buf->resize_array( _scales, 0);
 
   _weights = _buf->new_attribute( "weights", 'n', COM_DOUBLE, 1, "");
   _buf->resize_array( _weights, 0);
 
   _buf->init_done(false);
 }

Here is the call graph for this function:

virtual ~FaceOffset_3 ( )

inlinevirtual

Definition at line 53 of file FaceOffset_3.h.

53 {}

Member Function Documentation

void adjust_wavefrontal_motion ( COM::Attribute * disps )

protected

int append_obscure_ends	(	const COM::Attribute *	maxtrnangv,
		const COM::Attribute *	mintrnangv,
		const std::vector< std::map< Edge_ID, double > > &	edge_maps,
		const COM::Attribute *	maxranglev,
		const COM::Attribute *	minranglev,
		const std::vector< std::map< Edge_ID, double > > &	rangle_maps
	)

protected

PROP_BEGIN_NAMESPACE void balance_mass ( )

protected

Definition at line 33 of file cons_diff.C.

References Propagation_3::_buf, _faceareas, _facenormals, Propagation_3::_panes, _scales, Propagation_3::_surf, _tangranks, _vcenters, _vnormals, _weights, Rocblas::add(), COM_DOUBLE, COM_NC, Rocblas::copy(), Rocblas::copy_scalar(), distribute_volume_e2n(), Rocblas::dot(), E2N_ONE, i, j, k, and Rocblas::sub().

Referenced by obtain_constrained_directions().

                {
   COM::Attribute *disps = 
     _buf->new_attribute( "disps_buf_bm", 'n', COM_DOUBLE, 3, "");
   _buf->resize_array( disps, 0);
   COM::Attribute *pos = 
     _buf->new_attribute( "pos_buf_bm", 'n', COM_DOUBLE, 3, "");
   _buf->resize_array( pos, 0);
   COM::Attribute *fvol_buf = 
     _buf->new_attribute( "vfol_buf_bm", 'e', COM_DOUBLE, 1, "");
   _buf->resize_array( fvol_buf, 0);
   _buf->init_done( false);
 
   // Use face-area as buffer for facial-volume
   COM::Attribute *fvol = _faceareas;
 
   // Compute volume change at each vertex.
   const COM::Attribute *coors = _buf->attribute( COM::COM_NC);
   Rocblas::add( coors, _vcenters, pos);
 
   // Obtain density-vector for vertices.
   int normalize=false;
   SURF::Rocsurf::compute_element_normals( _facenormals, &normalize, pos);
   _surf->elements_to_nodes( _facenormals, disps, SURF::E2N_ONE, NULL, NULL, 1);
   
   // Use _scales for vortex area and _weights for vertex-volume.
   Rocblas::dot( disps, _vnormals, _scales);
   double zero=0;
   Rocblas::copy_scalar( &zero, disps);
 
   for ( int k=0;k<2;++k) {
     if ( k) {
       SURF::Rocsurf::compute_swept_volumes( pos, disps, fvol);
 
       Rocblas::sub( fvol_buf, fvol, fvol);
     }
     else {
       SURF::Rocsurf::compute_bounded_volumes( pos, coors, fvol);
       Rocblas::copy( fvol, fvol_buf);
     }
     distribute_volume_e2n( fvol, _tangranks, _weights);
 
 
     std::vector< COM::Pane*>::iterator it = _panes.begin();
     // Loop through the panes and its real vertices
     for (int i=0, local_npanes = _panes.size(); i<local_npanes; ++i, ++it) { 
       COM::Pane *pane = *it;
       
       const double *vs = reinterpret_cast<const double*>
         ( pane->attribute(_weights->id())->pointer());
       const double *as = reinterpret_cast<const double*>
         ( pane->attribute(_scales->id())->pointer());
       const Vector_3 *ds_m = reinterpret_cast<const Vector_3*>
         ( pane->attribute(_vnormals->id())->pointer());
       Vector_3 *ds = reinterpret_cast<Vector_3*>
         ( pane->attribute(disps->id())->pointer());
 
       // Loop through all real nodes of the pane
       for ( int j=0, jn=pane->size_of_real_nodes(); j<jn; ++j) {
         ds[j] += vs[j]/as[j]*ds_m[j]; // Omit the factor 3 here
       }
     }
   }
 
   Rocblas::add( _vcenters, disps, _vcenters);
   
   _buf->delete_attribute( fvol_buf->name());
   _buf->delete_attribute( pos->name());
   _buf->delete_attribute( disps->name());
   _buf->init_done( false);
 }

Here is the call graph for this function:

Here is the caller graph for this function:

int build_ridge_neighbor_list	(	const COM::Attribute *	maxtrnangv,
		const COM::Attribute *	mintrnangv,
		const std::vector< std::map< Edge_ID, double > > &	edge_maps,
		const COM::Attribute *	maxranglev,
		const COM::Attribute *	minranglev,
		const std::vector< std::map< Edge_ID, double > > &	rangle_maps,
		std::vector< ObscendSet > &	obscends
	)

protected

Build the list of ridge nighbors and obtain a list of weak-ended vertices Return the number of obscended vertices.

Definition at line 570 of file detect_features.C.

References _edges, Propagation_3::_panes, _ridgeneighbors, _strong, Propagation_3::_surf, _tangranks, _tol_angle_strong, _tol_cos_osta, _tol_eangle, _tol_turn_angle, NTS::abs(), COM_assertion, COM_NC, Rocblas::copy_scalar(), eval_angle(), Element_node_enumerator::size_of_edges(), and RidgeNeighbor::vid.

Referenced by filter_and_identify_ridge_edges().

                                                               {
 
   COM_assertion( _panes.size()==1); // Precondition
   
   // Loop through the panes and ridge edges
   Manifold::PM_iterator pm_it=_surf->pm_begin();
   COM::Pane *pane = _panes[0];
 
   // Build ACH list
   typedef std::map< int, std::list< RidgeNeighbor> >  ACH_Lists;
   ACH_Lists ach; // ACH list
 
   const std::set< Edge_ID> &eset_pn = _edges[0];  // List of ridge edges
   for ( std::set< Edge_ID>::const_iterator eit=eset_pn.begin(),
           eend=eset_pn.end(); eit!=eend; ++eit) {
     // Make sure that eid_opp is not border
     Edge_ID eid = *eit, eid_opp = pm_it->get_opposite_real_edge( eid);
     COM_assertion( eset_pn.find(eid_opp)!=eset_pn.end() &&
                    !eid.is_border());
     Element_node_enumerator ene( pane, eid.eid()); 
     int lid=eid.lid(), neighbor = (lid+1)%ene.size_of_edges();
 
     // Fill in ridge-neighbor into _ridgeneighbors.
     int vid = ene[lid], neighbor_id=ene[neighbor];
     ach[vid].push_back( RidgeNeighbor( neighbor_id, eid));
   }
 
   // initialize storage
   int zero=0;
   Rocblas::copy_scalar( &zero, _ridgeneighbors);
   obscends.clear(); obscends.resize( _panes.size());
 
   // Check ACH list and build rdgngbs
   RidgeNeighbor *rdgngbs = reinterpret_cast<RidgeNeighbor*>
     ( pane->attribute(_ridgeneighbors->id())->pointer());
   const Point_3 *pnts = reinterpret_cast<Point_3*>
     (pane->attribute(COM_NC)->pointer());
   const char *tranks = reinterpret_cast<char*>
     ( pane->attribute(_tangranks->id())->pointer());
   const double *minvs = reinterpret_cast<double*>
     ( pane->attribute(mindianglev->id())->pointer());
   const double *maxvs = reinterpret_cast<double*>
     ( pane->attribute(maxdianglev->id())->pointer());
   const double *minta = reinterpret_cast<double*>
     (pane->attribute(mintrnangv->id())->pointer());
   const double *maxta = reinterpret_cast<double*>
     ( pane->attribute(maxtrnangv->id())->pointer());
   char *strong = reinterpret_cast<char*>
     ( pane->attribute(_strong->id())->pointer());
   const std::map< Edge_ID, double> &diangle_pn = diangle_maps[0];
   const std::map< Edge_ID, double> &emap_pn = edge_maps[0];
   ObscendSet &obscend_pn = obscends[0];
 
   for ( ACH_Lists::iterator ach_it = ach.begin(); ach_it!=ach.end(); ++ach_it) {
     // Check angle defect and turning angle
     int vid= ach_it->first, index = (vid-1)*2;
     int nedges = ach_it->second.size();
     std::list< RidgeNeighbor>::iterator ait = ach_it->second.begin();
 
     if (nedges>2) {
       // Insert into rdgngbs
       for (; ait!=ach_it->second.end(); ) {
         // Remove non-joint halfedges
         Edge_ID eid = ait->eid;
         double fangle = diangle_pn.find( eid)->second;
         double feg = emap_pn.find( eid)->second;
 
         // If disjoint, then add edge into obscend_pn.
         if ( fangle < _tol_eangle &&
              ( strong[vid-1] || tranks[vid-1]==1 && 
                ( (-feg < _tol_cos_osta || feg != minta[vid-1]) &&
                  ( feg < _tol_cos_osta || feg != maxta[vid-1]) ||
                  ( feg<0 || fangle != maxvs[vid-1]) && 
                  ( feg>0 || fangle != -minvs[vid-1]))) ||
              fangle < _tol_angle_strong && 
              ( strong[vid-1]==2 || tranks[vid-1]==1 &&
                ( std::abs(feg) < _tol_cos_osta ||
                  feg != minta[vid-1] && feg != maxta[vid-1]) && 
                ( feg<0 || fangle != maxvs[vid-1]) && 
                ( feg>0 || fangle != -minvs[vid-1]))) {
           obscend_pn[ ait->eid] = std::make_pair( vid, ait->vid);
           ait = ach_it->second.erase( ait);
         }
         else
           ++ait;
       }
 
       rdgngbs[index].vid = 0; rdgngbs[index+1].vid = -1; 
     }
 
     // Classify dangling and semi-joint halfedges
     ait = ach_it->second.begin();
     if ( nedges == 1) {
       // Insert into obscure list to indicate dangling edges
       rdgngbs[index] = *ait; rdgngbs[index+1].vid=0;
       Edge_ID eid = ait->eid;
       double fangle = diangle_pn.find( eid)->second;
 
       if ( fangle < _tol_angle_strong || tranks[vid-1])
         obscend_pn[ eid] = std::make_pair( ach_it->first, ait->vid);
     }
     else if ( nedges == 2) {
       std::list< RidgeNeighbor>::iterator ait2 = ait; ++ait2;
       int v1 = ait->vid, v2 = ait2->vid;
       double fangle1 = diangle_pn.find( ait->eid)->second;
       double fangle2 = diangle_pn.find( ait2->eid)->second;
 
       if ( (fangle1<_tol_angle_strong || fangle2<_tol_angle_strong) &&
            ( !tranks[vid-1] || eval_angle
              ( pnts[v1-1]-pnts[vid-1], pnts[vid-1]-pnts[v2-1]) > _tol_turn_angle)) {
         // Do not insert the edges into rdgngbs but insert into obscend_pn
         obscend_pn[ ait->eid] = std::make_pair( vid, v1);
         obscend_pn[ ait2->eid] = std::make_pair( vid, v2);
       }
       else if (!rdgngbs[index+1].vid)
       { rdgngbs[index] = *ait; rdgngbs[index+1] = *ait2; }
     }
   }
     
   int  n_obscended = obscends[0].size();
   std::cerr << "There are " << n_obscended 
             << " obscure-ended edges" << std::endl;
   return n_obscended;
 }

Here is the call graph for this function:

Here is the caller graph for this function:

std::pair<double,double> comp_wavefrontal_motion	(	double	det,
		double	w,
		double	delta,
		const Vector_3 &	disp_nz,
		const Vector_3 &	nrm_nz
	)

protected

void compute_anisotropic_vertex_centers ( const COM::Attribute * disps_buf )

protected

Definition at line 85 of file AnisotropicSmoothing.C.

References _As, _bmedial, _facecenters, _feature_layer, Propagation_3::_panes, _strong, Propagation_3::_surf, _tangranks, _vcenters, _weights, COM_NC, Rocblas::copy_scalar(), Rocblas::div(), eigen_analyze_vertex(), eval_angle(), get_tangents(), i, j, k, min(), Element_node_enumerator::next(), nj, and Element_node_enumerator::size_of_edges().

Referenced by time_stepping().

                                                                      {
 
   const double zero = 0., eps = 1.e-100;
 
   Rocblas::copy_scalar( &zero, _vcenters);
   Rocblas::copy_scalar( &eps, _weights);
   
   // Loop through the panes and its real faces
   std::vector< COM::Pane*>::iterator it = _panes.begin();
   Manifold::PM_iterator pm_it=_surf->pm_begin();
   for ( int i=0, local_npanes = _panes.size(); 
         i<local_npanes; ++i, ++it, ++pm_it) { 
     COM::Pane *pane = *it;
 
     const Point_3 *pnts = reinterpret_cast<const Point_3*>
       (pane->attribute(COM_NC)->pointer());
     const Vector_3 *disps = nodal_disps?reinterpret_cast<const Vector_3*>
       (pane->attribute(nodal_disps->id())->pointer()):NULL;
     const Point_3 *fcnts = reinterpret_cast<const Point_3*>
       ( pane->attribute(_facecenters->id())->pointer());
     Vector_3 *vcnts = reinterpret_cast<Vector_3*>
       ( pane->attribute(_vcenters->id())->pointer());
     double   *ws = reinterpret_cast<double*>
       ( pane->attribute(_weights->id())->pointer());
 
 #if ANISOTROPIC    
     const Vector_3 *As = reinterpret_cast<const Vector_3*>
       ( pane->attribute(_As->id())->pointer());
     const Vector_3 *bs = reinterpret_cast<const Vector_3*>
       ( pane->attribute(_bmedial->id())->pointer());
 
     const char *tranks = reinterpret_cast<const char*>
       ( pane->attribute(_tangranks->id())->pointer());
 
     const char *strong = reinterpret_cast<char*>
       ( pane->attribute(_strong->id())->pointer());
 #endif
 
     // Loop through real elements of the current pane
     Element_node_enumerator ene( pane, 1);
     for ( int j=0, nj=pane->size_of_real_elements(); j<nj; ++j, ene.next()) {
       int ne = ene.size_of_edges();
       const Point_3  &cnt = fcnts[j];
       Vector_3 anipnt;
 
 #if ANISOTROPIC
       bool has_feature=false;
       Vector_3 Ae[3], mn;   // Matrix at an edge
 
       // First, compute medial quadric for the face
       Ae[0] = Ae[1] = Ae[2] = mn = Vector_3(0,0,0);
 
       for ( int k=0; k<ne; ++k) {
         int vindex = ene[k]-1;
         has_feature |= tranks[vindex]!=2;
         
         // Obtain the quadric at the face 
         const Vector_3 *As_v=&As[3*vindex];
         Ae[0] += As_v[0];
         Ae[1] += As_v[1];
         Ae[2] += As_v[2];
             
         mn += bs[vindex];
       }
 
       // Perform eigen-analysis
       Vector_3 lambdas = mn;
       int nrank = eigen_analyze_vertex( Ae, lambdas);
       
       Vector_3 vs_e[2];
       double   scales_e[2];
 
       mn = Ae[0]; nrank = 1;
       get_tangents( Ae, lambdas, mn, nrank, vs_e, scales_e);
 
       Vector_3 mn_face;
       Vector_3 vs_f[2];
       double   scales_f[2];
       if ( has_feature) {
         // Compute eigen-vector of weighted-sum of the tensors.
         Ae[0] = Ae[1] = Ae[2] = Vector_3(0,0,0);
         for ( int k=0; k<ne; ++k) {
           int vindex = ene[k]-1;
 
           double w=1;
           switch (tranks[ vindex]) {
           case 0: {
             w = 1.e-12;
             break;
           }
           case 1: {
             if ( strong[vindex]) { 
               w = 1.e-6;
               break;
             }
           }
           default: w = 1.;
           }
           
           // Obtain the quadric at the face 
           const Vector_3 *As_v=&As[3*vindex];
           Ae[0] += w*As_v[0];
           Ae[1] += w*As_v[1];
           Ae[2] += w*As_v[2];
         }
         
         // Perform eigen-analysis
         eigen_analyze_vertex( Ae, lambdas);
         mn_face = Ae[0];
         
         get_tangents( Ae, lambdas, mn_face, 1, vs_f, scales_f);
       }
       else {
         mn_face = Ae[0];
 
         scales_f[0] = scales_e[0]; scales_f[1] = scales_e[1]; 
         vs_f[0] = vs_e[0]; vs_f[1] = vs_e[1]; 
       }
 #endif
 
       int uindex=ene[ne-1]-1, vindex=ene[0]-1;
       // Loop through all vertices of current face.
       for ( int k=0; k<ne; ++k) {
         int windex = ene[(k+1==ne)?0:k+1]-1;
 
         Point_3 ps[3] = { pnts[uindex], pnts[vindex], pnts[windex] };
         if ( disps) { 
           ps[0] += disps[uindex]; ps[1] += disps[vindex]; 
           ps[2] += disps[windex]; 
         }
 
 #if ANISOTROPIC 
         double scales_v[2];
         Vector_3 vs_v[2];
 
         // Note: is_border node.
         if ( _feature_layer && 
              ( tranks[uindex]==2 || tranks[windex]==2)) {
           scales_v[0] = scales_e[0]; scales_v[1] = scales_e[1]; 
           vs_v[0] = vs_e[0]; vs_v[1] = vs_e[1];
         }
         else  {
           scales_v[0] = scales_f[0]; scales_v[1] = scales_f[1]; 
           vs_v[0] = vs_f[0]; vs_v[1] = vs_f[1];
         }
         
         const Vector_3 diff = cnt-ps[1];
         anipnt = scales_v[0]*(diff*vs_v[0])*vs_v[0]+
           scales_v[1]*(diff*vs_v[1])*vs_v[1];
         if ( disps) anipnt += disps[vindex];
 #else
         anipnt = cnt-ps[1];
 #endif
 
 #if MIN_ANGLE_WEIGHTED
         double w = std::min( eval_angle( ps[1]-ps[0], ps[2]-ps[0]),
                              eval_angle( ps[1]-ps[2], ps[0]-ps[2]));
 #else
         double  w=1;
 #endif
 
         vcnts[vindex] += w*anipnt;
         ws[vindex] += w;
 
         uindex=vindex; vindex=windex;
       }
     }
   }
 
   _surf->reduce_on_shared_nodes( _vcenters, Manifold::OP_SUM);
   _surf->reduce_on_shared_nodes( _weights, Manifold::OP_SUM);
   Rocblas::div( _vcenters, _weights, _vcenters);
 }

Here is the call graph for this function:

Here is the caller graph for this function:

void compute_directions	(	COM::Attribute *	b_offset,
		bool	in_principle
	)

protected

Compute mean normals.

Definition at line 447 of file FaceOffset_3.C.

References _eigvalues, _eigvecs, Propagation_3::_panes, _tangranks, _vnormals, i, j, and obtain_directions().

Referenced by time_stepping().

                                                                         {
   // Loop through the panes and its real vertices
   std::vector< COM::Pane*>::iterator it = _panes.begin();
   for (int i=0, local_npanes = _panes.size(); i<local_npanes; ++i, ++it) { 
     COM::Pane *pane = *it;
 
     // Obtain pointers
     Vector_3 *evs = reinterpret_cast<Vector_3*>
       ( pane->attribute(_eigvecs->id())->pointer());
     const Vector_3 *lambdas = reinterpret_cast<const Vector_3*>
       ( pane->attribute(_eigvalues->id())->pointer());
     Vector_3 *vnrms = reinterpret_cast<Vector_3*>
       ( pane->attribute(_vnormals->id())->pointer());
     Vector_3 *voffset = bo_attr ? reinterpret_cast<Vector_3*>
       ( pane->attribute(bo_attr->id())->pointer()) : NULL;
     const char *trank = princ ? reinterpret_cast<char*>
       ( pane->attribute(_tangranks->id())->pointer()) : NULL;
 
     // Loop through all real nodes of the pane
     for ( int j=0, jn=pane->size_of_real_nodes(); j<jn; ++j) {
       obtain_directions( &evs[3*j], lambdas[j], princ ? 3-trank[j] : 3, 
                          vnrms[j], voffset?voffset+j:NULL);
     }
   }
 }

Here is the call graph for this function:

Here is the caller graph for this function:

void compute_eigenvectors	(	Vector_3	A[3],
		Vector_3 &	lambdas
	)

staticprotected

Definition at line 572 of file quadric_analysis.C.

References COM_assertion_msg, dsyevq3(), cimg_library::cimg::info(), and swap().

Referenced by denoise_surface(), and eigen_analyze_vertex().

                                                         {
 
   double abuf[3][3] = { {A[0][0], A[0][1], A[0][2]},
                         {A[1][0], A[1][1], A[1][2]},
                         {A[2][0], A[2][1], A[2][2]}};
   double ebuf[3][3];
   
   int info = dsyevq3( abuf, ebuf, &lambdas[0]);
   COM_assertion_msg( info==0, "Computation of eigenvectos failed");
 
   std::swap( ebuf[0][1], ebuf[1][0]);
   std::swap( ebuf[0][2], ebuf[2][0]);
   std::swap( ebuf[1][2], ebuf[2][1]);
 
   const Vector_3 *buf = (const Vector_3*)&ebuf[0][0];
 
   if ( lambdas[0]<lambdas[1]) {
     if ( lambdas[1]<lambdas[2]) { // l2>l1>l0
       std::swap( lambdas[0], lambdas[2]);
       A[0] = buf[2]; A[1] = buf[1]; A[2] = buf[0];
     }
     else { 
       double t = lambdas[0];
       lambdas[0] = lambdas[1]; 
       A[0] = buf[1];
 
       if ( t<lambdas[2]) { // l1>=l2>l0
         lambdas[1] = lambdas[2]; lambdas[2] = t;
         A[1] = buf[2]; A[2] = buf[0];
       }
       else { // l1>l0>=l2
         lambdas[1] = t;
         A[1] = buf[0]; A[2] = buf[2];
       }
     }
   }
   else if ( lambdas[0]<lambdas[2]) { // l2>l0>=l1
     double t = lambdas[0];
     lambdas[0] = lambdas[2]; lambdas[2] = lambdas[1]; lambdas[1] = t;
     A[0] = buf[2]; A[1] = buf[0]; A[2] = buf[1];
   }
   else {
     A[0] = buf[0];
     if ( lambdas[1]<lambdas[2]) { // l0>=l2>l1
       std::swap( lambdas[1], lambdas[2]);
       A[1] = buf[2]; A[2] = buf[1];
     }
     else {// l0>=l1>=l2
       A[1] = buf[1]; A[2] = buf[2]; 
     }
   }
 }

Here is the call graph for this function:

Here is the caller graph for this function:

PROP_BEGIN_NAMESPACE void compute_quadrics	(	double	dt,
		const COM::Attribute *	spds,
		COM::Attribute *	rhs,
		COM::Attribute *	predicted
	)

protected

Compute quadrics for every vertex.

Definition at line 34 of file quadric_analysis.C.

References _As, _bmedial, Propagation_3::_buf, _dir_thres_weak, _eigvalues, _eigvecs, _faceareas, _facecenters, _facenormals, _is_strtd, Propagation_3::_panes, _smoother, Propagation_3::_surf, _tangranks, _vcenters, _vnormals, _weights, _wght_scheme, angle(), COM_assertion_msg, COM_DOUBLE, COM_NC, copy, Rocblas::copy(), Rocblas::copy_scalar(), Rocblas::div(), E2N_ANGLE, E2N_AREA, E2N_ONE, eigen_analyze_vertex(), eval_angle(), i, j, k, Element_node_enumerator::next(), nj, obtain_face_offset(), pi, s, Element_node_enumerator::size_of_edges(), SMOOTHER_ANISOTROPIC, SMOOTHER_LAPLACIAN, and sqrt().

Referenced by time_stepping().

                                                                   {
 
   COM_assertion_msg( spds || dt==0, 
                      "If no speed function, then time step must be zero.");
 
   const double zero = 0., eps = 1.e-100;
   if ( bo_attr) Rocblas::copy_scalar( &zero, bo_attr);
   if ( dt==0)   bo_attr = NULL;
 
   // Use the following buffer spaces: _As for matrix A; _eigvalues to 
   // store bm; _sumangles to store sum of angles
   COM::Attribute *A_attr = _As;
   COM::Attribute *bm_attr = _eigvalues;
   COM::Attribute *w_attr = _weights;
 
   Rocblas::copy_scalar( &zero, A_attr);
   Rocblas::copy_scalar( &zero, bm_attr);
   Rocblas::copy_scalar( &eps, w_attr);
   Rocblas::copy_scalar( &zero, _vcenters);
 
   COM::Attribute *sa_attr = 
     _buf->new_attribute( "sumangles", 'n', COM_DOUBLE, 1, "");
   _buf->resize_array( sa_attr, 0);
   
   bool shear= _smoother==SMOOTHER_LAPLACIAN || _is_strtd ||
     _smoother!=SMOOTHER_ANISOTROPIC && (dt==0 || bo_attr==0);
 
   // Loop through the panes and its real faces
   std::vector< COM::Pane*>::iterator it = _panes.begin();
   Manifold::PM_iterator pm_it=_surf->pm_begin();
   for ( int i=0, local_npanes = _panes.size(); 
         i<local_npanes; ++i, ++it, ++pm_it) { 
     COM::Pane *pane = *it;
 
     // Obtain nodal coordinates of current pane, assuming contiguous layout
     const Point_3 *pnts = reinterpret_cast<Point_3*>
       (pane->attribute(COM_NC)->pointer());
     Vector_3 *fnrms = reinterpret_cast<Vector_3*>
       ( pane->attribute(_facenormals->id())->pointer());
     Point_3 *fcnts = reinterpret_cast<Point_3*>
       ( pane->attribute(_facecenters->id())->pointer());
 
     double   *ws = reinterpret_cast<double*>
       ( pane->attribute(w_attr->id())->pointer());
     Vector_3 *vcnts = reinterpret_cast<Vector_3*>
       ( pane->attribute(_vcenters->id())->pointer());
 
     // Obtain the pointer to speed function
     const double *spds_ptr = spds?reinterpret_cast<const double*>
       (pane->attribute( spds->id())->pointer()):NULL;
 
     // Obtain pointers
     Vector_3 *As = reinterpret_cast<Vector_3*>
       ( pane->attribute(A_attr->id())->pointer());
     Vector_3 *bs_m = reinterpret_cast<Vector_3*>
       ( pane->attribute(bm_attr->id())->pointer());
     Vector_3 *bs_o = bo_attr ? reinterpret_cast<Vector_3*>
       ( pane->attribute(bo_attr->id())->pointer()) : NULL;
     Vector_3 *pred_dirs = predicted ? reinterpret_cast<Vector_3*>
       ( pane->attribute(predicted->id())->pointer()) : NULL;
     double   *sa = reinterpret_cast<double*>
       ( pane->attribute(sa_attr->id())->pointer());
     double   *areas = reinterpret_cast<double*>
       ( pane->attribute(_faceareas->id())->pointer());
 
     // Loop through real elements of the current pane
     Element_node_enumerator ene( pane, 1); 
     for ( int j=0, nj=pane->size_of_real_elements(); j<nj; ++j, ene.next()) {
       // If speed is uniform within a face, then unit normal is the same
       // after propagation.
       Point_3  &cnt = fcnts[j];
 
       int ne = ene.size_of_edges();
       int uindex=ene[ne-1]-1, vindex=ene[0]-1, windex=0;
 
       Vector_3 disps[9], nrms[9];
 
       // Compute signed distance from current face.
       obtain_face_offset( pnts, spds_ptr, spds, dt, pred_dirs,
                           ene, fnrms[j], cnt, disps, nrms);
 
       double   delta=0;
       // Loop through all vertices of current face.
       for ( int k=0; k<ne; ++k, uindex=vindex, vindex=windex) {
         windex = ene[(k+1==ne)?0:k+1]-1;
 
         // Evaluate angle at each vertex. 
         Vector_3 ts[2] = { pnts[uindex]-pnts[vindex], 
                            pnts[windex]-pnts[vindex]};
         Vector_3 ns_nz = nrms[k];
 
         if ( bo_attr) { // Update delta
           delta = -disps[k]*ns_nz;
         }
 
         double angle = eval_angle( ts[0], ts[1]);
         Vector_3 nrm = ns_nz;
         double w=1;
         if ( _wght_scheme!=SURF::E2N_ONE) {
           switch ( _wght_scheme) {
           case SURF::E2N_ANGLE: {
             double s = std::sqrt((ts[0]*ts[0])*(ts[1]*ts[1]));
             if ( s==0) w = 0; 
             else w = angle; 
 
             break;
           }
           case SURF::E2N_AREA: {
             w = areas[j]; break;
           }
           default: ; // w=1;
           }
           nrm *= w;
         }
 
         // Obtain the addresses of the linear system and weight
         Vector_3 *A_v = &As[3*vindex];
 
         // Update A, b_o, b_m, and as_v for the current vertex
         A_v[0] += ns_nz[0]*nrm;
         A_v[1] += ns_nz[1]*nrm;
         A_v[2] += ns_nz[2]*nrm;
 
         if ( delta!=0) bs_o[vindex] += delta*nrm;
         bs_m[vindex] -= nrm;
         sa[vindex] += angle;
 
         if ( shear) {
           Vector_3 diff = cnt-pnts[vindex];
           ws[vindex] += 1.;
           vcnts[vindex] += diff;
         }
       } // for k (edges)
     } // for j (elements)
   } // for i (panes)
   
   // Reduce on shared nodes for A, b_m, b_o, and r
   _surf->reduce_on_shared_nodes( A_attr, Manifold::OP_SUM);
   _surf->reduce_on_shared_nodes( bm_attr, Manifold::OP_SUM);
   if ( bo_attr) _surf->reduce_on_shared_nodes( bo_attr, Manifold::OP_SUM);
   _surf->reduce_on_shared_nodes( sa_attr, Manifold::OP_SUM);
 
   // Copy bm_attr into _vnormals and _bmedial
   Rocblas::copy( bm_attr, _vnormals);
   Rocblas::copy( bm_attr, _bmedial);
 
   // Reduce to maxabs so that shared nodes have exactly the same value
   _surf->reduce_on_shared_nodes( A_attr, Manifold::OP_MAXABS);
   _surf->reduce_on_shared_nodes( sa_attr, Manifold::OP_MAXABS);
 
   if ( shear) {
     _surf->reduce_on_shared_nodes( w_attr, Manifold::OP_SUM);
     _surf->reduce_on_shared_nodes( _vcenters, Manifold::OP_SUM);
     // Obtain angle-weighted average of mass center at each vertex.
     Rocblas::div( _vcenters, w_attr, _vcenters);
   }
 
   // Loop through the panes and its real vertices to determine trank
   it = _panes.begin(); 
   for (int i=0, local_npanes = _panes.size(); i<local_npanes; ++i, ++it) { 
     COM::Pane *pane = *it;
 
     // Obtain pointers
     Vector_3 *As = reinterpret_cast<Vector_3*>
       ( pane->attribute(A_attr->id())->pointer());
     Vector_3 *evs = reinterpret_cast<Vector_3*>
       ( pane->attribute(_eigvecs->id())->pointer());
     Vector_3 *bs_m = reinterpret_cast<Vector_3*>
       ( pane->attribute(bm_attr->id())->pointer());
     double   *sa = reinterpret_cast<double*>
       ( pane->attribute(sa_attr->id())->pointer());
 
     // Feature information
     char      *tranks = reinterpret_cast<char*>
       ( pane->attribute(_tangranks->id())->pointer());
 
     // Loop through all real nodes of the pane
     for ( int j=0, jn=pane->size_of_real_nodes(); j<jn; ++j) {
       // Skip vertices with zero right-hand side, such as those at
       // edge- and face-centers for quadratic elements
       if ( bs_m[j].is_null()) continue;
 
       Vector_3 *A_v = &As[3*j];
       Vector_3 *es = &evs[3*j];    // eigen-vectors of current vertex.
       std::copy( A_v, A_v+3, es);
       
       // Perform eigen-analysis for the vertex using medial quadric
       // nrank is the rank of the normal (primary) space of A.
       char nrank = eigen_analyze_vertex( es, bs_m[j], 2*pi-sa[j]);
 
       // Disable features if f-angle is greater than 1.
       if ( _dir_thres_weak>0.99) nrank=1;
       tranks[j] = 3-nrank;
     }
   }
 
   _buf->delete_attribute( sa_attr->name());
   _buf->init_done( false);
 
   // Reduce tangranks to ensure shared nodes have the same value across panes
   _surf->reduce_on_shared_nodes( _tangranks, Manifold::OP_MIN);
 }

Here is the call graph for this function:

Here is the caller graph for this function:

void compute_volume_vector	(	const COM::Attribute *	disps,
		COM::Attribute *	bs
	)

protected

Compute volumn-vector.

Vector_3 compute_weighted_normal	(	Vector_3	es[3],
		const Vector_3 &	lambdas,
		int	trank,
		const Vector_3 &	mean_nrm,
		const Vector_3 &	face_nrm
	)

protected

void denoise_surface	(	const COM::Attribute *	disps_buf,
		COM::Attribute *	normal_motion
	)

protected

Definition at line 260 of file AnisotropicSmoothing.C.

References _As, Propagation_3::_cnstr_nodes, _eigvalues, _faceareas, _facecenters, Propagation_3::_panes, Propagation_3::_surf, _tangranks, _vnormals, _weak, _weights, COM_NC, compute_eigenvectors(), Rocblas::copy_scalar(), Rocblas::div(), get_nonuniform_weight(), i, j, k, nj, and offset().

Referenced by time_stepping().

                                               {
   
   const double zero = 0., eps = 1.e-100;
 
   Rocblas::copy_scalar( &zero, normal_motion);
   Rocblas::copy_scalar( &eps, _weights);
 
   // Loop through the panes and its real faces
   std::vector< COM::Pane*>::iterator it = _panes.begin();
   Manifold::PM_iterator pm_it=_surf->pm_begin();
   for ( int i=0, local_npanes = _panes.size(); 
         i<local_npanes; ++i, ++it, ++pm_it) { 
     COM::Pane *pane = *it;
     
     const Vector_3 *As = reinterpret_cast<const Vector_3*>
       ( pane->attribute(_As->id())->pointer());
     const Vector_3 *vnrms = reinterpret_cast<const Vector_3*>
       ( pane->attribute(_vnormals->id())->pointer());
 
     const Point_3 *pnts = reinterpret_cast<const Point_3*>
       (pane->attribute(COM_NC)->pointer());
     const Vector_3 *disps = nodal_disps ? reinterpret_cast<const Vector_3*>
       (pane->attribute(nodal_disps->id())->pointer()) : NULL;
     const Point_3 *fcnts = reinterpret_cast<const Point_3*>
       ( pane->attribute(_facecenters->id())->pointer());
     const char *tranks = reinterpret_cast<const char*>
       ( pane->attribute(_tangranks->id())->pointer());
     const Vector_3 *lambdas = reinterpret_cast<const Vector_3*>
       ( pane->attribute(_eigvalues->id())->pointer());
 
     Vector_3 *nmotion = normal_motion ? reinterpret_cast<Vector_3*>
       ( pane->attribute(normal_motion->id())->pointer()) : NULL;
     double *ws = normal_motion ? reinterpret_cast<double*>
       ( pane->attribute(_weights->id())->pointer()) : NULL;
     const double *farea = normal_motion ? reinterpret_cast<double*>
       ( pane->attribute(_faceareas->id())->pointer()) : NULL;
     const char *weak = reinterpret_cast<const char*>
       ( pane->attribute(_weak->id())->pointer());
     const int *cnstrs = _cnstr_nodes ? reinterpret_cast<const int*>
       ( pane->attribute(_cnstr_nodes->id())->pointer()) : NULL;
 
     // Loop through real elements of the current pane
     Element_node_enumerator ene( pane, 1);
     for ( int j=0, nj=pane->size_of_real_elements(); j<nj; ++j, ene.next()) {
       int ne = ene.size_of_edges();
       const Point_3  &cnt = fcnts[j];
 
       Vector_3 Ae[3];   // Matrix at an edge
       Ae[0] = Ae[1] = Ae[2] = Vector_3(0,0,0);
       bool is_weak = false;
 
       for ( int k=0; k<ne; ++k) {
         int vindex = ene[k]-1;
         if ( weak[vindex]) { is_weak = true; break; }
       }
       if ( !is_weak) continue;
 
 
       for ( int k=0; k<ne; ++k) {
         int vindex = ene[k]-1;
         
         // Obtain the quadric at the face 
         const Vector_3 *As_v=&As[3*vindex];
 
         double w = get_nonuniform_weight
           ( lambdas[vindex], tranks[vindex]-(cnstrs && cnstrs[vindex]));
 
         Ae[0] += w*As_v[0];
         Ae[1] += w*As_v[1];
         Ae[2] += w*As_v[2];
       }
 
       // Perform eigen-analysis to obtain diffused face normal
       Vector_3 ls;
       compute_eigenvectors( Ae, ls);
       Vector_3 mn_face = Ae[0];
 
       for ( int k=0; k<ne; ++k) {
         int vindex = ene[k]-1;
         
         if ( weak[vindex]) {
           Vector_3 diff = (cnt-pnts[vindex]);
           if ( disps) diff -= disps[vindex];
 
           double offset = diff*mn_face, cos_a=mn_face*vnrms[vindex];
           if (cos_a<0) offset = -offset;
 
           // Use eigenvalue to control amount of dissipation
           double w = farea[j];
           double a = lambdas[vindex][1]/lambdas[vindex][0];
           nmotion[vindex] += a*(w*offset)*vnrms[vindex];
 
           ws[vindex] += w;
         }
       }
     }
   }
 
   _surf->reduce_on_shared_nodes( normal_motion, Manifold::OP_SUM);
   _surf->reduce_on_shared_nodes( _weights, Manifold::OP_SUM);
   Rocblas::div( normal_motion, _weights, normal_motion);
 }

Here is the call graph for this function:

Here is the caller graph for this function:

void distribute_volume_e2n	(	const COM::Attribute *	fvol,
		const COM::Attribute *	tranks,
		COM::Attribute *	vvol
	)

protected

Definition at line 105 of file cons_diff.C.

References Propagation_3::_panes, Propagation_3::_surf, Rocblas::copy_scalar(), i, j, k, Element_node_enumerator::next(), nj, and Element_node_enumerator::size_of_edges().

Referenced by balance_mass().

                                            {
 
   double zero=0;
   Rocblas::copy_scalar( &zero, vvol);
 
   // Loop through the panes and its real faces
   std::vector< COM::Pane*>::iterator it = _panes.begin();
   Manifold::PM_iterator pm_it=_surf->pm_begin();
   for ( int i=0, local_npanes = _panes.size(); 
         i<local_npanes; ++i, ++it, ++pm_it) { 
     COM::Pane *pane = *it;
 
     const double *fv = reinterpret_cast<const double*>
       ( pane->attribute(fvol->id())->pointer());
     double   *vv = reinterpret_cast<double*>
       ( pane->attribute(vvol->id())->pointer());
 
     // Loop through real elements of the current pane
     Element_node_enumerator ene( pane, 1); 
     for ( int j=0, nj=pane->size_of_real_elements(); j<nj; ++j, ene.next()) {
       int ne = ene.size_of_edges();
 
       for (int k=0; k<ne; ++k) {
         int index = ene[k]-1;
         vv[index] += fv[j]; // Omit the factor 1/3 here
       }
     }
   }
   
   // Reduce on shared nodes
   _surf->reduce_on_shared_nodes( vvol, Manifold::OP_SUM);
  
 }

Here is the call graph for this function:

Here is the caller graph for this function:

int eigen_analyze_vertex	(	Vector_3	A_io[3],
		Vector_3 &	b_io,
		double	angle_defect = `0`
	)

protected

Definition at line 240 of file quadric_analysis.C.

References _dir_thres_strong, _dir_thres_weak, _tol_ad, NTS::abs(), and compute_eigenvectors().

Referenced by compute_anisotropic_vertex_centers(), and compute_quadrics().

                                            {
 
   // A_io will be replaced by its eigenvectors.
   Vector_3 *es = A_io;
 
   // Create a backup of b_io, as b_io will be replaced by eigenvalues.
   const Vector_3 b = b_io;
   Vector_3 &lambdas = b_io;
 
   // Compute the eigenvalues and eigenvectors of the matrix. 
   // Note that lambdas will be in decreasing order!
   compute_eigenvectors( es, lambdas);
   
   // Classify offset intersection based on acos(-b.norm()/rv) and lambdas
   if ( std::abs(angle_defect) >= _tol_ad)
     return 3; // Sharp corner
   else if (lambdas[1] < _dir_thres_weak*lambdas[0])
     return 1; // Ridge
   else if ( lambdas[2] < _dir_thres_strong*lambdas[1])
     return 2; // Flat
   else        // Vertices that are neither sharp nor flat, but have
     return 4; // ambiguous ridge directions
 }

Here is the call graph for this function:

Here is the caller graph for this function:

static double eval_angle	(	const Vector_3 &	v1,
		const Vector_3 &	v2
	)

inlinestaticprotected

Definition at line 323 of file FaceOffset_3.h.

References cimg_library::acos(), s, sqrt(), and Vector_3< Type >::squared_norm().

Referenced by build_ridge_neighbor_list(), compute_anisotropic_vertex_centers(), and compute_quadrics().

                                                                     {
     double sqrnrm = v1.squared_norm()*v2.squared_norm();
     if ( sqrnrm==0) return 0;
     double s=v1*v2/std::sqrt( sqrnrm);
     if (s>1) s=1; else if (s<-1) s =-1;
     return std::acos(s);
   }

Here is the call graph for this function:

Here is the caller graph for this function:

static double eval_angle	(	const Vector_2 &	v1,
		const Vector_2 &	v2
	)

inlinestaticprotected

Definition at line 331 of file FaceOffset_3.h.

References cimg_library::acos(), s, sqrt(), and Vector_2< Type >::squared_norm().

                                                                     {
     double sqrnrm = v1.squared_norm()*v2.squared_norm();
     if ( sqrnrm==0) return 0;
     double s=v1*v2/std::sqrt( sqrnrm);
     if (s>1) s=1; else if (s<-1) s =-1;
     return std::acos(s);
   }

Here is the call graph for this function:

PROP_BEGIN_NAMESPACE void filter_and_identify_ridge_edges ( bool filter_curves )

protected

Filter out isolated ridge vertices and identify ridge edges.

Definition at line 33 of file detect_features.C.

References Propagation_3::_buf, _ctangranks, _edges, _eigvecs, _facenormals, Propagation_3::_panes, _ridges, _strong, Propagation_3::_surf, _tangranks, _tol_angle_strong, _tol_angle_weak, _tol_cos_osta, cimg_library::acos(), build_ridge_neighbor_list(), COM_assertion, COM_CHAR, COM_DOUBLE, COM_NC, COMMPI_Initialized(), Rocblas::copy(), Rocblas::copy_scalar(), filter_obscended_ridge(), i, insert_boundary_edges(), j, k, max(), n, Element_node_enumerator::next(), nj, Vector_3< Type >::normalize(), pi, reclassify_ridge_vertices(), and Element_node_enumerator::size_of_edges().

Referenced by time_stepping().

 {
   // Calculate boundary normals
   _surf->update_bd_normals( _facenormals, false);
 
   COM::Attribute *maxdianglev = 
     _buf->new_attribute( "FO_maxdianglev", 'n', COM_DOUBLE, 1, "");
   _buf->resize_array( maxdianglev, 0);
 
   double t=-pi;
   Rocblas::copy_scalar( &t, maxdianglev);
 
   COM::Attribute *mindianglev = 
     _buf->new_attribute( "FO_mindianglev", 'n', COM_DOUBLE, 1, "");
   _buf->resize_array( mindianglev, 0);
   t=pi;
   Rocblas::copy_scalar( &t, mindianglev);
 
   COM::Attribute *maxtrnangv = 
     _buf->new_attribute( "FO_maxtrnangv", 'n', COM_DOUBLE, 1, "");
   _buf->resize_array( maxtrnangv, 0);
 
   t=-2;
   Rocblas::copy_scalar( &t, maxtrnangv);
 
   COM::Attribute *mintrnangv = 
     _buf->new_attribute( "FO_mintrnangv", 'n', COM_DOUBLE, 1, "");
   _buf->resize_array( mintrnangv, 0);
   t=2;
   Rocblas::copy_scalar( &t, mintrnangv);
 
   COM::Attribute *neighbor_features = 
     _buf->new_attribute( "FO_neighbor_features", 'n', COM_CHAR, 1, "");
   _buf->resize_array( neighbor_features, 0);
 
   COM::Attribute *qsedges = 
     _buf->new_attribute( "FO_qsedges", 'e', COM_CHAR, 1, "");
   _buf->resize_array( qsedges, 0);
   _buf->init_done( false);
 
   int zero=0;
   Rocblas::copy_scalar( &zero, _strong);
 
   // Clear ridge edges
   _edges.clear(); _edges.resize( _panes.size());
 
   bool to_filter = filter_curves && _panes.size() == 1 && !COMMPI_Initialized();
 
   // Loop through the panes and its real faces
   std::vector< COM::Pane*>::iterator it = _panes.begin();
   Manifold::PM_iterator pm_it=_surf->pm_begin();
 
   int nedges=0, nvertices=0;
   std::vector<std::map<Edge_ID, double> > edge_maps( _panes.size());
   std::vector<std::map<Edge_ID, double> > diangle_maps( _panes.size());
 
   for ( int i=0, local_npanes = _panes.size();
         i<local_npanes; ++i, ++it, ++pm_it) { 
     COM::Pane *pane = *it;
     std::map< Edge_ID, double> &emap_pn = edge_maps[i];
     std::map< Edge_ID, double> &diangle_pn = diangle_maps[i];
 
     nvertices += pane->size_of_real_nodes();
 
     // Obtain nodal coordinates of current pane, assuming contiguous layout
     const Point_3 *pnts = reinterpret_cast<Point_3*>
       (pane->attribute(COM_NC)->pointer());
 
     const Vector_3 *es = reinterpret_cast<const Vector_3*>
       ( pane->attribute(_eigvecs->id())->pointer());
     char *tranks = reinterpret_cast<char*>
       ( pane->attribute(_tangranks->id())->pointer());
     char *strong = reinterpret_cast<char*>
       ( pane->attribute(_strong->id())->pointer());
     double *minvs = reinterpret_cast<double*>
       ( pane->attribute(mindianglev->id())->pointer());
     double *maxvs = reinterpret_cast<double*>
       ( pane->attribute(maxdianglev->id())->pointer());
     double *minta = reinterpret_cast<double*>
       ( pane->attribute(mintrnangv->id())->pointer());
     double *maxta = reinterpret_cast<double*>
       ( pane->attribute(maxtrnangv->id())->pointer());
     Vector_3 *fnrms = reinterpret_cast<Vector_3*>
       ( pane->attribute(_facenormals->id())->pointer());
 
     // Loop through real elements of the current pane
     Element_node_enumerator ene( pane, 1);
     for ( int j=0, nj=pane->size_of_real_elements(); j<nj; ++j, ene.next()) {
       int ne=ene.size_of_edges(), uindex=ene[ne-1]-1, vindex=ene[0]-1, windex=0;
       nedges += ne;
         
       for ( int k=0; k<ne; ++k, uindex=vindex, vindex=windex) {
         windex = ene[ (k+1)%ne]-1;
 
         Edge_ID eid(j+1,k), eid_opp = pm_it->get_opposite_real_edge( eid);
 
         // Consider border edges as flat and skip them when identifying ridges
         if ( pm_it->is_physical_border_edge( eid_opp)) 
         { ++nedges; continue; }
 
         Vector_3 tng = pnts[ windex] - pnts[vindex];
         tng.normalize();
         double feg = es[3*vindex+2]* tng;
 
         // Safeguard strong angles.
         // Important: We assume this computation gives identical results 
         // for eid and eid_opp near the thresholds.
         const Vector_3 &n1 = fnrms[j];
         const Vector_3 &n2 = eid_opp.is_border() ? 
           pm_it->get_bd_normal( eid_opp) : fnrms[eid_opp.eid()-1];
         double costheta = n1*n2;
         if ( costheta>1) costheta = 1;
         else if (costheta<-1) costheta = -1;
 
         double fangle = std::acos( costheta);
         if ( fangle >  _tol_angle_strong) {
           // Upgrade vertices.
           if ( tranks[vindex]==2) tranks[vindex] = 1;
           if ( tranks[windex]==2) tranks[windex] = 1;
 
           if ( fangle>pi*0.505 && fangle !=pi)
             strong[vindex] = strong[windex] = 2;
           else
             strong[vindex] = strong[windex] = 1;
 
           std::map< Edge_ID, double>::iterator rit = diangle_pn.find(eid);
           // Make sure the fangle is used consistently.
           if ( rit == diangle_pn.end()) diangle_pn[eid] = fangle;
           else fangle = rit->second;
  
           if ( feg>0) {
             if ( fangle > maxvs[vindex]) maxvs[vindex] = fangle; 
             if ( feg > maxta[vindex]) maxta[vindex] = feg;
             emap_pn[ eid] = feg;
           }
           else {
             if ( -fangle<minvs[vindex]) minvs[vindex] = -fangle;   
             if ( feg < minta[vindex]) minta[vindex] = feg;
 
             emap_pn[ eid] = feg;
           }
         }
         else if ( fangle > _tol_angle_weak) {
           if ( feg > 0) {
             if ( fangle > maxvs[vindex]) maxvs[vindex] = fangle;
             if ( feg > maxta[vindex])    maxta[vindex] = feg;
           }
           else {
             if ( -fangle < minvs[vindex]) minvs[vindex] = -fangle;
             if ( feg < minta[vindex])    minta[vindex] = feg;
           }
 
           // Save the current edge for comparison
           diangle_pn[eid] = fangle; emap_pn[ eid] = feg;
         }
       } // for k (edges)
     } // for j (elements)
   } // for i (panes)
 
   if ( to_filter)
     std::cout << "There are " << nedges/2 << " edges and " 
               << nvertices << " vertices " << std::endl;
   
   // Reduce on shared nodes for mindianglev and maxdianglev
   _surf->reduce_on_shared_nodes( mindianglev, Manifold::OP_MIN);
   _surf->reduce_on_shared_nodes( maxdianglev, Manifold::OP_MAX);
 
   _surf->reduce_on_shared_nodes( mintrnangv, Manifold::OP_MIN);
   _surf->reduce_on_shared_nodes( maxtrnangv, Manifold::OP_MAX);
   _surf->reduce_on_shared_nodes( _strong, Manifold::OP_MAX);
 
   // Count the number of incident strongly attached edges.
   it = _panes.begin();
   pm_it=_surf->pm_begin();
   Rocblas::copy_scalar( &zero, neighbor_features);
 
   for ( int i=0, local_npanes = _panes.size(); i<local_npanes; 
         ++i, ++it, ++pm_it) { 
     COM::Pane *pane = *it;
     std::map< Edge_ID, double> &emap_pn = edge_maps[i];
     std::map< Edge_ID, double> &diangle_pn = diangle_maps[i];
 
     const double *minvs = reinterpret_cast<double*>
       ( pane->attribute(mindianglev->id())->pointer());
     const double *maxvs = reinterpret_cast<double*>
       ( pane->attribute(maxdianglev->id())->pointer());
     const double *minta = reinterpret_cast<double*>
       ( pane->attribute(mintrnangv->id())->pointer());
     double *maxta = reinterpret_cast<double*>
       ( pane->attribute(maxtrnangv->id())->pointer());
 
     char *tranks = reinterpret_cast<char*>
       ( pane->attribute(_tangranks->id())->pointer());
     char *nnf = reinterpret_cast<char*>
       ( pane->attribute(neighbor_features->id())->pointer());
 
     COM_assertion( emap_pn.size() == diangle_pn.size());
     // Loop through ebuf and put ridges edges onto emap_pn
     for ( std::map< Edge_ID, double>::iterator eit=emap_pn.begin(), 
             rit=diangle_pn.begin(), eend=emap_pn.end(); eit!=eend; ++eit, ++rit) {
       Edge_ID eid = eit->first;
       double feg = eit->second, fangle = rit->second;
 
       COM_assertion(! eid.is_border());
       Element_node_enumerator ene( pane, eid.eid()); 
       int vi = eid.lid(), vindex = ene[vi]-1;
 
       // Determine the number of strongly attached edges
       if ( !tranks[vindex] || tranks[vindex] == char(-1) || fangle > _tol_angle_strong ||
            ( feg < 0 && fangle == -minvs[vindex] || 
              feg > 0 && fangle == maxvs[vindex]) &&
            ( feg < -_tol_cos_osta && feg == minta[vindex] || 
              feg > _tol_cos_osta && feg == maxta[vindex])) {
         ++nnf[vindex];
       }
     } // for eit
   } // for i
    
   _surf->reduce_on_shared_nodes( neighbor_features, Manifold::OP_SUM);
 
   it = _panes.begin();
   pm_it=_surf->pm_begin();
 
   for ( int i=0, local_npanes = _panes.size(); i<local_npanes; 
         ++i, ++it, ++pm_it) { 
     COM::Pane *pane = *it;
     std::map< Edge_ID, double> &emap_pn = edge_maps[i];
     std::map< Edge_ID, double> &diangle_pn = diangle_maps[i];
     std::set< Edge_ID> &eset_pn = _edges[i];
 
     const double *minvs = reinterpret_cast<double*>
       ( pane->attribute(mindianglev->id())->pointer());
     const double *maxvs = reinterpret_cast<double*>
       ( pane->attribute(maxdianglev->id())->pointer());
     const double *minta = reinterpret_cast<double*>
       ( pane->attribute(mintrnangv->id())->pointer());
     double *maxta = reinterpret_cast<double*>
       ( pane->attribute(maxtrnangv->id())->pointer());
     
     char *strong = reinterpret_cast<char*>
       ( pane->attribute(_strong->id())->pointer());
     char *tranks = reinterpret_cast<char*>
       ( pane->attribute(_tangranks->id())->pointer());
     char *qses = reinterpret_cast<char*>
       ( pane->attribute(qsedges->id())->pointer());
     char *nnf = reinterpret_cast<char*>
       ( pane->attribute(neighbor_features->id())->pointer());
  
     COM_assertion( emap_pn.size() == diangle_pn.size());
     // Loop through ebuf and put ridges edges onto emap_pn
     for ( std::map< Edge_ID, double>::const_iterator eit=emap_pn.begin(), 
             rit=diangle_pn.begin(), eend=emap_pn.end(); eit!=eend; ++eit, ++rit) {
       Edge_ID eid = eit->first, eid_opp = pm_it->get_opposite_real_edge( eid);
       COM_assertion( !eid.is_border());
       double feg = eit->second, fangle = rit->second, feg_opp;
 
       Element_node_enumerator ene( pane, eid.eid());
       int vi=eid.lid(), ne = ene.size_of_edges();
       int vindex = ene[vi]-1, windex = ene[ (vi+1)%ne]-1;
 
       if ( (nnf[vindex] && nnf[windex]) &&
            ( strong[vindex] || !tranks[vindex] || tranks[vindex] == char(-1) ||
              feg < 0 && fangle == -minvs[vindex] || 
              feg > 0 && fangle == maxvs[vindex] ||
              feg < -_tol_cos_osta && feg == minta[vindex] ||
              feg > _tol_cos_osta && feg == maxta[vindex]) ||
            nnf[vindex]==2 && nnf[windex]==2 && 
            ( (feg_opp=emap_pn.find( eid_opp)->second)< 0 && 
              fangle == -minvs[windex] ||
              feg_opp > 0 && fangle == maxvs[windex]) &&
            ( feg_opp < -_tol_cos_osta && feg_opp == minta[windex] ||
              feg_opp > _tol_cos_osta && feg_opp == maxta[windex])) {
         eset_pn.insert( eid);
         qses[ eid.eid()-1] |= (1<<eid.lid());
       }
     } // for eit
   } // for i
 
   // Update border-edge count
   _surf->update_bdedge_bitmap( qsedges);
 
   // Select the edges whose both half-edges are quasi-strong.
   it = _panes.begin();
   pm_it=_surf->pm_begin();
 
   for ( int i=0, local_npanes = _panes.size(); i<local_npanes; 
         ++i, ++it, ++pm_it) {
     std::set< Edge_ID> &eset_pn = _edges[i];
 
     // Select quasi-strong edges
     for ( std::set< Edge_ID>::iterator eit=eset_pn.begin(); 
           eit!=eset_pn.end();) {
       Edge_ID eid = *eit, eid_opp = pm_it->get_opposite_real_edge( eid);
 
       if ( !eid_opp.is_border() && eset_pn.find(eid_opp)==eset_pn.end() ||
            eid_opp.is_border() && !pm_it->get_bdedge_flag( eid_opp)) {
         std::set< Edge_ID>::iterator to_delete = eit;
         ++eit;
         eset_pn.erase( to_delete);
       }
       else 
         ++eit;
     }
   }
 
   // Note: neighbor_features is used as buffers in reclassify_ridge_vertices
   reclassify_ridge_vertices( 0, 1, neighbor_features, 
                              _tangranks, to_filter);
   
   // Currently the filtering step does not yet work in parallel.
   // Filter out obscure-ended ridges. Currently supports only sequential meshes
   if ( to_filter) for ( ;;) {
     std::vector< ObscendSet > obscends;
 
     int nobscended=build_ridge_neighbor_list( maxtrnangv,mintrnangv, edge_maps,
                                               maxdianglev,mindianglev, 
                                               diangle_maps, obscends);
     if ( nobscended==0) break;
 
     bool dropped = filter_obscended_ridge( edge_maps, diangle_maps, obscends);
     if ( !dropped) break;
 
     reclassify_ridge_vertices( 0, 0, neighbor_features, _tangranks, 1);
   }
   // Reclassify ridge edges to upgrade potentially missing corners.
   reclassify_ridge_vertices( 1, 0, neighbor_features, _tangranks, to_filter);
 
   Rocblas::copy( _tangranks, _ctangranks);
   // Update vertices on constraint boundary.
   if ( insert_boundary_edges( _ctangranks))
     // Reclassify ridge edges to upgrade potentially missing corners.
     reclassify_ridge_vertices(1, 1, neighbor_features, _ctangranks, to_filter);
 
   int nredges=0; // Number of ridge edges
   // Export ridge edges into list
   it = _panes.begin();
   pm_it=_surf->pm_begin();
   for ( int i=0, local_npanes = _panes.size(); 
         i<local_npanes; ++i, ++it, ++pm_it) { 
     COM::Pane *pane = *it;
     std::set< Edge_ID> &eset_pn = _edges[i];
 
     // Allocate memory for _ridges and save them. 
     COM::Attribute *att_rdg = pane->attribute(_ridges->id());
 
     int n=eset_pn.size();
     att_rdg->allocate( 2, std::max(n,att_rdg->size_of_items()), 0);
     int *rdgs = reinterpret_cast<int*>
       ( pane->attribute(_ridges->id())->pointer()), ii=0;
 
     for ( std::set< Edge_ID>::const_iterator eit=eset_pn.begin();
           eit!=eset_pn.end(); ++eit) {
       // Make sure that eid_opp is not border
       Edge_ID eid = *eit, eid_opp = pm_it->get_opposite_real_edge( eid);
 
       // Obtain vertex indices.
       Element_node_enumerator ene( pane, eid.eid()); 
       int vindex = eid.lid(), neighbor = (vindex+1)%ene.size_of_edges();
       vindex = ene[vindex]; neighbor=ene[neighbor];
      
       if ( vindex<neighbor || eid_opp.is_border())
       { rdgs[ii++] = vindex; rdgs[ii++] = neighbor; }
     }
     
     att_rdg->set_size( ii/2);
     nredges += att_rdg->size_of_items();
   }
 
   if ( to_filter)
     std::cout << "Found " << nredges << " ridge edges " << std::endl;
 
   // Deallocate mindianglev and maxdianglev
   _buf->delete_attribute( qsedges->name());
   _buf->delete_attribute( neighbor_features->name());
   _buf->delete_attribute( mintrnangv->name());
   _buf->delete_attribute( maxtrnangv->name());
   _buf->delete_attribute( mindianglev->name());
   _buf->delete_attribute( maxdianglev->name());
   _buf->init_done( false);
 }

Here is the call graph for this function:

Here is the caller graph for this function:

bool filter_obscended_ridge	(	const std::vector< std::map< Edge_ID, double > > &	edge_maps,
		const std::vector< std::map< Edge_ID, double > > &	diangl_maps,
		const std::vector< ObscendSet > &	obscends
	)

protected

Filter out weak-ended curves.

Definition at line 702 of file detect_features.C.

References Propagation_3::_panes, _ridgeneighbors, _strong, Propagation_3::_surf, _tangranks, _tol_kangle, _tol_kstrong, RidgeNeighbor::eid, i, remove_obscure_curves(), and RidgeNeighbor::vid.

Referenced by filter_and_identify_ridge_edges().

                                                                  {
   // Loop through the panes and ridge edges
   std::vector< COM::Pane*>::iterator it = _panes.begin();
   Manifold::PM_iterator pm_it=_surf->pm_begin();
   int local_npanes = _panes.size();
   std::vector< ObscendSet > todelete(local_npanes);
 
   for ( int i=0; i<local_npanes; ++i, ++it, ++pm_it) { 
     const ObscendSet &obscend_pn = obscends[i];
     ObscendSet &todelete_pn = todelete[i];
     if ( obscend_pn.empty()) continue;
 
     COM::Pane *pane = *it;
     const char *tranks = reinterpret_cast<char*>
       ( pane->attribute(_tangranks->id())->pointer());
     char *strong = reinterpret_cast<char*>
       ( pane->attribute(_strong->id())->pointer());
     const RidgeNeighbor *rdgngbs = reinterpret_cast<RidgeNeighbor*>
       ( pane->attribute(_ridgeneighbors->id())->pointer());
     const std::map< Edge_ID, double> &diangle_pn = diangle_maps[i];
 
     for ( ObscendSet::const_iterator rit=obscend_pn.begin(),
             rend=obscend_pn.end(); rit!=rend; ++rit) {
       // Loop through edges to count
       int cur=rit->second.first, next=rit->second.second, start=cur;
 
       Edge_ID eid = rit->first;
       bool is_dangling = ( rdgngbs[ (cur-1)*2].vid==next && 
                            !rdgngbs[ (cur-1)*2+1].vid ||
                            !rdgngbs[ (cur-1)*2].vid && 
                            !rdgngbs[ (cur-1)*2+1].vid);
       bool is_not_joint1 = ( rdgngbs[ (cur-1)*2].vid!=next && 
                              rdgngbs[ (cur-1)*2+1].vid!=next), is_not_joint2=false;
 
       if ( !is_dangling && !is_not_joint1) continue;
       bool is_corner=false, is_strong1=strong[cur-1] || !tranks[cur-1];
       bool is_strong2=strong[next-1] || !tranks[next-1];
 
       int count_ks = 0;
       for ( ;!(is_corner = !tranks[next-1]) && next && next!=start; ) {
         count_ks += diangle_pn.find( eid)->second>_tol_kangle;
         
         int  index = (next-1)*2;
         if ( rdgngbs[ index].vid == cur) { 
           cur = next; next = rdgngbs[index+1].vid; eid = rdgngbs[index+1].eid;
           is_strong2=strong[next-1] || !tranks[next-1];
           is_not_joint2 = false; 
         }
         else {
           Edge_ID eid_opp = pm_it->get_opposite_real_edge( eid);
 
           is_not_joint2 = (rdgngbs[ index+1].vid!=cur) && 
             (obscend_pn.find(eid_opp) != obscend_pn.end());
 
           if ( rdgngbs[ index+1].vid == cur) { 
             cur = next; next = rdgngbs[index].vid; eid = rdgngbs[ index].eid; 
             is_strong2=strong[next-1] || !tranks[next-1];
           }
           else
           { cur = next; next = 0; }
         }
       }
 
       // Check whether the edge is quasi-obscure.
       if ( (is_dangling || is_not_joint1 && is_not_joint2) && count_ks<=_tol_kstrong ||
            is_strong1 && is_strong2 && count_ks==0) {
         todelete_pn.insert( *rit);
       }
     }
   }
   
   return remove_obscure_curves( todelete);
 }

Here is the call graph for this function:

Here is the caller graph for this function:

double get_courant_constant ( ) const

inline

Get the constant in front of CFL condition (0<=c<=1)

Definition at line 87 of file FaceOffset_3.h.

References _courant.

87 { return _courant; }

FaceOffset_3::_courant

double _courant

Definition: FaceOffset_3.h:378

double get_eigen_threshold ( ) const

inline

Get the threshold of eigenvalues for splitting primary and null spaces.

Definition at line 80 of file FaceOffset_3.h.

References _eig_thres.

80 { return _eig_thres; }

FaceOffset_3::_eig_thres

double _eig_thres

Definition: FaceOffset_3.h:394

double get_fangle_strong_radians ( ) const

inline

Get the threshold for strong-feature angle.

Definition at line 114 of file FaceOffset_3.h.

References _tol_angle_strong.

115 { return _tol_angle_strong; }

FaceOffset_3::_tol_angle_strong

double _tol_angle_strong

Definition: FaceOffset_3.h:385

double get_fangle_turn_radians ( ) const

inline

Set the threshold for weak-feature angle.

Definition at line 118 of file FaceOffset_3.h.

References _tol_cos_osta, and cimg_library::acos().

119 { return std::acos(_tol_cos_osta)*2; }

FaceOffset_3::_tol_cos_osta

double _tol_cos_osta

Definition: FaceOffset_3.h:391

cimg_library::acos

CImg< _cimg_Tfloat > acos(const CImg< T > &instance)

Definition: CImg.h:6051

Here is the call graph for this function:

double get_fangle_weak_radians ( ) const

inline

Get the threshold for weak-feature angle.

Definition at line 110 of file FaceOffset_3.h.

References _tol_angle_weak.

111 { return _tol_angle_weak; }

FaceOffset_3::_tol_angle_weak

double _tol_angle_weak

Definition: FaceOffset_3.h:384

double get_nonuniform_weight	(	const Vector_3 &	lambdas,
		int	trank
	)

protected

Definition at line 47 of file AnisotropicSmoothing.C.

Referenced by denoise_surface().

                                                            {
 
   const double eps=1.e-6;
   double w = 1/lambdas[0];
   
   if (trank == 2) return w;
   else if ( trank==1 ) return eps*w;
   else return eps*eps*w;
 }

Here is the caller graph for this function:

int get_normal_diffuion ( ) const

inline

Get number of iterations for normal diffusion.

Definition at line 73 of file FaceOffset_3.h.

References _nrm_dfsn.

73 { return _nrm_dfsn; }

FaceOffset_3::_nrm_dfsn

char _nrm_dfsn

Definition: FaceOffset_3.h:375

void get_primary_component	(	const Vector_3 &	nrm,
		const Vector_3	es[],
		int	trank,
		Vector_3 &	prim
	)		const

inlineprotected

Definition at line 254 of file FaceOffset_3.h.

References i.

                                                                {
     prim = Vector_3(0,0,0);
 
     for ( int i=0; i<3-trank; ++i) prim += nrm*es[i]*es[i];
   }

void get_tangents	(	const Vector_3 *	es,
		const Vector_3 &	lambdas,
		const Vector_3 &	mn,
		int	nrank,
		Vector_3	vs[2],
		double	scales[2],
		int *	is_strong = `NULL`
	)		const

protected

Definition at line 58 of file AnisotropicSmoothing.C.

References _dir_thres_strong, is_acute_corner(), is_acute_ridge(), max(), min(), and sqrt().

Referenced by compute_anisotropic_vertex_centers().

                                                                       {
   vs[0] = es[1];
   vs[1] = es[2];
 
   // Bound from two ends appear to be the most robust.
   const double max_scale=0.07, min_scale=0.005;
 
   if ( scales) {
     scales[0] = std::min( max_scale, std::max(lambdas[1]/lambdas[0], min_scale));
     scales[1] = std::min( max_scale, std::max(lambdas[2]/lambdas[0], min_scale));
 
     scales[0] = std::sqrt( scales[0]);
     scales[1] = std::sqrt( scales[1]);
     
     if ( is_strong) {
       bool acute_ridge = FaceOffset_3::is_acute_ridge( es, lambdas, mn, nrank);
       bool acute_corner = FaceOffset_3::is_acute_corner( es, lambdas, mn, nrank);
 
       *is_strong = acute_ridge || acute_corner || 
         lambdas[1]>_dir_thres_strong*lambdas[0];
     }
   }
 }

Here is the call graph for this function:

Here is the caller graph for this function:

bool get_wavefrontal_expansion ( ) const

inline

Obtain the wavefrontal switch.

Definition at line 64 of file FaceOffset_3.h.

References _wf_expn.

64 { return _wf_expn; }

FaceOffset_3::_wf_expn

bool _wf_expn

Definition: FaceOffset_3.h:373

int get_weighting_scheme ( ) const

inline

Definition at line 91 of file FaceOffset_3.h.

References _wght_scheme.

91 { return _wght_scheme; }

FaceOffset_3::_wght_scheme

char _wght_scheme

Definition: FaceOffset_3.h:374

int insert_boundary_edges ( COM::Attribute * tr_attr )

protected

Definition at line 415 of file detect_features.C.

References Propagation_3::_buf, Propagation_3::_cnstr_bndry_edges, _edges, Propagation_3::_panes, Propagation_3::_surf, COMMPI_Initialized(), i, iend, j, k, MPI_SUM, Element_node_enumerator::next(), nj, and Element_node_enumerator::size_of_edges().

Referenced by filter_and_identify_ridge_edges().

                                               {
   std::vector< COM::Pane*>::iterator it = _panes.begin(), iend=_panes.end();
   Manifold::PM_iterator pm_it=_surf->pm_begin();
 
   int inserted=0;
   // Loop through the panes and its real elements
   for ( int i=0; it!=iend; ++i, ++it, ++pm_it) {
     COM::Pane *pane = *it;
     std::set< Edge_ID> &eset_pn = _edges[i];
 
     char *tranks = reinterpret_cast<char*>
       ( pane->attribute(tr_attr->id())->pointer());
     const char *val_bndry=(const char*)pane->attribute
       (_cnstr_bndry_edges->id())->pointer();
 
     // Loop through real elements of the current pane
     Element_node_enumerator ene( pane, 1);
     for ( int j=0, nj=pane->size_of_real_elements(); j<nj; ++j, ene.next()) {
       for ( int k=0, ne=ene.size_of_edges(); k<ne; ++k) {
         Edge_ID eid(j+1,k), eid_opp = pm_it->get_opposite_real_edge( eid);
         
         // If a vertex has mixed geometric and topological edges, 
         // set it to be a corner
         if ( pm_it->is_physical_border_edge( eid_opp) || 
              (val_bndry && (val_bndry[j] & (1<<k)))) {
           int vindex = ene[k]-1;
 
           if (tranks[ vindex]==2 || eset_pn.find(eid)==eset_pn.end()) { 
             if ( tranks[ vindex]==1) tranks[ vindex] = 0;
             eset_pn.insert(eid); ++inserted; 
           }
         }
       } // for k
     } // for j
   } // for it
 
   int inserted_total = inserted;
   if ( COMMPI_Initialized()) {
     MPI_Allreduce( &inserted, &inserted_total, 1, MPI_INT, MPI_SUM, 
                    _buf->get_communicator());
   }
 
   return inserted_total;
 }

Here is the call graph for this function:

Here is the caller graph for this function:

bool is_acute_corner	(	const Vector_3 *	es,
		const Vector_3 &	lambdas,
		const Vector_3 &	mn,
		int	nrank
	)		const

inlineprotected

Definition at line 357 of file FaceOffset_3.h.

References NTS::abs().

Referenced by get_tangents(), and is_strong().

                                                              {
     
     return nrank==3 && ( std::abs(es[0]*mn)<=std::abs(es[1]*mn) ||
                          std::abs(es[0]*mn)<=std::abs(es[2]*mn));
   }

Here is the call graph for this function:

Here is the caller graph for this function:

bool is_acute_ridge	(	const Vector_3 *	es,
		const Vector_3 &	lambdas,
		const Vector_3 &	mn,
		int	nrank
	)		const

inlineprotected

Definition at line 352 of file FaceOffset_3.h.

References NTS::abs().

Referenced by get_tangents(), and is_strong().

                                                             {
     return  nrank!=1 && std::abs(es[0]*mn)<=std::abs(es[1]*mn);
   }

Here is the call graph for this function:

Here is the caller graph for this function:

bool is_strong	(	const Vector_3 *	es,
		const Vector_3 &	lambdas,
		const Vector_3 &	mn,
		int	nrank
	)		const

protected

Definition at line 36 of file AnisotropicSmoothing.C.

References _dir_thres_strong, is_acute_corner(), and is_acute_ridge().

                                                 {
   if ( lambdas[1]>_dir_thres_strong*lambdas[0]) return true;
   
   return is_acute_ridge( es, lambdas, mn, nrank) || 
     is_acute_corner( es, lambdas, mn, nrank);
 }

Here is the call graph for this function:

void mark_weak_vertices ( )

protected

Mark vertices that are weak.

Definition at line 356 of file quadric_analysis.C.

References _eig_thres, _eig_weak_lb, _eig_weak_ub, _eigvalues, Propagation_3::_panes, Propagation_3::_surf, _tangranks, _weak, Rocblas::copy_scalar(), i, j, and nj.

Referenced by time_stepping().

                                       {
   char zero_c=0;
   Rocblas::copy_scalar( &zero_c, _weak);
 
   std::vector< COM::Pane*>::iterator it = _panes.begin();
   Manifold::PM_iterator pm_it=_surf->pm_begin();
   for ( int i=0, local_npanes = _panes.size(); 
         i<local_npanes; ++i, ++it, ++pm_it) { 
     COM::Pane *pane = *it;
     char *tranks = reinterpret_cast<char*>
       ( pane->attribute(_tangranks->id())->pointer());
     char *weak = reinterpret_cast<char*>
       ( pane->attribute(_weak->id())->pointer());
     const Vector_3 *lambdas = reinterpret_cast<const Vector_3*>
       ( pane->attribute(_eigvalues->id())->pointer());
 
     for ( int j=0, nj=pane->size_of_real_nodes(); j<nj; ++j) {
       weak[j] = tranks[j]==2 && 
         lambdas[j][1] >= _eig_weak_lb*lambdas[j][0] &&
         lambdas[j][1] <= _eig_weak_ub*lambdas[j][0] ||
         tranks[j]<=1 && _eig_thres*lambdas[j][0] > lambdas[j][1];
     }
   }
 
   _surf->reduce_on_shared_nodes( _weak, Manifold::OP_MIN);
 }

Here is the call graph for this function:

Here is the caller graph for this function:

PROP_BEGIN_NAMESPACE void nulaplacian_smooth	(	const COM::Attribute *	vert_normals,
		const COM::Attribute *	tangranks,
		const COM::Attribute *	disps,
		COM::Attribute *	vcenters,
		COM::Attribute *	buf,
		COM::Attribute *	vert_weights_buf,
		const COM::Attribute *	edge_weights = `NULL`
	)

protected

Redistribute smooth vertices by NuLaplacian smoothing.

Definition at line 36 of file NuLaplacian.C.

References _eigvecs, Propagation_3::_panes, Propagation_3::_surf, NTS::abs(), COM_assertion_msg, COM_NC, Rocblas::copy_scalar(), Vector_3< Type >::cross_product(), i, j, k, Element_node_enumerator::next(), nj, Vector_3< Type >::normalize(), and Element_node_enumerator::size_of_edges().

Referenced by time_stepping().

                                                       {
   const double alpha = 0.03;
 
   // First, sum up weights for each vertex.
   double eps = 1.e-100;
   Rocblas::copy_scalar( &eps, vert_weights_buf);
 
   // Loop through the panes and its real faces
   std::vector< COM::Pane*>::iterator it = _panes.begin();
   for ( int i=0, local_npanes = _panes.size(); i<local_npanes; ++i, ++it) { 
     COM::Pane *pane = *it; 
     
     // Obtain address for buffer space for vertex weights
     double *vws = reinterpret_cast<double*>
       (pane->attribute(vert_weights_buf->id())->pointer());
     // Obtain address for edge weights
     const double *ews = edge_weights ? reinterpret_cast<const double*>
       (pane->attribute(edge_weights->id())->pointer()) : NULL;
 
     Element_node_enumerator ene( pane, 1); 
     // Loop through the real elements of the pane
     for ( int j=0, nj=pane->size_of_real_elements(); j<nj; ++j, ene.next()) {
       int ne = ene.size_of_edges();     
 
       COM_assertion_msg( ne==4, "NuLapacian supports only quadrilateral meshes");
 
       // Loop through the vertices of the element to compute weights
       for ( int k=0; k<ne; ++k) {
         vws[ene[k]-1] += ews ? ews[j*4+k] + ews[j*4+(k+3)%4] : 1;
       }
     } 
   } 
   _surf->reduce_on_shared_nodes( vert_weights_buf, Manifold::OP_SUM);
 
   // Zero out buf
   double zero=0;
   Rocblas::copy_scalar( &zero, buf);
 
   // Loop through the panes and its real faces
   it = _panes.begin(); 
   for ( int i=0, local_npanes = _panes.size(); i<local_npanes; ++i, ++it) { 
     COM::Pane *pane = *it; 
     
     // Obtain nodal coordinates of current pane, assuming contiguous layout
     const Vector_3 *pnts = reinterpret_cast<Vector_3*>
       (pane->attribute(COM_NC)->pointer());
     const Vector_3 *ds = reinterpret_cast<Vector_3*>
       (pane->attribute(disps->id())->pointer());
     const Vector_3 *nrms = reinterpret_cast<Vector_3*>
       (pane->attribute(vert_normals->id())->pointer());
     const Vector_3 *vcnt = reinterpret_cast<Vector_3*>
       (pane->attribute(vcenters->id())->pointer());
     
     // Obtain address for output displacements
     Vector_3 *dbuf = reinterpret_cast<Vector_3*>
       (pane->attribute(buf->id())->pointer());
     // Obtain address for buffer space for vertex weights
     const double *vws = reinterpret_cast<const double*>
       (pane->attribute(vert_weights_buf->id())->pointer());
     // Obtain address for edge weights
     const double *ews = edge_weights ? reinterpret_cast<const double*>
       (pane->attribute(edge_weights->id())->pointer()) : NULL;
 
     // pnts_buf stores difference between center and perturbed position.
     Vector_3 pnts_buf[4];
     int      is_regular[4];
 
     Element_node_enumerator ene( pane, 1); 
     // Loop through the real elements of the pane
     for ( int j=0, nj=pane->size_of_real_elements(); j<nj; ++j, ene.next()) {
       int ne = ene.size_of_edges();     
 
       COM_assertion_msg( ne==4, 
                          "NuLapacian now supports only quadrilateral meshes");
 
       // Compute face center of the face perturbed along normal direction.
       Vector_3 cnt(0,0,0);
       for ( int k=0; k<ne; ++k) {
         int uindex = ene[k]-1;
 
         is_regular[k] = std::abs(vws[uindex]-4)<1.e-6;
         pnts_buf[k] = pnts[uindex] + ds[uindex];
 
         cnt += pnts_buf[k];
       }
       cnt *= 0.25; // Divide cnt by ne
 
       // pnts_buf contains the displacements from the face center to the vertex
       for ( int k=0; k<ne; ++k) pnts_buf[k] -= cnt;
       
       // Loop through the edges of the element
       int uindex=ene[ne-1]-1, vindex=ene[0]-1;
       for ( int k=0; k<ne; ++k, uindex=vindex, vindex=ene[(k==ne)?0:k]-1) {
         // Compute normal of the perturbed edge as the cross product 
         // of the vertex normal and edge tangent
         Vector_3 tng = pnts[uindex]-pnts[vindex]+ vcnt[uindex]-vcnt[vindex];
         double w = ews ? ews[j*4+k] : 0.5;
 
         int k_1 = (k+3)%4;
         if ( is_regular[k_1]) {
           // Accumulate buf and vert_weights for vertex uindex
           Vector_3 nrm_e=Vector_3::cross_product(tng, nrms[uindex]).normalize();
           
           dbuf[uindex] -= w * ( pnts_buf[k_1] * nrm_e * nrm_e +
                                 alpha * pnts_buf[k_1]);
         }
         
         if ( is_regular[k]) {
           // Accumulate buf and vert_weights for vertex vindex
           Vector_3 nrm_e=Vector_3::cross_product(tng, nrms[vindex]).normalize();
 
           dbuf[vindex] -= w * ( pnts_buf[k] * nrm_e * nrm_e +
                                 alpha * pnts_buf[k]);
         }
       }
     } 
   }
 
   // Reduce on shared nodes for buf and vert_weights_buf
   _surf->reduce_on_shared_nodes( buf, Manifold::OP_SUM);
 
   // Loop through all real vertices to update displacement vector
   it = _panes.begin();
   for (int i=0, local_npanes = _panes.size(); i<local_npanes; ++i, ++it) { 
     COM::Pane *pane = *it;
 
     const char *tranks = reinterpret_cast<const char*>
       ( pane->attribute( tangranks->id())->pointer());
     const double *vws = reinterpret_cast<const double*>
       ( pane->attribute( vert_weights_buf->id())->pointer());
 
     Vector_3 *vcnt = reinterpret_cast<Vector_3*>
       ( pane->attribute( vcenters->id())->pointer());
     Vector_3 *dbuf = reinterpret_cast<Vector_3*>
       ( pane->attribute( buf->id())->pointer());
 
     const Vector_3 *es = reinterpret_cast<const Vector_3*>
       ( pane->attribute(_eigvecs->id())->pointer());
 
     // Loop through all real nodes of the pane
     for ( int j=0, nj=pane->size_of_real_nodes(); j<nj; ++j) {
       // If vertex is unmasked and is regular
       if ( std::abs(vws[j]-4) < 1.e-6) {
         // Assign the displacement to be the weighted average of projections
         vcnt[j] = dbuf[j] / vws[j];
         
         if ( tranks[j] == 1) // Project onto the direction
           vcnt[j] = (vcnt[j]*es[3*j+2])*es[3*j+2];
         else if ( tranks[j] == 2) // Project onto the plane
           vcnt[j] -= (vcnt[j]*es[3*j])*es[3*j];
         else
           vcnt[j] = Vector_3(0,0,0);
       }
     }
   }
 }

Here is the call graph for this function:

Here is the caller graph for this function:

void numerical_dissipation ( COM::Attribute * nodal_buffer )

protected

Introduce numerical dissipation into equation.

bool obtain_constrained_directions	(	COM::Attribute *	disps,
		COM::Attribute *	vcenters
	)

protected

Decompose propagation directions based on constraints.

Definition at line 506 of file FaceOffset_3.C.

References _As, _bmedial, _boffset, Propagation_3::_cnstr_bndry_nodes, Propagation_3::_cnstr_nodes, _conserv, _eig_thres, _eigvalues, _eigvecs, Propagation_3::_panes, Propagation_3::_surf, _tangranks, NTS::abs(), balance_mass(), Vector_3< Type >::cross_product(), Mesquite::det(), Propagation_3::get_constraint_directions(), i, j, k, Vector_3< Type >::normalize(), and obtain_directions().

Referenced by time_stepping().

                                                        {
 
   // Loop through the panes and its real faces
   std::vector< COM::Pane*>::iterator it = _panes.begin();
   Manifold::PM_iterator pm_it=_surf->pm_begin();
   
   const Vector_3 null_vec(0,0,0);
 
   // Loop through the panes and its real vertices
   for (int i=0, local_npanes = _panes.size(); i<local_npanes; ++i, ++it) { 
     COM::Pane *pane = *it;
 
     // Obtain pointers 
     Vector_3 *evs = reinterpret_cast<Vector_3*>
       ( pane->attribute(_eigvecs->id())->pointer());
     Vector_3 *lambdas = reinterpret_cast<Vector_3*>
       ( pane->attribute(_eigvalues->id())->pointer());
     Vector_3 *ds = disps_buf ? reinterpret_cast<Vector_3*>
       ( pane->attribute(disps_buf->id())->pointer()) : NULL;
     const Vector_3 *bs = reinterpret_cast<Vector_3*>
       ( pane->attribute(_boffset->id())->pointer());
     const Vector_3 *bms = reinterpret_cast<Vector_3*>
       ( pane->attribute(_bmedial->id())->pointer());
 
     // Feature information
     const char *tranks = reinterpret_cast<char*>
       ( pane->attribute(_tangranks->id())->pointer());
     Vector_3 *vcnts = vcenters ? reinterpret_cast<Vector_3*>
       ( pane->attribute(vcenters->id())->pointer()) : NULL;
     const char *val_bndry_nodes = reinterpret_cast<const char*>
       ( pane->attribute(_cnstr_bndry_nodes->id())->pointer());
 
     Vector_3 *As = reinterpret_cast<Vector_3*>
       ( pane->attribute(_As->id())->pointer());
     const int *cnstrs = _cnstr_nodes ? reinterpret_cast<int*>
       ( pane->attribute(_cnstr_nodes->id())->pointer()) : NULL;
 
     // Loop through all real nodes of the pane
     for ( int j=0, jn=pane->size_of_real_nodes(); j<jn; ++j) {
       // Skip vertices with zero right-hand side, such as those at
       // edge- and face-centers for quadratic elements
       if ( lambdas[j][0]==0) continue;
 
       Vector_3 *es_v = &evs[3*j];   // eigen-vectors of current vertex.
       int nrank = 3-tranks[j];
 
       Vector_3 V[2];
       int     ndirs = 0;
       if ( cnstrs) get_constraint_directions( cnstrs[j], ndirs, V);
       
       if ( ndirs==3 ) { // Fixed point
         lambdas[j] = es_v[0] = es_v[1] = es_v[2] = null_vec;
 
         if (ds) ds[j] = null_vec;
         if (vcnts) vcnts[j] = null_vec;
         continue;
       }
 
       const Vector_3 *A_v = &As[3*j];
       double eps = _eig_thres*lambdas[j][0];
       // Propagate along the constrained direction.
       double cos_fangle = 0.98481, sin_fangle = 0.17365; // (10 deg)
       
       if ( ds) {
         Vector_3 b(0, 0, 0); 
         for ( int k=0; k<nrank; ++k) {
           if ( lambdas[j][k]>eps)
             b += ds[j]*es_v[k]*es_v[k];
         }
         ds[j] = b;
       }
 
       if ( vcnts) { // Project tangent motion within null space
         Vector_3 t(0, 0, 0);
         for ( int k=1; k<3; ++k) {
           if ( k>=nrank || lambdas[j][k]<eps)
             // Note that vcnts may have component in domain space.
             t += vcnts[j]*es_v[k]*es_v[k];
         }
         vcnts[j] = t;
       }
       if ( ndirs == 0) continue;
 
       // Evaluate surface normal
       Vector_3 nrm_v;
       if ( ds && ds[j]!=null_vec) {
         nrm_v = ds[j]; nrm_v.normalize(); 
       }
       else {
         if ( ds) ds[j] = null_vec;
 
         nrm_v = bms[j]; 
         obtain_directions( es_v, lambdas[j], tranks[j], nrm_v); 
       }
 
       // Propagate along the constrained direction.
       if ( ndirs == 2) {
         // Split planar constraints
         Vector_3 nrm = Vector_3::cross_product( V[0],V[1]);
         double  det = std::abs(nrm*nrm_v);
 
         if ( det < cos_fangle || val_bndry_nodes[j]) {
           // Split into two directions.
           V[0] = nrm_v; V[0] -= V[0]*nrm*nrm; V[0].normalize();
           V[1] = Vector_3::cross_product( V[0], nrm);
         
           if (vcnts) vcnts[j] = vcnts[j]*V[1]*V[1];
         }
         else if (vcnts) { // Plane is tangent space
           vcnts[j] -= vcnts[j]*nrm*nrm;
 
           if ( ds) ds[j] = null_vec; 
           continue; 
         }
       }
       else if ( vcnts && !val_bndry_nodes[j] &&
                 ( tranks[j]==1 && std::abs(V[0]*es_v[2]) > cos_fangle ||
                   tranks[j]==2 && std::abs(V[0]*nrm_v) < sin_fangle)) {
         // ndir == 1
         vcnts[j] = vcnts[j]*V[0]*V[0];
         if ( ds) ds[j] = null_vec;
         continue;
       }
 
       if ( ds && !ds[j].is_null()) {
         double VAV=Vector_3(V[0]*A_v[0],V[0]*A_v[1],V[0]*A_v[2])*V[0];
         double Vb = V[0]*bs[j];
 
         if ( VAV<eps) 
           ds[j] = null_vec;
         else {
           ds[j] = (-Vb/VAV)*V[0];
           if ( ndirs==1 && vcnts) vcnts[j] = null_vec;
         }
       }
     }
   }
 
   if ( _conserv && !disps_buf) balance_mass();
 
   return true;
 }

Here is the call graph for this function:

Here is the caller graph for this function:

void obtain_directions	(	Vector_3	es[3],
		const Vector_3 &	lambdas,
		int	nrank,
		Vector_3 &	b_medial,
		Vector_3 *	b_offset = `NULL`
	)		const

protected

Definition at line 475 of file FaceOffset_3.C.

References _eig_thres, k, and Vector_3< Type >::normalize().

Referenced by compute_directions(), and obtain_constrained_directions().

 {
   // Use smaller threshold for acute angles.
   double eps = lambdas[0]*_eig_thres;
   
   if ( b_offset) {
     // Solve for x within primary space
     Vector_3 d_offset(0,0,0);
     for ( int k=0; k<nrank; ++k) {
       if ( lambdas[k]>eps) {
         d_offset += (es[k]* *b_offset/-lambdas[k])*es[k];
       }
     }
 
     *b_offset = d_offset;
   }
 
   // Solve for x within primary space
   Vector_3 d_medial(0,0,0);
   for ( int k=0; k<nrank; ++k) {
     if ( lambdas[k]>eps) {
       d_medial += (es[k] * b_medial/-lambdas[k])*es[k];
     }
   }
 
   b_medial = d_medial.normalize();
 }

Here is the call graph for this function:

Here is the caller graph for this function:

void obtain_face_offset	(	const Point_3 *	pnts,
		const double *	spd_ptr,
		const COM::Attribute *	spd,
		double	dt,
		const Vector_3 *	dirs,
		COM::Element_node_enumerator &	ene,
		Vector_3 &	ns_nz,
		Point_3 &	cnt,
		Vector_3 *	disps,
		Vector_3 *	ns
	)

protected

Definition at line 332 of file FaceOffset_3.C.

References NTS::abs(), COM_assertion_msg, Vector_3< Type >::cross_product(), d, i, Element_node_enumerator::id(), Attribute::is_elemental(), Attribute::is_nodal(), k, Vector_3< Type >::norm(), Vector_3< Type >::normalize(), Element_node_vectors_k_const< Value >::set(), Attribute::size_of_components(), and Element_node_enumerator::size_of_edges().

Referenced by compute_quadrics().

                                                    {
 
   const Vector_3 null_vec(0,0,0);
   int ne=ene.size_of_edges();
 
   // Evaluate current face normal 
   Element_node_vectors_k_const<Point_3> ps; ps.set( pnts, ene, 1);
   SURF::Generic_element_2 e(ne);
   Point_3 pnts_buf[9];
 
   // Compute new position for all vertices 
   int ncomp = spd?spd->size_of_components():0;
   COM_assertion_msg( !ncomp || ncomp==1 && spd->is_elemental() ||
                      ncomp==3 && spd->is_nodal(), 
                      "Speed must be scalar-elemental or 3-vector-nodal");
 
   // Copy coordinates into buffer.
   for ( int i=0; i<ne; ++i) pnts_buf[i] = ps[i];
 
   // If dirs is present, then compute the displacement
   if ( dt!=0 && ncomp==3) { // Add dt*speed to each point
     for ( int i=0; i<ne; ++i)
       pnts_buf[i] += disps[i] = dt*((const Vector_3*)spd_ptr)[ene[i]-1];
   }
   else if ( dt==0) {
     for ( int i=0; i<ne; ++i) disps[i] = null_vec;
   }
 
   // Offset face normal equal to current face normal 
   e.interpolate_to_center((Vector_3*)pnts_buf, (Vector_3*)&cnt);
 
   // Compute normal directions.
   Vector_2 nc(0.5,0.5);
   Vector_3 J[2];
   e.Jacobian( pnts_buf, nc, J);
   ns_nz = Vector_3::cross_product( J[0], J[1]).normalize();
 
   if ( ne==3) {
     ns[0] = ns[1] = ns[2] = ns_nz;
   }
   else {
     for ( int k=0; k<ne; ++k) {
       Vector_3 ts[] = { pnts_buf[(k+1)%ne]-pnts_buf[k], 
                         pnts_buf[(k+ne-1)%ne]-pnts_buf[k]};
       ns[k] = Vector_3::cross_product( ts[0], ts[1]);
 
       // If the angle is far from 180 or 0 degrees, then use triangle's normal
       double l = ns[k].norm();
       if ( l > std::abs(0.05*ts[0]*ts[1])) ns[k] /= l;
       else ns[k] = ns_nz;
     }
   }
 
   // Compute points for normal motion.
   if (dt!=0 && ncomp==1) {
     double l = dt*spd_ptr[ene.id()-1];
     Vector_3 d = l*ns_nz;
     (Vector_3&)cnt += l*ns_nz; 
 
     if ( ne==3) {
       for ( int k=0; k<ne; ++k) pnts_buf[k] += (disps[k] = d);
     }
     else {
       for ( int k=0; k<ne; ++k) pnts_buf[k] += (disps[k] = l * ns[k]);
     }
   }
 
   // If wave-frontal motion, recalculate the directions and displacements.
   if ( dirs) {
     for ( int i=0; i<ne; ++i) {
       int vindex= ene[i]-1;
       
       if ( dirs[vindex].is_null()) {
         disps[i] = null_vec;
         pnts_buf[i] = ps[i] + disps[i];
       }
       else {
         Vector_3 dir = dirs[vindex];
         dir.normalize();
 
         Vector_3 tng = ( pnts_buf[(i+ne-1)%ne] - pnts_buf[i]) +
           (pnts_buf[(i+1)%ne] - pnts_buf[i]);   
 
         if ( tng * dir < 0) { // Expansion
           double l = (disps[i] * ns[i]);
           if ( dir * ns[i]<0) l=-l;
 
           (disps[i] = dir) *= l;
 
           pnts_buf[i] = ps[i] + disps[i];
         }
       }
     }
 
     // Reevaluate normal directions.
     e.Jacobian( pnts_buf, nc, J);
     ns_nz = Vector_3::cross_product( J[0], J[1]).normalize();
 
     if ( ne==3) {
       ns[0] = ns[1] = ns[2] = ns_nz;
     }
     else {
       for ( int k=0; k<ne; ++k) {
         ns[k] = Vector_3::cross_product
           ( pnts_buf[(k+1)%ne]-pnts_buf[k],
             pnts_buf[(k+ne-1)%ne]-pnts_buf[k]).normalize();
       }
     }
   }
 }

Here is the call graph for this function:

Here is the caller graph for this function:

static Vector_2 proj	(	const Vector_3 &	v,
		const Vector_3 &	d1,
		const Vector_3 &	d2
	)

inlinestaticprotected

Definition at line 318 of file FaceOffset_3.h.

                                                                 {
     return Vector_2( v*d1, v*d2);
   }

void reclassify_ridge_vertices	(	const bool	upgrade_corners,
		const bool	upgrade_ridge,
		COM::Attribute *	neighbor_feas,
		COM::Attribute *	tr_attr,
		bool	to_filter
	)

protected

Definition at line 461 of file detect_features.C.

References Propagation_3::_buf, _edges, Propagation_3::_panes, Propagation_3::_surf, _tol_cos_osta, NTS::abs(), COM_assertion, COM_DOUBLE, COM_NC, Rocblas::copy_scalar(), j, nj, and Vector_3< Type >::normalize().

Referenced by filter_and_identify_ridge_edges().

                                            {
 
   COM::Attribute *vecsums_buf = NULL;
   if ( upgrade_corners) {
     vecsums_buf = _buf->new_attribute( "FO_vecsums_buf", 'n', COM_DOUBLE, 3, "");
     _buf->resize_array( vecsums_buf, 0);
   }
 
   int zero=0;
   Rocblas::copy_scalar( &zero, neighbor_feas);
 
   Vector_3 *pnts=NULL, *vss=NULL;
   // Loop through the panes and the candidate edges to count incident edges
   std::vector< COM::Pane*>::iterator it = _panes.begin();
   Manifold::PM_iterator pm_it=_surf->pm_begin();
   for ( int ii=0, local_npanes = _panes.size(); 
         ii<local_npanes; ++ii, ++it, ++pm_it) { 
     COM::Pane *pane = *it;
     char *nnf = reinterpret_cast<char*>
       ( pane->attribute(neighbor_feas->id())->pointer());
     if ( upgrade_corners) {
       pnts = reinterpret_cast<Vector_3*>
         (pane->attribute(COM_NC)->pointer());
       vss = reinterpret_cast<Vector_3*>
         (pane->attribute(vecsums_buf->id())->pointer());
     }
  
     std::set< Edge_ID> &eset_pn = _edges[ii];
     for ( std::set< Edge_ID>::const_iterator eit=eset_pn.begin(),
             eend = eset_pn.end();  eit!=eend; ++eit) {
       Edge_ID eid = *eit;  COM_assertion(!eid.is_border());
       Element_node_enumerator ene( pane, eid.eid());
       int lid = eid.lid(), vindex = ene[lid]-1;
       ++nnf[vindex];
 
       Edge_ID eid_opp = pm_it->get_opposite_real_edge( eid);
       bool is_border_edge = pm_it->is_physical_border_edge( eid_opp);
       if ( is_border_edge)
         ++nnf[ ene[(lid+1)%ene.size_of_edges()]-1];
 
       if ( upgrade_corners) {
         int windex = ene[(lid+1)%ene.size_of_edges()]-1;
         Vector_3 diff = pnts[ windex] - pnts[vindex];
         vss[vindex] += diff.normalize();
       
         if ( is_border_edge) vss[windex] -= diff;
       }
     }
   }
   _surf->reduce_on_shared_nodes( neighbor_feas, Manifold::OP_SUM);
 
   if ( upgrade_corners)
     _surf->reduce_on_shared_nodes( vecsums_buf, Manifold::OP_SUM);
 
   double tol_sqsin = 4*( 1 - _tol_cos_osta*_tol_cos_osta);
   int ncrns = 0; // Number of corners
   // Loop through all vertices to upgrade all the vertices incident on
   // one or more than two ridge edges.  
   it = _panes.begin(); 
   for ( int ii=0, local_npanes = _panes.size(); ii<local_npanes; ++ii, ++it) {
     COM::Pane *pane = *it;
     char *nnf = reinterpret_cast<char*>
       ( pane->attribute(neighbor_feas->id())->pointer());
     char *tranks = reinterpret_cast<char*>
       ( pane->attribute(tr_attr->id())->pointer());
     if ( upgrade_corners)
       vss = reinterpret_cast<Vector_3*>
         (pane->attribute(vecsums_buf->id())->pointer());
 
     for ( int j=0, nj=pane->size_of_real_nodes(); j<nj; ++j) {
       if ( upgrade_corners && ( nnf[j]>=3 || nnf[j] == 2 &&
                                 vss[j].squared_norm() >= tol_sqsin))
         tranks[j] = 0;   // Upgrade joints to corners
       else if ( nnf[j] == 0 && std::abs(tranks[j]) == 1)
         tranks[j] = 2;  // Downgrade ridge vertices without neighbors
       else if ( upgrade_ridge && tranks[j] == 2 && nnf[j] >= 2 || 
                 upgrade_corners && tranks[j] == -1)
         tranks[j] = 1; // Upgrade ridge vertices with two or more neighbors
 
       ncrns += !tranks[j];
     }
   }
 
   // Reduce tangranks to ensure shared nodes have the same value across panes
   _surf->reduce_on_shared_nodes( tr_attr, Manifold::OP_MIN);
   
   if ( upgrade_corners) {
     _buf->delete_attribute( vecsums_buf->name());
     _buf->init_done( false);
 
     if ( to_filter)
       std::cout << "Found " << ncrns << " corner vertices" << std::endl;
   }
 }

Here is the call graph for this function:

Here is the caller graph for this function:

int remove_obscure_curves ( const std::vector< ObscendSet > & obscends )

protected

Definition at line 780 of file detect_features.C.

References _edges, Propagation_3::_panes, _ridgeneighbors, Propagation_3::_surf, _tangranks, RidgeNeighbor::eid, i, and RidgeNeighbor::vid.

Referenced by filter_obscended_ridge().

                                                                 {
   // Loop through the panes and ridge edges
   Manifold::PM_iterator pm_it=_surf->pm_begin();
 
   // Remove quasi-obscure edges.
   int dropped = 0;
   for ( int i=0, local_npanes = _panes.size(); 
         i<local_npanes; ++i, ++pm_it) { 
     std::set< Edge_ID> &eset_pn = _edges[i];
     const ObscendSet &obscend_pn = obscends[i];
     COM::Pane *pane = _panes[i];
 
     char *tranks = reinterpret_cast<char*>
       ( pane->attribute(_tangranks->id())->pointer());
     const RidgeNeighbor *rdgngbs = reinterpret_cast<RidgeNeighbor*>
       ( pane->attribute(_ridgeneighbors->id())->pointer());
 
     std::cout << "There are " << obscend_pn.size() 
               << " obscure curves. " << std::endl;
 
     for ( ObscendSet::const_iterator rit=obscend_pn.begin(),
             rend=obscend_pn.end(); rit!=rend; ++rit) {
       dropped += eset_pn.erase( rit->first);
       eset_pn.erase( pm_it->get_opposite_real_edge( rit->first));
 
       int cur = rit->second.first, next = rit->second.second, start=cur;
       for ( ; tranks[next-1]; ) {
         int  index = (next-1)*2;
         Edge_ID eid;
 
         if ( rdgngbs[ index].vid == cur) 
         { cur = next; next = rdgngbs[index+1].vid; eid = rdgngbs[ index+1].eid; }
         else if ( rdgngbs[ index+1].vid == cur) 
         { cur = next; next = rdgngbs[index].vid; eid = rdgngbs[ index].eid; }
         else
         { cur = next; next = 0; }
         if ( next==0 || next==start) break;
 
         dropped += eset_pn.erase( eid);
         eset_pn.erase( pm_it->get_opposite_real_edge( eid));
       }
     }
   }
 
   if ( dropped)
     std::cerr << "Dropped " << dropped << " false candidate edges. " << std::endl;
 
   return dropped;
 }

Here is the caller graph for this function:

double rescale_displacements	(	COM::Attribute *	disps,
		COM::Attribute *	vcenters,
		int	depth = `0`
	)

protected

Reduce time step based on stability limit, and rescale displacements by the reduced time step and diffusing expansion.

Returns the relative time step.

Definition at line 654 of file FaceOffset_3.C.

References Propagation_3::_buf, _courant, _facenormals, Propagation_3::_panes, Propagation_3::_surf, _tangranks, Rocblas::add(), COM_assertion_msg, COM_NC, COMMPI_Comm_size(), COMMPI_Initialized(), Vector_3< Type >::cross_product(), i, j, k, MPI_LOR, MPI_MIN, Rocblas::mul_scalar(), Element_node_enumerator::next(), nj, Element_node_enumerator::size_of_edges(), and solve_quadratic_func().

Referenced by time_stepping().

                                                           {
   
   double alpha = HUGE_VAL;
 
   // Loop through the panes and its real faces
   std::vector< COM::Pane*>::iterator it = _panes.begin();
   Manifold::PM_iterator pm_it=_surf->pm_begin();
   for ( int i=0, local_npanes = _panes.size(); i<local_npanes; 
         ++i, ++it, ++pm_it) { 
     COM::Pane *pane = *it;
     
     const Point_3 *pnts = reinterpret_cast<Point_3*>
       (pane->attribute(COM_NC)->pointer());
     const Vector_3 *ds = reinterpret_cast<const Vector_3*>
       ( pane->attribute(disps->id())->pointer());    
     const Vector_3 *vcnts = reinterpret_cast<const Vector_3*>
       ( pane->attribute(vcenters->id())->pointer());
     const Vector_3 *fnrms = reinterpret_cast<Vector_3*>
       ( pane->attribute(_facenormals->id())->pointer());
     const char      *tranks = reinterpret_cast<char*>
       ( pane->attribute(_tangranks->id())->pointer());
 
     // Loop through real elements of the current pane
     Element_node_enumerator ene( pane, 1); 
 
     for ( int j=0, nj=pane->size_of_real_elements(); j<nj; ++j, ene.next()) {
       int ne = ene.size_of_edges();
       int uindex=ene[ne-1]-1, vindex=ene[0]-1, windex=0;
 
       // Loop through all vertices of current face.
       for ( int k=0; k<ne; ++k, uindex=vindex, vindex=windex ) {
         windex = ene[(k+1==ne)?0:k+1]-1;
 
         if ( k!=0 && ne!=4 && (tranks[uindex] == 2 || tranks[vindex] == 2)) 
           continue; 
 
         Vector_3 tngs[2] = { pnts[uindex]-pnts[vindex],
                              pnts[windex]-pnts[vindex]};
 
         Vector_3 disp_v  = ds[vindex]+vcnts[vindex];
         Vector_3 dsps[2] = { ds[uindex]+vcnts[uindex]-disp_v, 
                              ds[windex]+vcnts[windex]-disp_v };
 
         Vector_3 nrm = Vector_3::cross_product( tngs[0], tngs[1]);
         Vector_3 dt  = Vector_3::cross_product( dsps[0], tngs[1])+
           Vector_3::cross_product( tngs[0], dsps[1]);
         Vector_3 dd = Vector_3::cross_product( dsps[0], dsps[1]);
   
         double a = dd*nrm;
         double b = dt*nrm;
         double c = nrm*nrm;
 
         std::pair<double, double>  roots=solve_quadratic_func(a,b,c);
         if ( roots.first<alpha && roots.first>0) alpha = roots.first;
         if ( roots.second<alpha && roots.second>0) alpha = roots.second;
 
         // Solve for time-step constraints at edges.
         Edge_ID eid(j+1,k), eid_opp = pm_it->get_opposite_real_edge( eid);
         bool is_border_edge = pm_it->is_physical_border_edge( eid_opp);
 
         if ( tranks[uindex] != 2 && tranks[vindex] != 2 && !is_border_edge) {
           const Vector_3 &n1 = fnrms[j];
           const Vector_3 &n2 = eid_opp.is_border() ? 
             pm_it->get_bd_normal( eid_opp) : fnrms[eid_opp.eid()-1];
 
           nrm = n1+n2;
           double a = dd*nrm;
           double b = dt*nrm;
           double c = nrm*nrm;
 
           roots=solve_quadratic_func(a,b,c);
           if ( roots.first<alpha && roots.first>0) alpha = roots.first;
           if ( roots.second<alpha && roots.second>0) alpha = roots.second;
         }
       }
     }
   }
 
   // Reduce alpha on all processors.
   MPI_Comm comm = _buf->get_communicator();
   int  needs_rescale = alpha*_courant<=1;
 
   if ( COMMPI_Initialized() && COMMPI_Comm_size( comm)>0) {
     double local = alpha;
     MPI_Allreduce(&local, &alpha, 1, MPI_DOUBLE, MPI_MIN, comm);
 
     int local_nr = needs_rescale;
     MPI_Allreduce( &local_nr, &needs_rescale, 1, MPI_INT, MPI_LOR, comm);    
   }
 
   // We need to scale back displacement
   if ( _courant>0 && alpha*_courant<=1) {
     double scale = 0.9*_courant*alpha;
 
     Rocblas::mul_scalar( disps, &scale, disps);
     Rocblas::mul_scalar( vcenters, &scale, vcenters);
 
     const int maxdepth=6;
     if ( depth>maxdepth) {
       std::cerr << "Rocprop: Mesh too distorted. Stopping" << std::endl;
       COM_assertion_msg( depth<=maxdepth, "Too many sub-iterations.");
       abort();
     }
 
     // Recursively call the function to recompute the scale.
     return scale*rescale_displacements( disps, vcenters, depth+1);
   }
   else {
     // Add tangential components to normal components.
     Rocblas::add( disps, vcenters, disps);
     
     return 1.;
   }
 }

Here is the call graph for this function:

Here is the caller graph for this function:

void reset_default_parameters ( bool b = false )

Definition at line 44 of file FaceOffset_3.C.

References _conserv, _courant, _dir_thres_strong, _dir_thres_weak, _eig_thres, _eig_weak_lb, _eig_weak_ub, _feature_layer, _inv_sqrt_c, _nrm_dfsn, _smoother, _tol_ad, _tol_angle_strong, _tol_angle_weak, _tol_cos_osta, _tol_eangle, _tol_kangle, _tol_kstrong, _tol_turn_angle, _wf_expn, _wght_scheme, cos, E2N_ANGLE, E2N_AREA, pi, SMOOTHER_LAPLACIAN, sqrt(), Propagation_3::square(), and cimg_library::tan().

Referenced by FaceOffset_3().

                                                                   {
 
   // Two default feature modes:
   if ( feature_preserve) {
     // Feature-preserving mode: (Recommended for structured meshes)
     _wght_scheme = SURF::E2N_ANGLE;
     _tol_angle_weak = 18.*pi/180;
   }
   else {
     // Feature-smearing mode: (Recommended for unstructured meshes)
     _wght_scheme = SURF::E2N_AREA;
     _tol_angle_weak = 30.*pi/180.;
   }
 
   _dir_thres_weak = square(tan(0.5*_tol_angle_weak)); // 0.003;
   _tol_angle_strong = 65.*pi/180;
   _dir_thres_strong = 0.7;
 
   _tol_kstrong    = 5;
   _tol_kangle     = (49./180.)*pi;
   _tol_eangle     = (25/180.)*pi;
 
   _tol_turn_angle = pi/4.5;
   _tol_ad = pi/3;
   _tol_cos_osta   = std::cos(pi/9.);
 
   _nrm_dfsn = 1; _eig_thres = 0.003;
   _eig_weak_lb = 1.e-4; _eig_weak_ub = 0.025;
 
   _courant = 0.5; _inv_sqrt_c = 1/std::sqrt(_courant);
   _smoother = SMOOTHER_LAPLACIAN; _conserv = 0; 
   _wf_expn = true; _feature_layer=0;
 }

Here is the call graph for this function:

Here is the caller graph for this function:

void set_conserve ( int b )

inline

Set mass conservation.

Definition at line 125 of file FaceOffset_3.h.

References _conserv.

125 { _conserv = b; }

FaceOffset_3::_conserv

bool _conserv

Definition: FaceOffset_3.h:399

void set_courant_constant ( double c )

inline

Set the constant in front of CFL condition (0<=c<1)

Definition at line 83 of file FaceOffset_3.h.

References _courant, _inv_sqrt_c, COM_assertion, and sqrt().

84 { _courant = c; _inv_sqrt_c = 1/std::sqrt(c); COM_assertion(0<=c && c<1); }

COM_assertion

#define COM_assertion(EX)

Error checking utility similar to the assert macro of the C language.

Definition: roccom_assertion.h:62

sqrt

double sqrt(double d)

Definition: double.h:73

FaceOffset_3::_inv_sqrt_c

double _inv_sqrt_c

Definition: FaceOffset_3.h:379

FaceOffset_3::_courant

double _courant

Definition: FaceOffset_3.h:378

Here is the call graph for this function:

void set_eigen_threshold ( double eps )

inline

Set the threshold of eigenvalues for splitting primary and null spaces.

Definition at line 76 of file FaceOffset_3.h.

References _eig_thres, and COM_assertion.

77 { _eig_thres = eps; COM_assertion( eps>0 && eps<1); }

COM_assertion

#define COM_assertion(EX)

Error checking utility similar to the assert macro of the C language.

Definition: roccom_assertion.h:62

FaceOffset_3::_eig_thres

double _eig_thres

Definition: FaceOffset_3.h:394

void set_fangle_strong ( double psi )

inline

Set the threshold for weak-feature angle.

Definition at line 100 of file FaceOffset_3.h.

References _tol_angle_strong, and pi.

                                       { 
     _tol_angle_strong = psi * pi/180;
   }

void set_fangle_turn ( double psi )

inline

Set the threshold for weak-feature angle.

Definition at line 105 of file FaceOffset_3.h.

References _tol_cos_osta, cos, and pi.

                                     { 
     _tol_cos_osta = std::cos(psi * pi/360);
   }

void set_fangle_weak ( double psi )

inline

Set the threshold for weak-feature angle.

Definition at line 94 of file FaceOffset_3.h.

References _dir_thres_weak, _tol_angle_weak, pi, Propagation_3::square(), and cimg_library::tan().

                                     {
     _tol_angle_weak = psi * pi/180;
     _dir_thres_weak = square(tan(0.5*_tol_angle_weak));
   }

Here is the call graph for this function:

void set_feature_layer ( bool b )

inline

Set whether or not to use feature layer.

Definition at line 70 of file FaceOffset_3.h.

References _feature_layer.

70 { _feature_layer = b; }

FaceOffset_3::_feature_layer

bool _feature_layer

Definition: FaceOffset_3.h:401

void set_normal_diffusion ( char i )

inline

Set number of iterations for normal diffusion.

Definition at line 67 of file FaceOffset_3.h.

References _nrm_dfsn, and i.

67 { _nrm_dfsn = i; }

i

blockLoc i

Definition: read.cpp:79

FaceOffset_3::_nrm_dfsn

char _nrm_dfsn

Definition: FaceOffset_3.h:375

void set_smoother ( int smoother )

inline

Set the choice of mesh smoother.

Definition at line 122 of file FaceOffset_3.h.

References _smoother.

122 { _smoother = smoother; }

FaceOffset_3::_smoother

char _smoother

Definition: FaceOffset_3.h:376

void set_wavefrontal_expansion ( bool b )

inline

Set the wavefrontal switch.

Definition at line 61 of file FaceOffset_3.h.

References _wf_expn.

61 { _wf_expn = b; }

FaceOffset_3::_wf_expn

bool _wf_expn

Definition: FaceOffset_3.h:373

void set_weighting_scheme ( int w )

inline

Weighting scheme.

Definition at line 90 of file FaceOffset_3.h.

References _wght_scheme.

90 { _wght_scheme = w; }

FaceOffset_3::_wght_scheme

char _wght_scheme

Definition: FaceOffset_3.h:374

double sign ( double x )

inlineprotected

Definition at line 312 of file FaceOffset_3.h.

                          {
     if ( x==0) return 0;
     if ( x>0) return 1;
     return -1;
   }

static void solve	(	const FT &	a1,
		const FT &	a2,
		const FT &	b1,
		const FT &	b2,
		const FT &	c1,
		const FT &	c2,
		FT &	x,
		FT &	y
	)

inlinestaticprotected

Definition at line 271 of file FaceOffset_3.h.

References denom.

Referenced by solve().

   {
     FT denom = a1*b2-b1*a2;
     
     x = - (b1*c2-c1*b2)/denom;
     
     y = (a1*c2-c1*a2)/denom;
   }

Here is the caller graph for this function:

static void solve	(	const FT &	a1,
		const FT &	a2,
		const FT &	a3,
		const FT &	b1,
		const FT &	b2,
		const FT &	b3,
		const FT &	c1,
		const FT &	c2,
		const FT &	c3,
		const FT &	d1,
		const FT &	d2,
		const FT &	d3,
		FT &	x,
		FT &	y,
		FT &	z
	)

inlinestaticprotected

Definition at line 284 of file FaceOffset_3.h.

References denom.

   {
     FT denom = b2*c1*a3-b1*c2*a3+c3*b1*a2+b3*c2*a1-c1*b3*a2-b2*c3*a1;
     
     x = - (b2*c3*d1-b2*c1*d3+c1*b3*d2+b1*c2*d3-c3*b1*d2-b3*c2*d1)/denom;
     
     z = (b2*d1*a3-b2*a1*d3+b1*a2*d3-b1*d2*a3-d1*b3*a2+a1*b3*d2)/denom;
     
     y = (a2*c3*d1-a2*c1*d3-c2*d1*a3+c2*a1*d3+d2*c1*a3-d2*c3*a1)/denom;
   }

static void solve	(	const Vector_3	A[3],
		const Vector_3 &	q,
		Vector_3 &	x
	)

inlinestaticprotected

Definition at line 300 of file FaceOffset_3.h.

References solve().

                                    {
     solve( A[0][0], A[0][1], A[0][2], A[1][0], A[1][1], A[1][2], 
            A[2][0], A[2][1], A[2][2], q[0], q[1], q[2], 
            x[0], x[1], x[2]);
   }

Here is the call graph for this function:

std::pair< double, double > solve_quadratic_func	(	double	a,
		double	b,
		double	c
	)

protected

Definition at line 772 of file FaceOffset_3.C.

References NTS::abs(), Mesquite::det(), max(), and sqrt().

Referenced by rescale_displacements().

                                                                 {
   double max_abc = std::max( std::max(std::abs(a), std::abs(b)), std::abs(c));
 
   if ( max_abc > 0) 
   { a /= max_abc; b /= max_abc; c /= max_abc; }
 
   if ( a == 0) {
     if ( b == 0) {
       return std::make_pair(HUGE_VAL, HUGE_VAL);
     }
 
     return std::make_pair(-c/b, HUGE_VAL);   
   }
   else {
     double det = b*b - 4*a*c;
     if ( det < 0) { 
       return std::make_pair(HUGE_VAL, HUGE_VAL);
     }
     
     double sqrt_d = std::sqrt( det);
 
     if ( b>0) {
       double temp = -b - sqrt_d;
       return std::make_pair(2*c/temp, 0.5*temp/a);
     }
     else {
       double temp = -b + sqrt_d;
       return std::make_pair( 0.5*temp/a, 2*c/temp);
     }
   }
 }

Here is the call graph for this function:

Here is the caller graph for this function:

double time_stepping	(	const COM::Attribute *	spds,
		double	dt,
		COM::Attribute *	disps,
		int *	smoothed = `NULL`
	)

virtual

Main entry of the algorithm.

At input, the value of smoothed specifies whether smoothing is desired.

Implements Propagation_3.

Definition at line 160 of file FaceOffset_3.C.

References Propagation_3::_bnd_set, _boffset, Propagation_3::_buf, Propagation_3::_cnstr_nodes, Propagation_3::_cnstr_set, _courant, _eigvalues, _eigvecs, _faceareas, _facecenters, _facenormals, _is_strtd, Propagation_3::_rank, _ridges, _scales, _smoother, Propagation_3::_surf, _tangranks, _vcenters, Propagation_3::_verb, _vnormals, _weights, _wf_expn, _wght_scheme, Rocblas::add(), Propagation_3::bound_nodal_motion(), COM_assertion_msg, compute_anisotropic_vertex_centers(), compute_directions(), compute_quadrics(), Rocblas::copy(), Rocblas::copy_scalar(), denoise_surface(), E2N_ANGLE, filter_and_identify_ridge_edges(), i, mark_weak_vertices(), Rocblas::mul_scalar(), nulaplacian_smooth(), obtain_constrained_directions(), rescale_displacements(), SMOOTHER_ANISOTROPIC, SMOOTHER_LAPLACIAN, and update_vertex_centers().

                                                    {
 
   COM_assertion_msg( spds || dt==0, 
                      "If no speed function, then time step must be zero.");
 
   // Inherit speeds and displacement onto the buffer window
   const COM::Attribute *spds_buf = spds ?
     ( (spds->window() == _buf) ? spds : 
       _buf->inherit( const_cast<COM::Attribute*>(spds), "speeds_inherited", 
                      COM::Pane::INHERIT_USE, true, NULL, 0)) : NULL;
   COM::Attribute *disps_buf = 
     _buf->inherit( disps, "disps_inherited", COM::Pane::INHERIT_USE, 
                    true, NULL, 0);
   COM::Attribute *disps_tmp =
     _buf->inherit( disps, "disps_inherited2", COM::Pane::INHERIT_CLONE, 
                    true, NULL, 0);
   _buf->init_done( false);
 
   // Compute the face areas and store them.
   SURF::Rocsurf::compute_element_areas( _faceareas);
 
   // First, compute the quadrics to determine the offset intersection 
   // with given time step. 
   bool with_nodal_velo = dt>0 && spds->is_nodal() && spds->size_of_components()==3;
 
   // If wavefrontal, then perform two setps.
   for ( int i=0; i<1+(_wf_expn && dt>0); ++i) {
     compute_quadrics( dt, spds_buf, _boffset, i ? disps_buf : NULL);
 
     if ( with_nodal_velo)
       Rocblas::mul_scalar( spds, &dt, disps_buf);
     else
       Rocblas::copy( _boffset, disps_buf);
 
     // Compute mean normals and offset directions during first iteration
     if ( i==0 && _wf_expn && dt>0) {
       compute_directions( !with_nodal_velo ? disps_buf : NULL, false);
 
       // Enforce constraints.
       obtain_constrained_directions( disps_buf, NULL);
     }
   }
   
   // Filter out isolated ridge vertices and identify ridge edges.
   filter_and_identify_ridge_edges(!_is_strtd && _wght_scheme==SURF::E2N_ANGLE && dt==0);
   mark_weak_vertices();
 
   // Recompute mean normals and offset directions with only primary component
   compute_directions( (dt!=0 && !with_nodal_velo) ? disps_buf : NULL, true);
 
   if ( dt==0)
     denoise_surface( NULL, disps_buf);
   else {
     double zero=0.;
     Rocblas::copy_scalar( &zero, disps_tmp);
     denoise_surface( disps_buf, disps_tmp);
     Rocblas::add( disps_buf, disps_tmp, disps_buf);
   }
 
   bool needs_smoothing=(!smoothed || *smoothed);
   if ( smoothed) *smoothed = true;
 
   if ( !needs_smoothing) {
     Rocblas::copy( disps_buf, _vcenters);
     if ( smoothed && *smoothed) *smoothed=false;
   }
   else if ( _is_strtd) {
     update_vertex_centers();  // Update vertex centers for ridge vertices
     nulaplacian_smooth( _vnormals, _tangranks, disps_buf,
                         _vcenters, disps_tmp, _weights);
   }
   else if ( _smoother == SMOOTHER_ANISOTROPIC) {
     // Compute anisotropic vertex centers.
     compute_anisotropic_vertex_centers( disps_buf);
   }
   update_vertex_centers();  // Update vertex centers for ridge vertices
 
   // Recompute face normals if propagation under nodal velocity
   if ( with_nodal_velo) {
     _surf->compute_normals( _facenormals);
     _surf->update_bd_normals( _facenormals, false);
   }
 
   // Constrain the displacements and vertex-centers.
   obtain_constrained_directions( disps_buf, _vcenters);
   if (_bnd_set) {
     bound_nodal_motion( disps_buf);
     _surf->reduce_on_shared_nodes( disps_buf, Manifold::OP_MAXABS);
   }
 
 #if EXPORT_DEBUG_INFO
   COM::Attribute *ctypes = disps->window()->attribute( "cnstr_types_nodal");
   if ( _cnstr_set && ctypes) Rocblas::copy( _cnstr_nodes, ctypes);
 
   COM::Attribute *tranks = disps->window()->attribute( "tangranks");
   if ( tranks) Rocblas::copy( _tangranks, tranks);
 
   COM::Attribute *facenormals = disps->window()->attribute( "facenormals");
   if ( facenormals) Rocblas::copy( _facenormals, facenormals);
   
   COM::Attribute *facecenters = disps->window()->attribute( "facecenters");
   if ( facecenters) Rocblas::copy( _facecenters, facecenters);
     
   COM::Attribute *eigvecs = disps->window()->attribute( "eigvecs");
   if ( eigvecs) Rocblas::copy( _eigvecs, eigvecs);
   
   COM::Attribute *lambdas = disps->window()->attribute( "lambdas");
   if ( lambdas) Rocblas::copy( _eigvalues, lambdas);
 
   COM::Attribute *ridges = disps->window()->attribute( "ridges");
   if ( ridges) 
     disps->window()->inherit( _ridges, "ridges", COM::Pane::INHERIT_CLONE, 
                               true, NULL, 0);
  
   if ( dt>0) {
     COM::Attribute *disps_novis = disps->window()->attribute( "disps_novis");
     if ( disps_novis) Rocblas::copy( disps, disps_novis);
   }
 #endif
 
   // Second, reduce time step based on stability limit, and rescale
   // displacements by the reduced time step and wavefrontal expansion.
   double dt_sub = dt; 
   if ( _courant>0) {
     for ( ;;) {
       double scale = rescale_displacements( disps_buf, _vcenters);
       
       if ( scale==1) break;
 
       if ( dt==0 && _verb && _rank==0)
         std::cout << "Rocprop: Scaled back displacement with a factor of " 
                   << scale << std::endl;
 
       dt_sub *= scale;
       
       if ( !with_nodal_velo || _smoother == SMOOTHER_LAPLACIAN) break;
       // If nodal velocity
       Rocblas::mul_scalar( spds, &dt_sub, disps_buf);
       if ( smoothed) smoothed=false;
     }
   }
   else if ( dt==0) {
     Rocblas::copy( _vcenters, disps_buf);
   }
 
 #if EXPORT_DEBUG_INFO
   if ( dt>0) {
     COM::Attribute *scales = disps->window()->attribute( "scales");
     if ( scales) Rocblas::copy( _scales, scales);
   }
 #endif
 
   _surf->reduce_on_shared_nodes( disps_buf, Manifold::OP_MAXABS);
 
   // Deallocate spds_buf and disps_buf
   _buf->delete_attribute( disps_tmp->name());
   _buf->delete_attribute( disps_buf->name());
   if ( spds_buf && spds_buf!=spds)
     _buf->delete_attribute( spds_buf->name());
   _buf->init_done( false);
 
   // Finally, return the current time step.
   return dt_sub;
 }

Here is the call graph for this function:

void update_face_volumes	(	const COM::Attribute *	fnormal,
		const COM::Attribute *	vdsps,
		COM::Attribute *	fvol
	)

protected

void update_vertex_centers ( )

protected

Definition at line 267 of file quadric_analysis.C.

References Propagation_3::_buf, Propagation_3::_cnstr_bndry_nodes, _ctangranks, _edges, Propagation_3::_panes, Propagation_3::_surf, _tangranks, _vcenters, COM_DOUBLE, COM_NC, i, j, and Element_node_enumerator::size_of_edges().

Referenced by time_stepping().

                         {
   COM::Attribute *vcenters_buf =
     _buf->new_attribute( "FO_vcenters_buf", 'n', COM_DOUBLE, 3, "");
   COM::Attribute *vecdiff_buf =
     _buf->new_attribute( "FO_vecdiff_buf", 'n', COM_DOUBLE, 3, "");
   _buf->resize_array( vcenters_buf, 0);
   _buf->resize_array( vecdiff_buf, 0);
 
   // Loop through the panes and ridge edges
   std::vector< COM::Pane*>::iterator it = _panes.begin();
   Manifold::PM_iterator pm_it=_surf->pm_begin();
 
   for ( int i=0, local_npanes = _panes.size(); i<local_npanes; 
         ++i, ++it, ++pm_it) { 
     COM::Pane *pane = *it;
     const std::set< Edge_ID> &eset_pn = _edges[i];  // List of ridge edges
 
     const Point_3 *pnts = reinterpret_cast<Point_3*>
       (pane->attribute(COM_NC)->pointer());
     Vector_3 *vcs = reinterpret_cast<Vector_3*>
       (pane->attribute(vcenters_buf->id())->pointer());
     Vector_3 *vdiff = reinterpret_cast<Vector_3*>
       (pane->attribute(vecdiff_buf->id())->pointer());
 
     for ( std::set< Edge_ID>::const_iterator eit=eset_pn.begin(),
             eend=eset_pn.end(); eit!=eend; ++eit) {
       // Make sure that eid_opp is not border
       Edge_ID eid = *eit, eid_opp = pm_it->get_opposite_real_edge( eid);
       bool is_border_edge = pm_it->is_physical_border_edge( eid_opp);
 
       Element_node_enumerator ene( pane, eid.eid());
       int lid = eid.lid(), neighbor = (lid+1)%ene.size_of_edges();
       int vindex = ene[lid]-1, windex = ene[neighbor]-1;
       
       // For ridge and flat vertices, compute based on edge centers
       Vector_3 diff = pnts[ windex] - pnts[vindex];
       vcs[vindex] += diff;
       if ( vdiff[vindex].is_null()) vdiff[vindex] = diff; 
       else vdiff[vindex] -= diff;
 
       if ( is_border_edge) {
         vcs[windex] -= diff;
         if ( vdiff[windex].is_null()) vdiff[windex] = -diff; 
         else vdiff[windex] += diff;
       }
     }
   }
   
   _surf->reduce_on_shared_nodes( vcenters_buf, Manifold::OP_SUM);
   _surf->reduce_on_shared_nodes( vecdiff_buf, Manifold::OP_DIFF);
 
   const Vector_3 null_vec(0,0,0);
   it = _panes.begin();
   // Finally, copy from vcenters_buf to vcenters
   for ( int i=0, local_npanes = _panes.size(); 
         i<local_npanes; ++i, ++it) {
     COM::Pane *pane = *it;
 
     Vector_3 *vcs_buf = reinterpret_cast<Vector_3*>
       (pane->attribute(vcenters_buf->id())->pointer());
     Vector_3 *vcs = reinterpret_cast<Vector_3*>
       (pane->attribute(_vcenters->id())->pointer());
     Vector_3 *vdiff = reinterpret_cast<Vector_3*>
       (pane->attribute(vecdiff_buf->id())->pointer());
     const char *tranks = reinterpret_cast<const char*>
       ( pane->attribute(_tangranks->id())->pointer());
     const char *cranks = reinterpret_cast<const char*>
       ( pane->attribute(_ctangranks->id())->pointer());
     const char *val_bndry_nodes = reinterpret_cast<const char*>
       ( pane->attribute(_cnstr_bndry_nodes->id())->pointer());
     
     for (int j=0, jn=pane->size_of_real_nodes(); j<jn; ++j) if (cranks[j]!=2) {
       if ( cranks[j]==0)
         vcs[ j] = null_vec;
       else if ( val_bndry_nodes && val_bndry_nodes[j] && !vdiff[j].is_null()) {
         vcs_buf[j] = (vcs_buf[j]*vdiff[j]*vdiff[j])/vdiff[j].squared_norm();
 
         vcs[ j] = (1./4.)*vcs_buf[j];
       }
     }
   }
 
   _buf->delete_attribute( vecdiff_buf->name());
   _buf->delete_attribute( vcenters_buf->name());
   _buf->init_done( false);
 }