Rocmop Class Reference

A Roccom mesh optimization module. More...

#include <Rocmop.h>

Inheritance diagram for Rocmop:

Collaboration diagram for Rocmop:

Public Types
enum	{ SMOOTH_VOL_MESQ_WG = 0, SMOOTH_VOL_MESQ_NG, SMOOTH_SURF_MEDIAL, SMOOTH_NONE }
	Valid Smoothing Methods. More...
enum	{ SMOOTH_VOL_MESQ_WG = 0, SMOOTH_VOL_MESQ_NG, SMOOTH_SURF_MEDIAL, SMOOTH_NONE }
	Valid Smoothing Methods. More...
enum	{ SMOOTH_VOL_MESQ_WG = 0, SMOOTH_VOL_MESQ_NG, SMOOTH_SURF_MEDIAL, SMOOTH_NONE }
	Valid Smoothing Methods. More...
enum	{ WRAPPER_SHAPE = 0, WRAPPER_BOEING, WRAPPER_MAX }
	Valid Wrapper Choices. More...

Public Member Functions
void	read_config_file (const std::string &)
void	add_aspect_ratios (COM::Attribute *usr_att, COM::Attribute *buf_att=NULL)
void	read_config_file (const std::string &)
void	add_aspect_ratios (COM::Attribute *usr_att, COM::Attribute *buf_att=NULL)
Contructors and Destructors
	Rocmop ()
	Default Constructor. More...
virtual	~Rocmop ()
	Destructor. More...
	Rocmop ()
	Default Constructor. More...
virtual	~Rocmop ()
	Destructor. More...
	Rocmop ()
	Default Constructor. More...
virtual	~Rocmop ()
	Destructor. More...

Static Public Member Functions
Load and Unload
static void	load (const std::string &mname)
	Loads Rocmop onto Roccom with a given module name. More...
static void	unload (const std::string &mname)
	Unloads Rocmop from Roccom. More...
static void	load (const std::string &mname)
	Loads Rocmop onto Roccom with a given module name. More...
static void	unload (const std::string &mname)
	Unloads Rocmop from Roccom. More...
static void	load (const std::string &mname)
	Loads Rocmop onto Roccom with a given module name. More...
static void	unload (const std::string &mname)
	Unloads Rocmop from Roccom. More...

Protected Types
enum	{ MOP_COOKIE =762667 }
enum	{ MOP_COOKIE =762667 }
enum	{ MOP_COOKIE =762667 }

Protected Member Functions
void	perturb_stationary ()
	Randomly perturn stationary nodes on pane boundary, not on phys. surface. More...
bool	check_displacements (COM::Attribute *disp)
void	zero_displacements (COM::Attribute *disp)
void	perturb_stationary ()
	Randomly perturn stationary nodes on pane boundary, not on phys. surface. More...
bool	check_displacements (COM::Attribute *disp)
void	zero_displacements (COM::Attribute *disp)
Smoothing Methods
void	smooth_vol_mesq_wg (double pre_quality)
	Smooth a volume via Mesquite using ghost information. More...
void	smooth_vol_mesq_ng (double pre_quality)
	Smooths a volume using Mesquite with only shared node information. More...
void	smooth_surf_medial ()
	Smooths a surface using the medial quadric. More...
void	smooth_mesquite (std::vector< COM::Pane * > &allpanes, int ghost_level=0)
	Smooth the panes of the working window using MESQUITE. More...
void	smooth_vol_mesq_wg ()
	Smooth a volume via Mesquite using ghost information. More...
void	smooth_vol_mesq_ng ()
	Smooths a volume using Mesquite with only shared node information. More...
void	smooth_surf_medial ()
	Smooths a surface using the medial quadric. More...
void	smooth_mesquite (std::vector< COM::Pane * > &allpanes, int ghost_level=0)
	Smooth the panes of the buffer window using MESQUITE. More...
void	smooth_vol_mesq_wg ()
	Smooth a volume via Mesquite using ghost information. More...
void	smooth_vol_mesq_ng ()
	Smooths a volume using Mesquite with only shared node information. More...
void	smooth_surf_medial ()
	Smooths a surface using the medial quadric. More...
void	smooth_boeing (COM::Attribute *att, int *niter)
void	smooth_mesquite (std::vector< COM::Pane * > &allpanes, int ghost_level=0)
	Smooth the panes of the buffer window using MESQUITE. More...
Quality Checking
double	check_marked_elem_quality (std::vector< std::vector< bool > > &marked_elems, std::vector< COM::Pane * > &allpanes)
	Get the largest dihedral angle of marked real elements. More...
double	check_all_elem_quality (std::vector< const COM::Pane * > &allpanes, bool with_ghost=false)
	Get the largest dihedral angle of all real elements. More...
double	check_all_elem_quality (std::vector< COM::Pane * > &allpanes, bool with_ghosts=false)
	Get the largest dihedral angle of all real elements marked. More...
void	print_quality (std::string &s)
	Print the quality range of all elements, for debugging. More...
void	print_mquality (std::string &s, std::vector< std::vector< bool > > &to_check)
	Print the quality range of marked elements, for debugging. More...
double	check_marked_elem_quality (std::vector< std::vector< bool > > &marked_elems, std::vector< COM::Pane * > &allpanes)
	Get the largest dihedral angle of marked real elements. More...
double	check_all_elem_quality (std::vector< const COM::Pane * > &allpanes, bool with_ghost=false)
	Get the largest dihedral angle of all real elements. More...
double	check_all_elem_quality (std::vector< COM::Pane * > &allpanes, bool with_ghosts=false)
	Get the largest dihedral angle of all real elements marked. More...
void	print_legible (int verb, const char *msg)
	Single process print message if given verbosity level is exceeded. More...
void	constrain_displacements (COM::Attribute *w_disp)
	Contrain displacements to _maxdisp. More...
void	print_extremal_dihedrals (COM::Window *window)
	Print the min and max dihedral angles along with their locations. More...
void	print_quality (std::string &s, const std::string &s2)
	Print the quality range of all elements, for debugging. More...
void	print_mquality (std::string &s, std::vector< std::vector< bool > > &to_check)
	Print the quality range of marked elements, for debugging. More...
double	check_marked_elem_quality (std::vector< std::vector< bool > > &marked_elems, std::vector< COM::Pane * > &allpanes)
	Get the largest dihedral angle of marked real elements. More...
double	check_all_elem_quality (std::vector< const COM::Pane * > &allpanes, bool with_ghost=false)
	Get the largest dihedral angle of all real elements. More...
double	check_all_elem_quality (std::vector< COM::Pane * > &allpanes, bool with_ghosts=false)
	Get the largest dihedral angle of all real elements marked. More...
void	print_legible (int verb, const char *msg)
	Single process print message if given verbosity level is exceeded. More...
void	constrain_displacements (COM::Attribute *w_disp)
	Contrain displacements to _maxdisp. More...
void	obtain_extremal_dihedrals (const COM::Attribute *att, double *min, double *max)
	Obtain the min and max dihedral angles. More...
void	print_extremal_dihedrals (COM::Window *window)
	Print the min and max dihedral angles along with their locations. More...
void	print_quality (std::string &s, const std::string &s2)
	Print the quality range of all elements, for debugging. More...
void	print_mquality (std::string &s, std::vector< std::vector< bool > > &to_check)
	Print the quality range of marked elements, for debugging. More...

Protected Attributes
bool	_reorthog
	Reorthogonalize? More...
char	_wght_scheme
	Weighting scheme. More...
double	_eig_thres
	Eigenvalue thresholds. More...
double	_dir_thres
	Another threshold. More...
double	_saliency_crn
int	_rediter
	No.iterations for vertex redistribution. More...
const COM::Attribute *	_cnstr_types
	Stores types of nodal constraints. More...
COM::Attribute *	_cnstr_dirs
	Stores directions of nodal contraints. More...
std::vector< std::set< Edge_ID > >	_edges
	ridge edges More...
std::vector< std::vector< bool > >	_is_shared_node
	Is the node shared? More...
std::vector< std::vector< bool > >	_is_shared_elem
	Does the element contain shared nodes? More...
std::vector< std::vector< bool > >	_is_phys_bnd_node
	Is the node on the physical boundary? More...
std::vector< std::vector< bool > >	_is_phys_bnd_elem
	Does the element contain nodes on the phys. boundary? More...
std::vector< std::vector< bool > >	_is_pane_bnd_node
	Is the node on the pane boundary? More...
std::vector< std::vector< bool > >	_is_pane_bnd_elem
	Does the element contain nodes on the pane boundary? More...
const COM::Window *	_usr_window
	The user's window. More...
COM::Window *	_wrk_window
	The working window. More...
bool	_is_pmesh
	pmesh or mesh ? More...
int	_method
	Choice of smoothing method. More...
int	_verb
	Verbose level. More...
int	_lazy
	Check quality before smoothing? More...
float	_tol
	Smoother iterating tolerance. More...
int	_niter
	Maximum number of iterations for smoothing. More...
float	_ctol
	Subcycling tolerance. More...
int	_ncycle
	Max number of subcycles for convergence. More...
int	_monotone
	Impose non-decreasing quality? More...
int	_cookie
	For Roccom. More...
int	_invert
	If true (default false), then invert the mesh. More...
COM::Window *	_buf_window
	The buffer window. More...
std::vector < MAP::Pane_dual_connectivity * >	_dcs
float	_maxdisp
	Maximum displacement allowed. More...
int	_smoothfreq
int	N
	Smooth every _smoothfreq'th call. More...
double	disp_tally
	originally a static in smooth() More...
float	_disp_thresh
	originally a static in check_displacements More...
int	_wrapper
	Choice of Mesquite Smoothing Wrappers. More...
int	_invert_tets
	If true (default false), then invert tets. More...
int	_invert_hexes
	If true (default false), then invert hexes. More...

Smoothing control flow
void	smooth (const COM::Attribute *pmesh, COM::Attribute *disp)
	Smooth the mesh in a Rocmop managed buffer. More...
void	smooth_in_place (COM::Attribute *pmesh)
	Smooth a mesh in place.. More...
void	smooth (const COM::Attribute *pmesh, COM::Attribute *disp)
	Smooth the mesh in a Rocmop managed buffer. More...
void	smooth_in_place (COM::Attribute *pmesh)
	Smooth a mesh in place.. More...
void	smooth (const COM::Attribute *pmesh, COM::Attribute *disp, double *timestep=NULL)
	Smooth the mesh in a Rocmop managed buffer. More...
void	smoother_specific_init ()
	Perform smoother specific initialization. More...
void	perform_smoothing (double pre_quality)
	Perform smoothing on _wrk_window. More...
void	perform_iterative_smoothing (double pre_quality)
	Perform iterative smoothing. More...
void	perform_noniterative_smoothing ()
	Perform noniterative smoothing. More...
void	smoother_specific_init ()
	Perform smoother specific initialization. More...
void	perform_smoothing ()
	Perform smoothing on _buf_window. More...
void	perform_iterative_smoothing ()
	Perform iterative smoothing. More...
void	perform_noniterative_smoothing ()
	Perform noniterative smoothing. More...
void	smoother_specific_init ()
	Perform smoother specific initialization. More...
void	perform_smoothing ()
	Perform smoothing on _buf_window. More...
void	perform_iterative_smoothing ()
	Perform iterative smoothing. More...
void	perform_noniterative_smoothing ()
	Perform noniterative smoothing. More...
void	update_buf_real_nc ()
	Update real nodal coordinates of _buf_window from _usr_window. More...
int	check_input_pconn ()
	Check pconn block 3 of the input mesh. More...
void	get_usr_disp (const COM::Attribute *pmesh, COM::Attribute *disp, double *timestep)
	Get displacement in _usr_window based on nc in _buf_window. More...

Miscellaneous
void	set_value (const char *opt, const void *val)
	Set a Rocomp option. More...
void	determine_pane_border ()
	Determine which nodes and elements are on pane borders. More...
void	determine_physical_border (COM::Attribute *w_is_surface)
	Determine which nodes are on the physical border. More...
void	determine_physical_border ()
	Determine which nodes and elements are on the physical border. More...
void	determine_shared_border ()
	Determine which nodes and elements are shared. More...
void	mark_elems_from_nodes (std::vector< std::vector< bool > > &marked_nodes, std::vector< std::vector< bool > > &marked_elems)
	Mark the nodes which contain marked elems. More...
void	set_value (const char *opt, const void *val)
	Set a Rocomp option. More...
void	determine_pane_border ()
	Determine which nodes and elements are on pane borders. More...
void	determine_physical_border (COM::Attribute *w_is_surface)
	Determine which nodes are on the physical border. More...
void	determine_physical_border ()
	Determine which nodes and elements are on the physical border. More...
void	determine_shared_border ()
	Determine which nodes and elements are shared. More...
void	mark_elems_from_nodes (std::vector< std::vector< bool > > &marked_nodes, std::vector< std::vector< bool > > &marked_elems)
	Mark the nodes which contain marked elems. More...
void	set_value (const char *opt, const void *val)
	Set a Rocomp option. More...
void	determine_pane_border ()
	Determine which nodes and elements are on pane borders. More...
void	determine_physical_border (COM::Attribute *w_is_surface)
	Determine which nodes are on the physical border. More...
void	determine_physical_border ()
	Determine which nodes and elements are on the physical border. More...
void	determine_shared_border ()
	Determine which nodes and elements are shared. More...
void	mark_elems_from_nodes (std::vector< std::vector< bool > > &marked_nodes, std::vector< std::vector< bool > > &marked_elems)
	Mark the nodes which contain marked elems. More...
virtual SURF::Window_manifold_2 *	manifold ()
	Obtain a reference to the manifold. More...
int	validate_object () const
	Check that the object is valid. More...
void	invert_tets ()
	Repair inverted tets. More...
virtual SURF::Window_manifold_2 *	manifold ()
	Obtain a reference to the manifold. More...
int	validate_object () const
	Check that the object is valid. More...
void	invert_tets ()
	Repair inverted tets. More...
virtual SURF::Window_manifold_2 *	manifold ()
	Obtain a reference to the manifold. More...
int	validate_object () const
	Check that the object is valid. More...
void	invert_elements (int conn_type)
	Repair inverted tets or hexes. More...

Inter-Pane Communication
void	agree_int (int &val, MPI_Op op)
	Agree on an integer value across all panes. More...
void	agree_double (double &val, MPI_Op op)
	Agree on a double value across all panes. More...
void	agree_int (int &val, MPI_Op op)
	Agree on an integer value across all panes. More...
void	agree_double (double &val, MPI_Op op)
	Agree on a double value across all panes. More...
void	agree_int (int &val, MPI_Op op)
	Agree on an integer value across all panes. More...
void	agree_double (double &val, MPI_Op op)
	Agree on a double value across all panes. More...
static void	reduce_sum_on_shared_nodes (COM::Attribute *att)
	Perform a sum-reduction on the shared nodes for the given attribute. More...
static void	reduce_sum_on_shared_nodes (COM::Attribute *att)
	Perform a sum-reduction on the shared nodes for the given attribute. More...
static void	reduce_sum_on_shared_nodes (COM::Attribute *att)
	Perform a sum-reduction on the shared nodes for the given attribute. More...

Functions exclusive to medial quadric smoothing.
void	evaluate_face_normals ()
	Evaluate face normals (copied from FaceOffset_3.[hC]. More...
void	identify_ridge_edges ()
	Identify ridge edges. More...
void	redistribute_vertices_ridge ()
	Redistribute ridge vertices within their tangent spaces. More...
void	redistribute_vertices_smooth ()
	Redistribute smooth vertices within their tangent spaces. More...
void	compute_medial_quadric ()
	Compute medial quadric for every vertex. More...
int	eigen_analyze_vertex (Vector_3< double > A_io[3], Vector_3< double > &b_io, Vector_3< double > *nrm_nz, double angle_defect)
void	get_redist_safe_factor (COM::Attribute *c_attr, COM::Attribute *w_attr, int rank)
void	evaluate_face_normals ()
	Evaluate face normals (copied from FaceOffset_3.[hC]. More...
void	identify_ridge_edges ()
	Identify ridge edges. More...
void	redistribute_vertices_ridge ()
	Redistribute ridge vertices within their tangent spaces. More...
void	redistribute_vertices_smooth ()
	Redistribute smooth vertices within their tangent spaces. More...
void	compute_medial_quadric ()
	Compute medial quadric for every vertex. More...
int	eigen_analyze_vertex (Vector_3< double > A_io[3], Vector_3< double > &b_io, Vector_3< double > *nrm_nz, double angle_defect)
void	get_redist_safe_factor (COM::Attribute *c_attr, COM::Attribute *w_attr, int rank)
void	evaluate_face_normals ()
	Evaluate face normals (copied from FaceOffset_3.[hC]. More...
void	identify_ridge_edges ()
	Identify ridge edges. More...
void	redistribute_vertices_ridge ()
	Redistribute ridge vertices within their tangent spaces. More...
void	redistribute_vertices_smooth ()
	Redistribute smooth vertices within their tangent spaces. More...
void	compute_medial_quadric ()
	Compute medial quadric for every vertex. More...
int	eigen_analyze_vertex (Vector_3< double > A_io[3], Vector_3< double > &b_io, Vector_3< double > *nrm_nz, double angle_defect)
void	get_redist_safe_factor (COM::Attribute *c_attr, COM::Attribute *w_attr, int rank)
static void	get_constraint_directions (int type, const Vector_3< double > &dir, int &ndirs, Vector_3< double > dirs[2])
	Get orthonormals of the constraints. More...
static void	compute_eigenvectors (Vector_3< double > A[3], Vector_3< double > &lambdas)
static void	compute_eigenvectors (Vector_2< double > A[2], Vector_2< double > &lambdas)
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< double > A[3], const Vector_3< double > &q, Vector_3< double > &x)
static void	get_constraint_directions (int type, const Vector_3< double > &dir, int &ndirs, Vector_3< double > dirs[2])
	Get orthonormals of the constraints. More...
static void	compute_eigenvectors (Vector_3< double > A[3], Vector_3< double > &lambdas)
static void	compute_eigenvectors (Vector_2< double > A[2], Vector_2< double > &lambdas)
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< double > A[3], const Vector_3< double > &q, Vector_3< double > &x)
static void	get_constraint_directions (int type, const Vector_3< double > &dir, int &ndirs, Vector_3< double > dirs[2])
	Get orthonormals of the constraints. More...
static void	compute_eigenvectors (Vector_3< double > A[3], Vector_3< double > &lambdas)
static void	compute_eigenvectors (Vector_2< double > A[2], Vector_2< double > &lambdas)
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< double > A[3], const Vector_3< double > &q, Vector_3< double > &x)

Detailed Description

A Roccom mesh optimization module.

Definition at line 41 of file Rocmop.h.

Member Enumeration Documentation

anonymous enum

Valid Smoothing Methods.

Enumerator
SMOOTH_VOL_MESQ_WG
SMOOTH_VOL_MESQ_NG
SMOOTH_SURF_MEDIAL
SMOOTH_NONE

Definition at line 45 of file Rocmop.h.

       { SMOOTH_VOL_MESQ_WG = 0, SMOOTH_VOL_MESQ_NG, 
         SMOOTH_SURF_MEDIAL, SMOOTH_NONE};

anonymous enum

protected

Enumerator
MOP_COOKIE

Definition at line 445 of file Rocmop.h.

{ MOP_COOKIE=762667};

anonymous enum

Valid Smoothing Methods.

Enumerator
SMOOTH_VOL_MESQ_WG
SMOOTH_VOL_MESQ_NG
SMOOTH_SURF_MEDIAL
SMOOTH_NONE

Definition at line 46 of file Rocmop_1.h.

       { SMOOTH_VOL_MESQ_WG = 0, SMOOTH_VOL_MESQ_NG, 
         SMOOTH_SURF_MEDIAL, SMOOTH_NONE};

anonymous enum

protected

Enumerator
MOP_COOKIE

Definition at line 472 of file Rocmop_1.h.

{ MOP_COOKIE=762667};

anonymous enum

Valid Smoothing Methods.

Enumerator
SMOOTH_VOL_MESQ_WG
SMOOTH_VOL_MESQ_NG
SMOOTH_SURF_MEDIAL
SMOOTH_NONE

Definition at line 46 of file Rocmop_2.h.

       { SMOOTH_VOL_MESQ_WG = 0, SMOOTH_VOL_MESQ_NG, 
         SMOOTH_SURF_MEDIAL, SMOOTH_NONE};

anonymous enum

Valid Wrapper Choices.

Enumerator
WRAPPER_SHAPE
WRAPPER_BOEING
WRAPPER_MAX

Definition at line 50 of file Rocmop_2.h.

{ WRAPPER_SHAPE = 0, WRAPPER_BOEING, WRAPPER_MAX};

anonymous enum

protected

Enumerator
MOP_COOKIE

Definition at line 499 of file Rocmop_2.h.

{ MOP_COOKIE=762667};

Constructor & Destructor Documentation

Rocmop ( )

inline

Default Constructor.

Definition at line 55 of file Rocmop.h.

           : _reorthog(true), _wght_scheme(SURF::E2N_ANGLE),
    _eig_thres(1.e-12),_dir_thres(.1), _saliency_crn(.1),
    _rediter(1), _cnstr_types(NULL), _cnstr_dirs(NULL),
    _usr_window(NULL), _wrk_window(NULL), _is_pmesh(false),
    _method(SMOOTH_NONE), _verb(0), _lazy(0), _tol(0.0),
    _niter(1), _ctol(.99), _ncycle(1), _monotone(0), _cookie(MOP_COOKIE),
    _invert(0)
    {;}

~Rocmop ( )

virtual

Destructor.

Definition at line 76 of file Rocmop.C.

                { 
  if(_verb)
    cout << "Entering Rocmop::~Rocmop\n";

  if (_wrk_window) {
    delete _wrk_window; _wrk_window = NULL;
  }

  if(_verb > 1)
    cout << "Exiting Rocmop::~Rocmop\n";
}

Rocmop ( )

inline

Default Constructor.

Definition at line 56 of file Rocmop_1.h.

           : _reorthog(true), _wght_scheme(SURF::E2N_ANGLE),
             _eig_thres(1.e-12),_dir_thres(.1), _saliency_crn(.1),
             _rediter(1), _cnstr_types(NULL), _cnstr_dirs(NULL),
             _usr_window(NULL), _buf_window(NULL), _is_pmesh(false), 
             _method(SMOOTH_VOL_MESQ_WG), _verb(0), _lazy(0), _tol(165.0),
             _niter(1), _ctol(.99), _ncycle(1), 
             _cookie(MOP_COOKIE),_invert(0), _maxdisp(0.0), _smoothfreq(1),
             N(-9999), disp_tally(0.0), _disp_thresh(0.0)
  {;}

virtual ~Rocmop ( )

virtual

Destructor.

Rocmop ( )

inline

Default Constructor.

Definition at line 59 of file Rocmop_2.h.

           : _reorthog(true), _wght_scheme(SURF::E2N_ANGLE),
    _eig_thres(1.e-12),_dir_thres(.1), _saliency_crn(.1),
    _rediter(1), _cnstr_types(NULL), _cnstr_dirs(NULL),
    _usr_window(NULL), _buf_window(NULL), _is_pmesh(false), 
    _method(SMOOTH_VOL_MESQ_WG),_wrapper(WRAPPER_SHAPE), _verb(0), _lazy(0), _tol(165.0),
    _niter(1), _ctol(.99), _ncycle(1), 
    _cookie(MOP_COOKIE),_invert_tets(0), _invert_hexes(0),
    _maxdisp(0.0), _smoothfreq(1),
    _disp_thresh(0.0)
  {;}

virtual ~Rocmop ( )

virtual

Destructor.

Member Function Documentation

void add_aspect_ratios	(	COM::Attribute *	usr_att,
		COM::Attribute *	buf_att = NULL
	)

Definition at line 454 of file Rocmop_1.C.

References COM_assertion_msg, COM_DOUBLE, Angle_Metric_3::compute(), Rocblas::copy(), Aspect_Metric_3::getAspects(), i, Geo_Metric_Base_3::initialize(), j, ni, nj, and Rocmap::update_ghosts().

                                                     {
  
  COM::Window* usr_window = usr_att->window();
  COM::Window* buf_window = NULL;

  bool del_buffer = false;

  if(!buf_att){      
    del_buffer = true;
    // Create a buffer window by cloning from the user window
    std::string buf_name(usr_window->name()+"-Rocmop_add_aspect");
    buf_window = new COM::Window(buf_name, 
                                 usr_window->get_communicator());
    buf_window->inherit( usr_att, "", 
                         COM::Pane::INHERIT_CLONE, true, NULL,0);
    buf_window->init_done();
  }           
  else
    buf_window = _buf_window;

  std::vector<Pane*> allpanes;

  buf_window->panes(allpanes);

  // Create and resize window level buffer attributes
  COM::Attribute * w_buff_R =
    buf_window->new_attribute("Aspect Ratio:  R (circumradius)",'e',COM_DOUBLE,1,"");
  COM::Attribute * w_buff_r =
    buf_window->new_attribute("Aspect Ratio:  r (inradius)",'e',COM_DOUBLE,1,"");
  COM::Attribute * w_buff_l =
    buf_window->new_attribute("Aspect Ratio:  l (shortest edge)",'e',COM_DOUBLE,1,"");
  COM::Attribute * w_buff_rR =
    buf_window->new_attribute("Aspect Ratio:  3r/R",'e',COM_DOUBLE,1,"");
  COM::Attribute * w_buff_Rr =
    buf_window->new_attribute("Aspect Ratio:  R/r",'e',COM_DOUBLE,1,"");
  COM::Attribute * w_buff_Rl =
    buf_window->new_attribute("Aspect Ratio:  R/l",'e',COM_DOUBLE,1,"");
  COM::Attribute * w_buff_min =
    buf_window->new_attribute("Aspect Ratio:  min dih",'e',COM_DOUBLE,1,"");
  COM::Attribute * w_buff_max =
    buf_window->new_attribute("Aspect Ratio:  max dih",'e',COM_DOUBLE,1,"");
  
  buf_window->resize_array(w_buff_R,0);
  buf_window->resize_array(w_buff_r,0);
  buf_window->resize_array(w_buff_l,0);
  buf_window->resize_array(w_buff_rR,0);
  buf_window->resize_array(w_buff_Rr,0);
  buf_window->resize_array(w_buff_Rl,0);
  buf_window->resize_array(w_buff_min,0);
  buf_window->resize_array(w_buff_max,0);
  buf_window->init_done();

  // Obtain buffer attribute ids
  int R_id   = w_buff_R->id();
  int r_id   = w_buff_r->id();
  int l_id   = w_buff_l->id();
  int rR_id  = w_buff_rR->id();
  int Rr_id  = w_buff_Rr->id();
  int Rl_id  = w_buff_Rl->id();
  int min_id = w_buff_min->id();
  int max_id = w_buff_max->id();

   // Create and resize window level user attributes
  COM::Attribute * w_usr_R =
    usr_window->new_attribute("Aspect Ratio:  R",'e',COM_DOUBLE,1,"");
  COM::Attribute * w_usr_r =
    usr_window->new_attribute("Aspect Ratio:  r",'e',COM_DOUBLE,1,"");
  COM::Attribute * w_usr_l =
    usr_window->new_attribute("Aspect Ratio:  l",'e',COM_DOUBLE,1,"");
  COM::Attribute * w_usr_rR =
    usr_window->new_attribute("Aspect Ratio:  3r/R",'e',COM_DOUBLE,1,"");
  COM::Attribute * w_usr_Rr =
    usr_window->new_attribute("Aspect Ratio:  R/r",'e',COM_DOUBLE,1,"");
  COM::Attribute * w_usr_Rl =
    usr_window->new_attribute("Aspect Ratio:  R/l",'e',COM_DOUBLE,1,"");
  COM::Attribute * w_usr_min =
    usr_window->new_attribute("Aspect Ratio:  min dih",'e',COM_DOUBLE,1,"");
  COM::Attribute * w_usr_max =
    usr_window->new_attribute("Aspect Ratio:  max dih",'e',COM_DOUBLE,1,"");
  
  usr_window->resize_array(w_usr_R,0);
  usr_window->resize_array(w_usr_r,0);
  usr_window->resize_array(w_usr_l,0);
  usr_window->resize_array(w_usr_rR,0);
  usr_window->resize_array(w_usr_Rr,0);
  usr_window->resize_array(w_usr_Rl,0);
  usr_window->resize_array(w_usr_min,0);
  usr_window->resize_array(w_usr_max,0);
  usr_window->init_done();
    
  double R,r,l;
  double angles[] = {0.0,0.0};

  for(int i=0,ni=(int)allpanes.size();i<ni;++i){

    // get Pane level buffer attributes
    double *R_ptr = reinterpret_cast<double *>
      (allpanes[i]->attribute(R_id)->pointer());    
    double *r_ptr = reinterpret_cast<double *>
      (allpanes[i]->attribute(r_id)->pointer());    
    double *l_ptr = reinterpret_cast<double *>
      (allpanes[i]->attribute(l_id)->pointer());    
    double *rR_ptr = reinterpret_cast<double *>
      (allpanes[i]->attribute(rR_id)->pointer());    
    double *Rr_ptr = reinterpret_cast<double *>
      (allpanes[i]->attribute(Rr_id)->pointer());    
    double *Rl_ptr = reinterpret_cast<double *>
      (allpanes[i]->attribute(Rl_id)->pointer());    
    double *min_ptr = reinterpret_cast<double *>
      (allpanes[i]->attribute(min_id)->pointer());    
    double *max_ptr = reinterpret_cast<double *>
      (allpanes[i]->attribute(max_id)->pointer());    
    
    for(int j=0,nj=allpanes[i]->size_of_elements();j<nj;++j){
      
      Element_node_enumerator ene(allpanes[i],j+1);
      
      Aspect_Metric_3 aspect_met;       
      aspect_met.initialize(ene);
      aspect_met.getAspects(R,r,l);

      Angle_Metric_3 angle_met;
      angle_met.initialize(ene);
      angle_met.compute(angles);
     
      R_ptr[j] = R;
      r_ptr[j] = r;
      l_ptr[j] = l;
      COM_assertion_msg(R != 0.0,
                        "Cannot calculate aspect ratios with zero magnitude circumradius");
      rR_ptr[j] = (3.0*r)/R;
      COM_assertion_msg(r != 0.0,
                        "Cannot calculate aspect ratios with zero magnitude circumradius");
      Rr_ptr[j] = R/r;
      COM_assertion_msg(l != 0.0,
                        "Cannot calculate aspect ratios with zero magnitude circumradius");
      Rl_ptr[j] = R/l;

      min_ptr[j] = angles[0];
      max_ptr[j] = angles[1];
    }
  }

  Rocmap::update_ghosts(w_buff_R);
  Rocmap::update_ghosts(w_buff_r);
  Rocmap::update_ghosts(w_buff_l);
  Rocmap::update_ghosts(w_buff_rR);
  Rocmap::update_ghosts(w_buff_Rr);
  Rocmap::update_ghosts(w_buff_Rl);
  Rocmap::update_ghosts(w_buff_min);
  Rocmap::update_ghosts(w_buff_max);

  Rocblas::copy(w_buff_R,w_usr_R);
  Rocblas::copy(w_buff_r,w_usr_r);
  Rocblas::copy(w_buff_l,w_usr_l);
  Rocblas::copy(w_buff_rR,w_usr_rR);
  Rocblas::copy(w_buff_Rr,w_usr_Rr);
  Rocblas::copy(w_buff_Rl,w_usr_Rl);
  Rocblas::copy(w_buff_min,w_usr_min);
  Rocblas::copy(w_buff_max,w_usr_max);

  if(del_buffer)
    delete buf_window;
}

Here is the call graph for this function:

void add_aspect_ratios	(	COM::Attribute *	usr_att,
		COM::Attribute *	buf_att = NULL
	)

void agree_double ( double &  val, MPI_Op  op  )

inlineprotected

Agree on a double value across all panes.

Parameters

val	The local value.
op	The type of reduction to perform.

Definition at line 239 of file Rocmop_1.h.

References _usr_window, and COMMPI_Initialized().

                                           {
    if(COMMPI_Initialized()){
      double temp = val;
      MPI_Allreduce(&val, &temp, 1, MPI_DOUBLE, op, 
                    _usr_window->get_communicator());
      val = temp;
    }
  }

Here is the call graph for this function:

void agree_double ( double &  val, MPI_Op  op  )

inlineprotected

Agree on a double value across all panes.

Parameters

val	The local value.
op	The type of reduction to perform.

Definition at line 240 of file Rocmop.h.

References _usr_window, and COMMPI_Initialized().

                                           {
    if(COMMPI_Initialized()){
      double temp = val;
      MPI_Allreduce(&val, &temp, 1, MPI_DOUBLE, op, 
                    _usr_window->get_communicator());
      val = temp;
    }
  }

Here is the call graph for this function:

void agree_double ( double &  val, MPI_Op  op  )

inlineprotected

Agree on a double value across all panes.

Parameters

val	The local value.
op	The type of reduction to perform.

Definition at line 261 of file Rocmop_2.h.

References _usr_window, and COMMPI_Initialized().

                                           {
    if(COMMPI_Initialized()){
      double temp = val;
      MPI_Allreduce(&val, &temp, 1, MPI_DOUBLE, op, 
                    _usr_window->get_communicator());
      val = temp;
    }
  }

Here is the call graph for this function:

void agree_int ( int &  val, MPI_Op  op  )

inlineprotected

Agree on an integer value across all panes.

Parameters

val	The local value.
op	The type of reduction to perform.

Definition at line 225 of file Rocmop_1.h.

References _usr_window, and COMMPI_Initialized().

                                     {
    if(COMMPI_Initialized()){
      int temp = val;
      MPI_Allreduce(&val, &temp, 1, MPI_INT, op, 
                    _usr_window->get_communicator());
      val = temp;
    }
  }

Here is the call graph for this function:

void agree_int ( int &  val, MPI_Op  op  )

inlineprotected

Agree on an integer value across all panes.

Parameters

val	The local value.
op	The type of reduction to perform.

Definition at line 226 of file Rocmop.h.

References _usr_window, and COMMPI_Initialized().

                                     {
    if(COMMPI_Initialized()){
      int temp = val;
      MPI_Allreduce(&val, &temp, 1, MPI_INT, op, 
                    _usr_window->get_communicator());
      val = temp;
    }
  }

Here is the call graph for this function:

void agree_int ( int &  val, MPI_Op  op  )

inlineprotected

Agree on an integer value across all panes.

Parameters

val	The local value.
op	The type of reduction to perform.

Definition at line 247 of file Rocmop_2.h.

References _usr_window, and COMMPI_Initialized().

                                     {
    if(COMMPI_Initialized()){
      int temp = val;
      MPI_Allreduce(&val, &temp, 1, MPI_INT, op, 
                    _usr_window->get_communicator());
      val = temp;
    }
  }

Here is the call graph for this function:

double check_all_elem_quality ( std::vector< const COM::Pane * > &  allpanes, bool  with_ghost = false  )

protected

Get the largest dihedral angle of all real elements.

Parameters

allpanes

the set of panes containing the elements

Definition at line 785 of file Rocmop.C.

References Angle_Metric_3::compute(), i, Geo_Metric_Base_3::initialize(), k, ni, nk, and rank.

                                                       {
  if(_verb) 
    std::cout << "Exiting Rocmop::check_all_elem_quality" << std::endl;

  int rank =0;
  int ierr = MPI_Comm_rank( _wrk_window->get_communicator(),
                            &rank); assert( ierr == 0);
  double worst_angle = 0.0;
  double angles[] = {0.0, 0.0};
  for(int i=0,ni = allpanes.size(); i<ni; ++i){
    int nk=allpanes[i]->size_of_real_elements();
    if(with_ghosts)
      nk = allpanes[i]->size_of_elements();
    for(int k =0; k<nk; ++k){
      Element_node_enumerator ene(allpanes[i],k+1);
      Angle_Metric_3 am;
      am.initialize(ene);
      am.compute(angles);
      if(angles[1]>worst_angle)
        worst_angle = angles[1];
    }
  }
  return worst_angle;

  if(_verb > 1) 
    std::cout << "Exiting Rocmop::check_all_elem_quality" << std::endl;
}

Here is the call graph for this function:

double check_all_elem_quality ( std::vector< const COM::Pane * > &  allpanes, bool  with_ghost = false  )

protected

Get the largest dihedral angle of all real elements.

Parameters

allpanes

the set of panes containing the elements

double check_all_elem_quality ( std::vector< COM::Pane * > &  allpanes, bool  with_ghosts = false  )

protected

Get the largest dihedral angle of all real elements marked.

Parameters

allpanes

the set of panes containing the elements

Definition at line 814 of file Rocmop.C.

References Angle_Metric_3::compute(), i, Geo_Metric_Base_3::initialize(), k, ni, nk, and rank.

                                                       {
  if(_verb) 
    std::cout << "Entering Rocmop::check_all_elem_quality" << std::endl;

  int rank =0;
  int ierr = MPI_Comm_rank( _wrk_window->get_communicator(),
                            &rank); assert( ierr == 0);

  double worst_angle = 0.0;
  double angles[] = {0.0, 0.0};
  for(int i=0,ni = allpanes.size(); i<ni; ++i){
    int nk=allpanes[i]->size_of_real_elements();
    if(with_ghosts){
      nk = allpanes[i]->size_of_elements();
    }
    for(int k =0; k<nk; ++k){
      Element_node_enumerator ene(allpanes[i],k+1);
      Angle_Metric_3 am;
      am.initialize(ene);
      am.compute(angles);

      if(angles[1]>worst_angle)
        worst_angle = angles[1];
    }
  }
  return worst_angle;

  if(_verb > 1) 
    std::cout << "Exiting Rocmop::check_all_elem_quality" << std::endl;
}

Here is the call graph for this function:

double check_all_elem_quality ( std::vector< COM::Pane * > &  allpanes, bool  with_ghosts = false  )

protected

Get the largest dihedral angle of all real elements marked.

Parameters

allpanes

the set of panes containing the elements

double check_all_elem_quality ( std::vector< const COM::Pane * > &  allpanes, bool  with_ghost = false  )

protected

Get the largest dihedral angle of all real elements.

Parameters

allpanes

the set of panes containing the elements

double check_all_elem_quality ( std::vector< COM::Pane * > &  allpanes, bool  with_ghosts = false  )

protected

Get the largest dihedral angle of all real elements marked.

Parameters

allpanes

the set of panes containing the elements

bool check_displacements ( COM::Attribute *  disp )

protected

Definition at line 357 of file Rocmop_1.C.

References i, j, MPI_MAX, ni, nj, Vector_3< Type >::norm(), and nvc::norm().

{
#ifdef ROCPROF
  Profile_begin("Rocmop::check_disp");
#endif
  
  int disp_id = w_disp->id();
  double max_norm = 0.0;
  // move to class member function
  //static double disp_tally = 0.0;
  bool retval = true;
  //  std::ofstream Ouf;
  //  Ouf.open("displacements.txt",ios::app);
  std::vector<Pane*> allpanes;
  const_cast<COM::Window*>(_usr_window)->panes(allpanes);
    
  // Find the maximum displacement on this processor
  for(int i=0,ni = allpanes.size(); i<ni; ++i){
    
    Vector_3<double> *ptr = NULL;
    COM::Attribute* ptr_att = allpanes[i]->attribute(disp_id);
    void* void_ptr = ptr_att->pointer();
    ptr = reinterpret_cast<Vector_3<double>*>(void_ptr);

    for(int j=0,nj = allpanes[i]->size_of_real_nodes(); j<nj; ++j){     
      double norm = ptr[j].norm();
      // Re-assign them to zero to maintain original Rocmop behavior
      // just in case.  I'm pretty sure that any values in this array
      // already get wiped out when the Mesquite displacements are 
      // copied in.
      ptr[j][0] = 0.0;
      ptr[j][1] = 0.0;
      ptr[j][2] = 0.0;
      //      Ouf << j << ": " << norm << std::endl;
      if(norm > max_norm)
        max_norm = norm;      
    }
  }
  // Maximize the displacement across all processors
  agree_double(max_norm,MPI_MAX);

  // Update the tally
  disp_tally += max_norm;

  // Build a blurb about displacement
  std::ostringstream Ostr;
  Ostr << "Mesh Displacement: " << disp_tally;

  // Determine whether we have exceeded the threshold
  if(disp_tally  > _disp_thresh){
    disp_tally = 0.0;
    retval = true;
    // If we intend to smooth - also print out the displacement
    print_legible(0,Ostr.str().c_str());
  }
  else {
    retval = false;
    // If not smoothing, print out displacement if verbosity is higher
    print_legible(1,Ostr.str().c_str());
  }
#ifdef ROCPROF
  Profile_end("Rocmop::check_disp");
#endif
  //  Ouf.close();
  return(retval);
}

Here is the call graph for this function:

bool check_displacements ( COM::Attribute *  disp )

protected

int check_input_pconn ( )

protected

Check pconn block 3 of the input mesh.

Check if all ghost nodes are listed in pconn block 3 of _usr_window

Definition at line 886 of file Rocmop_2.C.

References COM_NC, COM_PCONN, COMMPI_Initialized(), i, j, k, ni, and rank.

{
  int rank =0;
  if(COMMPI_Initialized())
  {
    int ierr = MPI_Comm_rank(_usr_window->get_communicator(), &rank); 
    assert( ierr == 0);
  }

  std::vector<const Pane*> allusrpanes;
  _usr_window->panes(allusrpanes);

  for(int i = 0, ni = allusrpanes.size(); i < ni; ++i)
  {
    // Obtain pane level attributes
    const COM::Attribute *usr_nc = allusrpanes[i]->attribute(COM::COM_NC);
    const COM::Attribute *usr_pconn = allusrpanes[i]->attribute(COM::COM_PCONN);
    int minGhostNodeID = usr_nc->size_of_real_items() + 1;
    int maxGhostNodeID = usr_nc->size_of_items();
    int numGhostNodes  = maxGhostNodeID - minGhostNodeID + 1;
    int * pconnArray = (int *)usr_pconn->pointer();

    // Check pconn
    bool OORGhostNodes  = false;    // found out of range ghost node
    bool missGhostNodes = false;    // ghost nodes missed from pconn

    // Allocate marking array
    int * isMarked = new int [numGhostNodes];
    for (int j = 0; j < numGhostNodes; j++)
      isMarked[j] = 0;

    // Jump to pconn block 3
    int index = 0;
    int numNbrs = 0;
    int numItems = 0;

    for (int j = 0; j < 2; j++)
    {
      numNbrs = pconnArray[index];
      index++;

      for (int k = 0; k < numNbrs; k++)
      {
        index++;                         // step over proc id
        numItems = pconnArray[index];
        index++;
        index += numItems;
      }
    }

    // Begin marking
    numNbrs = pconnArray[index];
    index++;

    for (int j = 0; j < numNbrs; j++)
    {
      index++;                         // step over proc id
      numItems = pconnArray[index];
      index++;

      for (int k = 0; k < numItems; k++)
      {
        int lid = pconnArray[index];
        index++;

        if (lid < minGhostNodeID || lid > maxGhostNodeID)
        {
          std::cerr << "Rocmop: Error in check_input_pconn(): Rank "
                    << rank << ", Ghost node " << lid 
                    << " is out of ghost node range.\n" << std::endl;
          OORGhostNodes = true;
        }
        else if (isMarked[lid - minGhostNodeID] == 1)
        {
          if (_verb > 0)
            std::cerr << "Rocmop: Warning in check_input_pconn(): Rank "
                      << rank << ", Ghost node " << lid 
                      << " has > 1 owners.\n" << std::endl;
        }
        else
          isMarked[lid - minGhostNodeID] = 1;
      }
    }

    // Check if any node is not marked
    for (int j = 0; j < numGhostNodes; j++)
      if (isMarked[j] == 0)
      {
        missGhostNodes = true;
        std::cerr << "Rocmop: Error in check_input_pconn(): Rank "
                  << rank << ", Ghost node " << (j + minGhostNodeID)
                  << " is not listed in pconn block 3.\n" << std::endl;
      }

    // Free up memory
    delete isMarked;

    // Return the result
    if (OORGhostNodes == true || missGhostNodes == true)
      return 1;
  }

  return 0;
}

Here is the call graph for this function:

double check_marked_elem_quality ( std::vector< std::vector< bool > > &  marked_elems, std::vector< COM::Pane * > &  allpanes  )

protected

Get the largest dihedral angle of marked real elements.

Parameters

marked_elems	the elements whose quality is to be checked.
allpanes	the set of panes containing the elements

Definition at line 760 of file Rocmop.C.

References Angle_Metric_3::compute(), i, Geo_Metric_Base_3::initialize(), k, ni, and nk.

                                                                       {
  if(_verb) 
    std::cout << "Entering Rocmop::check_marked_elems" << std::endl;

  double worst_angle = 0.0;
  double angles[] = {0.0, 0.0};
  for(int i=0,ni = allpanes.size(); i<ni; ++i){
    for(int k =0,nk=allpanes[i]->size_of_real_elements(); k<nk; ++k){
      if(marked_elems[i][k]){
        Element_node_enumerator ene(allpanes[i],k+1);
        Angle_Metric_3 am;      
        am.initialize(ene);
        am.compute(angles);
        if(angles[1]>worst_angle)
          worst_angle = angles[1];
      }
    }
  }
  return worst_angle;

  if(_verb > 1) 
    std::cout << "Exiting Rocmop::check_marked_elems" << std::endl;
}

Here is the call graph for this function:

double check_marked_elem_quality ( std::vector< std::vector< bool > > &  marked_elems, std::vector< COM::Pane * > &  allpanes  )

protected

Get the largest dihedral angle of marked real elements.

Parameters

marked_elems	the elements whose quality is to be checked.
allpanes	the set of panes containing the elements

double check_marked_elem_quality ( std::vector< std::vector< bool > > &  marked_elems, std::vector< COM::Pane * > &  allpanes  )

protected

Get the largest dihedral angle of marked real elements.

Parameters

marked_elems	the elements whose quality is to be checked.
allpanes	the set of panes containing the elements

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

staticprotected

Definition at line 1191 of file smooth_medial.C.

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

Referenced by 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<double> *buf = (const Vector_3<double>*)&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:

static void compute_eigenvectors ( Vector_2< double >  A[2], Vector_2< double > &  lambdas  )

staticprotected

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

staticprotected

static void compute_eigenvectors ( Vector_2< double >  A[2], Vector_2< double > &  lambdas  )

staticprotected

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

staticprotected

static void compute_eigenvectors ( Vector_2< double >  A[2], Vector_2< double > &  lambdas  )

staticprotected

void compute_medial_quadric ( )

protected

Compute medial quadric for every vertex.

Definition at line 286 of file smooth_medial.C.

References _cnstr_dirs, _cnstr_types, _verb, _wght_scheme, _wrk_window, cimg_library::acos(), angle(), COM_NC, Rocblas::copy_scalar(), Vector_3< Type >::cross_product(), d, E2N_ANGLE, E2N_AREA, E2N_ONE, E2N_SPHERE, eigen_analyze_vertex(), half_pi, i, j, k, max(), ni, nj, Vector_3< Type >::norm(), Vector_3< Type >::normalize(), pi, REAL_PANE, rs(), s, sin, and sqrt().

Referenced by smooth_surf_medial().

                         {
 if(_verb > 1) 
   std::cout << "Entering Rocmop::compute_medial_quadric" << std::endl;
  // Use the following buffer spaces: _eigvecs for matrix A; 
  // _eigvalues to store bm; disps to store bo, _weights to store r
 const std::string att1("eigvecs"), att2("lambda"), att3("weights"),
    att4("facenormals"), att5("facecenters"), att6("tangranks"),
   att7("cntnranks"), att8("cntnvecs"), att9("vnormals"), att10("disps"),
   att11("awnormals"), att12("uwnormals");
  COM::Attribute *A_attr =  _wrk_window->attribute(att1);
  COM::Attribute *bm_attr = _wrk_window->attribute(att2);
  COM::Attribute *aw_attr = _wrk_window->attribute(att11);
  COM::Attribute *uw_attr = _wrk_window->attribute(att12);
  COM::Attribute *r_attr = _wrk_window->attribute(att3);

  int facenormals_id = _wrk_window->attribute(att4)->id();
  int facecenters_id = _wrk_window->attribute(att5)->id();
  int tangranks_id = _wrk_window->attribute(att6)->id();
  int cntnranks_id = _wrk_window->attribute(att7)->id();
  int cntnvecs_id = _wrk_window->attribute(att8)->id();
  int vnormals_id = _wrk_window->attribute(att9)->id();
  int awnormals_id = _wrk_window->attribute(att11)->id();
  int uwnormals_id = _wrk_window->attribute(att12)->id();
  int disps_id = _wrk_window->attribute(att10)->id();

  double zero = 0.;
  Rocblas::copy_scalar( &zero, A_attr);
  Rocblas::copy_scalar( &zero, bm_attr);
  Rocblas::copy_scalar( &zero, r_attr);

  // Initialize disps to 0
  std::vector< COM::Pane*> allpanes;
  _wrk_window->panes(allpanes);
  for(int i =0, local_npanes = allpanes.size();
      i<local_npanes; ++i){
    Vector_3<double> *ds = reinterpret_cast<Vector_3<double>*>
      ( allpanes[i]->attribute(disps_id)->pointer());
    for(int i =0, ni = allpanes[i]->size_of_real_nodes(); i<ni; ++i){
      ds[i] = Vector_3<double>(0.0,0.0,0.0);
    }
  }
  
  // Loop through the panes and its real faces, calculate A, b, and c = sum_i(wi)
  std::vector< COM::Pane*>::iterator it = allpanes.begin();
  SURF::Window_manifold_2::PM_iterator pm_it=_wm->pm_begin();
  for ( int i=0, local_npanes = allpanes.size(); 
        i<local_npanes; ++i, ++it, ++pm_it) { 
    COM::Pane *pane = *it;

    // Obtain nodal coordinates of current pane, assuming contiguous layout
    const Vector_3<double> *pnts = reinterpret_cast<Vector_3<double>*>
      (pane->attribute(COM_NC)->pointer());
    Vector_3<double> *nrms = reinterpret_cast<Vector_3<double>*>
      ( pane->attribute(facenormals_id)->pointer());

    // Obtain pointers to A_attr, bm_attr, and r_attr 
    Vector_3<double> *As = reinterpret_cast<Vector_3<double>*>
      ( pane->attribute(A_attr->id())->pointer());
    Vector_3<double> *bs_m = reinterpret_cast<Vector_3<double>*>
      ( pane->attribute(bm_attr->id())->pointer());
    Vector_3<double> *aws = reinterpret_cast<Vector_3<double>*>
      ( pane->attribute(aw_attr->id())->pointer());
    Vector_3<double> *uws = reinterpret_cast<Vector_3<double>*>
      ( pane->attribute(uw_attr->id())->pointer());
    double   *rs = reinterpret_cast<double*>
      ( pane->attribute(r_attr->id())->pointer());
    Vector_3<double> *cnts = reinterpret_cast<Vector_3<double>*>
      ( pane->attribute(facecenters_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.
      const Vector_3<double> &ns_nz = nrms[j];
      Vector_3<double> nrm = ns_nz;

      int ne = ene.size_of_edges();
      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;
        
        Vector_3<double> cnt = cnts[j];

        // Calculate the weight of the current face for this node.
        double w=1;

        
        Vector_3<double> ts[2] = { pnts[uindex]-pnts[vindex], 
                                   pnts[windex]-pnts[vindex]};
        
        double s = std::sqrt((ts[0]*ts[0])*(ts[1]*ts[1]));
        double angle=0;
        if ( s>=0) {
          double cosw = ts[0]*ts[1]/s; 
          if (cosw>1) cosw=1; else if ( cosw<-1) cosw=-1;
          angle = std::acos( cosw);
        }
        
        if ( _wght_scheme!=SURF::E2N_ONE) {         
          
          switch ( _wght_scheme) {
          case SURF::E2N_ANGLE: 
          case SURF::E2N_SPHERE: {
            
            if ( _wght_scheme==SURF::E2N_SPHERE)
              w = std::sin(angle)/s; 
            else if ( w>half_pi) 
              w = std::max(0.,pi-angle); 
            else
              w = angle;
            break;
          }
          case SURF::E2N_AREA: {
            w = Vector_3<double>::cross_product(ts[0],ts[1]).norm()/2.; break;
          }
          }
          nrm = w*ns_nz;
        }

        // Obtain the addresses of the linear system and weight
        Vector_3<double> *A_v = &As[3*vindex];
        Vector_3<double> &b_mv = bs_m[vindex];
        Vector_3<double> &aw = aws[vindex];
        Vector_3<double> &uw = uws[vindex];
        double   &r_v  = rs[vindex];

        // Update A, b_m, and r 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;
        b_mv   -= nrm;
        aw     += (std::sin(angle)/s)*ns_nz;
        uw     += ns_nz;
        r_v    += angle;

        // Update indices for next node
        uindex=vindex; vindex=windex;
      }

#if defined(BOUNDARY_TREATMENT) && BOUNDARY_TREATMENT
      // Loop through the halfedges of the face to handle border edges
      Halfedge h( &*pm_it, Edge_ID(j+1,0), SURF::REAL_PANE), h0=h;
      do {
        if ( h.is_border()) {
          // Incorporate normal plane passing through the border edge into 
          // its incident vertices with weight pi/2.
          Vector_3<double> nrm_nz = Vector_3<double>::cross_product( ns_nz, h.tangent()).normalize();
          Vector_3<double> nrm = half_pi*nrm_nz;

          // Update the origin and destination of h with the orthogonal pane
          for ( int k=0; k<2; ++k) {
            int vindex = ene[ (k?h:h.next()).id().lid()]-1;
            // Obtain the addresses of the linear system and weight
            Vector_3<double> *A_v = &As[3*vindex];
            Vector_3<double> &b_mv = bs_m[vindex]; 
            double   &r_v  = rs[vindex];

            // Update A, b_m, and r for the current vertex. 
            // No need to update as sigma is 0 for it.
            A_v[0] += nrm_nz[0]*nrm;
            A_v[1] += nrm_nz[1]*nrm;
            A_v[2] += nrm_nz[2]*nrm;
            b_mv   -= nrm;
            r_v    += half_pi;
          }
        }
      } while ( (h=h.next()) != h0);
#endif
    }
  }

  // Reduce on shared nodes for A, b_m, and r
  _wm->reduce_on_shared_nodes( A_attr, SURF::Window_manifold_2::OP_SUM);
  _wm->reduce_on_shared_nodes( bm_attr, SURF::Window_manifold_2::OP_SUM);
  _wm->reduce_on_shared_nodes( aw_attr, SURF::Window_manifold_2::OP_SUM);
  _wm->reduce_on_shared_nodes( uw_attr, SURF::Window_manifold_2::OP_SUM);
  _wm->reduce_on_shared_nodes( r_attr, SURF::Window_manifold_2::OP_SUM);

  // Loop through the panes and its real vertices
  it = allpanes.begin(); 
  for (int i=0, local_npanes = allpanes.size(); i<local_npanes; ++i, ++it) { 

    COM::Pane *pane = *it;

    // Obtain pointers to A_attr, bo_attr, bm_attr, and r_attr 
    Vector_3<double> *As = reinterpret_cast<Vector_3<double>*>
      ( pane->attribute(A_attr->id())->pointer());
    Vector_3<double> *bs_m = reinterpret_cast<Vector_3<double>*>
      ( pane->attribute(bm_attr->id())->pointer());
    double   *rs = reinterpret_cast<double*>
      ( pane->attribute(r_attr->id())->pointer());
    int      *tranks = reinterpret_cast<int*>
      ( pane->attribute(tangranks_id)->pointer());
    int      *cranks = reinterpret_cast<int*>
      ( pane->attribute(cntnranks_id)->pointer());
    Vector_3<double> *cvecs = reinterpret_cast<Vector_3<double>*>
      ( pane->attribute(cntnvecs_id)->pointer());
    int      *cnstrs = _cnstr_types?reinterpret_cast<int*>
      ( pane->attribute(_cnstr_types->id())->pointer()):NULL;
    Vector_3<double> *cnstr_dirs = _cnstr_dirs?reinterpret_cast<Vector_3<double>*>
      ( pane->attribute(_cnstr_dirs->id())->pointer()):NULL;
    Vector_3<double>  *vnrms = reinterpret_cast<Vector_3<double>*>
      ( pane->attribute(vnormals_id)->pointer());
    Vector_3<double>  *awnrms = reinterpret_cast<Vector_3<double>*>
      ( pane->attribute(awnormals_id)->pointer());
    Vector_3<double>  *uwnrms = reinterpret_cast<Vector_3<double>*>
      ( pane->attribute(uwnormals_id)->pointer());


    // Loop through all real nodes of the pane
    //std::cout << "Number of nodes = " << pane->size_of_real_nodes() << std::endl;
    for ( int j=0, jn=pane->size_of_real_nodes(); j<jn; ++j) {

      // Skip vertices with zero weight, such as those at
      // edge- and face-centers for quadratic elements
      if ( rs[j]==0) continue;
      Vector_3<double> *es = &As[3*j];    // eigen-vectors of current vertex.
      Vector_3<double> *cs = &cvecs[2*j]; // cnstr-tangent of current vertex.
      // Make copy of matrix A as it will be overwritten by eigenvectors
      Vector_3<double> A_v[3] = {es[0], es[1], es[2]};
      // Likewise, b will be overwritten by eigenvalues
      Vector_3<double> b_v = bs_m[j];
      awnrms[j].normalize();
      uwnrms[j].normalize();
      

      // Perform eigen-analysis for the vertex using medial quadric
      // orank is the rank of the orthogonal (primar) space of A.
      // The non-constrained solution is stored in vnrm.
      int orank = eigen_analyze_vertex( es, bs_m[j], &vnrms[j], 2*pi-rs[j]);
      // EDIT
      //orank = (orank == 3) ? 2 : orank;

      //COM_assertion_msg((orank <=3) && (orank >= 1), "Orank is invalid\n");

      // Store the rank of the tangent space.
      tranks[j] = 3-orank;

      // Copy eigenvectors into cvecs.
      for ( int k=tranks[j]-1; k>=0; --k)
        cs[k] = es[orank+k];

      // Compute normal estimate using medial quadric, 
      // and save it into vnrms (overwrite the non-constrained solution)
      Vector_3<double> d(0,0,0);
      // Set the tolerance to avoid perturbation due to small eigenvalues
      double tol = bs_m[j][0]*5.e-3;
      for ( int k=0; k<orank; ++k)
        d += (es[k]*b_v/-std::max(bs_m[j][k],tol))*es[k];
      vnrms[j] = d;
    }
  }
  if(_verb > 1) 
    std::cout << "Exiting Rocmop::compute_medial_quadric" << std::endl;
}

Here is the call graph for this function:

Here is the caller graph for this function:

void compute_medial_quadric ( )

protected

Compute medial quadric for every vertex.

void compute_medial_quadric ( )

protected

Compute medial quadric for every vertex.

void constrain_displacements ( COM::Attribute *  w_disp )

protected

Contrain displacements to _maxdisp.

Definition at line 1424 of file Rocmop_1.C.

References Rocblas::div_scalar(), i, j, MPI_MAX, ni, nj, Vector_3< Type >::norm(), and nvc::norm().

                                                         {
#ifdef ROCPROF
  Profile_begin("Rocmop::const_disp");
#endif
  
  if(_maxdisp > 0.0){
    
    int disp_id = w_disp->id();
    double max_norm = 0.0;

    std::vector<Pane*> allpanes;
    const_cast<COM::Window*>(_usr_window)->panes(allpanes);
    
    for(int i=0,ni = allpanes.size(); i<ni; ++i){
      
      Vector_3<double> *ptr = NULL;
      COM::Attribute* ptr_att = allpanes[i]->attribute(disp_id);
      void* void_ptr = ptr_att->pointer();
      ptr = reinterpret_cast<Vector_3<double>*>(void_ptr);

      for(int j=0,nj = allpanes[i]->size_of_real_nodes(); j<nj; ++j){   
        double norm = ptr[j].norm();
        if(norm > max_norm)
          max_norm = norm;      
      }
    }

    agree_double(max_norm,MPI_MAX);
    
    if(max_norm > _maxdisp){
      double div = max_norm/_maxdisp;
      Rocblas::div_scalar(const_cast<const COM::Attribute*>(w_disp),&div,w_disp);
    }
  }
#ifdef ROCPROF
  Profile_end("Rocmop::const_disp");
#endif
}

Here is the call graph for this function:

void constrain_displacements ( COM::Attribute *  w_disp )

protected

Contrain displacements to _maxdisp.

void determine_pane_border

(

)

Determine which nodes and elements are on pane borders.

void determine_pane_border

(

)

Determine which nodes and elements are on pane borders.

Definition at line 603 of file Rocmop.C.

References i.

                                  {
  if(_verb) 
    std::cout << "Entering Rocmop::determine_pane_border" << std::endl;
  
  std::vector<const COM::Pane*> allpanes;
  _wrk_window->panes(allpanes);
  int local_npanes = (int)allpanes.size();

  _is_pane_bnd_node.resize(local_npanes);
  _is_pane_bnd_elem.resize(local_npanes);

  for(int i=0; i< local_npanes; ++i){
    int size_of_real_nodes = allpanes[i]->size_of_real_nodes();
    int size_of_real_elems = allpanes[i]->size_of_real_elements();
    _is_pane_bnd_node[i].resize(size_of_real_nodes,0);

    std::vector<bool> is_isolated; // is a node isolated?
    MAP::Pane_boundary pb (allpanes[i]);
    pb.determine_border_nodes(_is_pane_bnd_node[i], is_isolated);
  }

  mark_elems_from_nodes(_is_pane_bnd_node,_is_pane_bnd_elem);
  
  if(_verb > 1) 
    std::cout << "Exiting Rocmop::determine_pane_border" << std::endl;    
}

void determine_pane_border

(

)

Determine which nodes and elements are on pane borders.

void determine_physical_border

(

COM::Attribute *

w_is_surface

)

Determine which nodes are on the physical border.

void determine_physical_border

(

COM::Attribute *

w_is_surface

)

Determine which nodes are on the physical border.

Definition at line 473 of file smooth_mesquite_ng.C.

References COM_assertion_msg, COM_compatible_types(), COM_INT, COM_PCONN, else, i, j, k, max(), nj, nk, Rocmap::reduce_maxabs_on_shared_nodes(), total_npanes, and Rocmap::update_ghosts().

                                                                {
  COM_assertion_msg( w_is_surface, "Unexpected NULL pointer");
  COM_assertion_msg( COM_compatible_types( w_is_surface->data_type(), COM_INT),
                     "Surface-list must have integer type");
  COM_assertion_msg( w_is_surface->is_nodal() == 1,
                     "Surface-list must be nodal");
  COM_assertion_msg( w_is_surface->size_of_components() == 1,
                     "Surface-list must have a single component");
  int w_is_surface_id = w_is_surface->id();

  // This function is implemented as follows:
  // 0 Declare and resize data structures.
  // 
  // Loop through panes
  //   1 Get the pane boundary nodes and faces using Rocmap::Pane_boundary
  //   2 Identify adjacent panes, and shared nodes.  Also, create an id mappings
  //     from a node's id to its index in the shared node array and vice versa.
  //   3 Identify "possibly shared faces", those faces whose nodes are all
  //     shared with the same adjacent pane. Make a list of these faces.
  // 4 Communicate the size of the local face list to send to / recieve from
  //   adjacent panes.
  // 5 Create buffers for the incoming face list and communicate faces lists.
  // 6 For every face received, remove any matching faces in the list of 
  //   boundary faces.
  // 7 Mark all nodes still contained in boundary faces as physical boundary
  //   nodes.


  // 0 Declare and resize data structures.

  // vector corresponds to adjacent panes (reused amongst local panes)
  std::vector<bool> is_border; // is a node a border node?
  std::vector<bool> is_shared; // is a node a shared node?
  std::vector<bool> is_isolated; // is a node isolated?
  std::vector<MAP::Facet_ID> border_facet_list;

  // Space in which false pconns are built. (reused amongs local panes)
  // outside vectors correspond to adjacent panes
  // inside vectors correspond to connectivity information
  //std::vector<std::vector<int> > false_pconn;

  // Per-local-pane data structures
  // outside vectors correspond to the local panes
  // inside vectors correspond to adjacent panes
  std::vector<std::set<MAP::Facet_ID> > border_facets; // List of all border faces.
  std::vector<std::vector<Corners2Face_Map> > maybe_shared; // Possibly shared facets
  std::vector<std::vector<Id_Map > > lid2ind; // map from local id to pconn index
  std::vector<std::vector<Id_Map > > ind2lid; // map from pconn index to local id
  std::vector<std::vector<std::set<int> > > adj_pane_set; // nodes shared with adjacent panes
  std::vector<std::vector<int> > adj_pane_id; // ids of adjacent panes
  std::vector<std::vector<int> > adj_pane_recv; // amount of data coming from adjacent panes.

  // The amount of data to be sent by each pane
  std::vector<int> send_sizes;

  std::set<MAP::Facet_ID>::iterator bf_it, bf_it2;
  Id_Map::iterator idm_it,idm_it2; // iterator for lid2ind or ind2lid
  Corners2Face_Map::iterator ms_it,ms_it2;
  std::set<int >::iterator aps_it, aps_it2; //iterator for adj_pane_map

  // get the pconn offset
  int pconn_offset = MAP::Pane_connectivity::pconn_offset();

  // Resize per-local-pane data structures
  std::vector<COM::Pane*> allpanes;
  COM::Window * wrk_window = w_is_surface->window();
  wrk_window->panes(allpanes);
  int total_npanes = (int)allpanes.size(); 
  border_facets.resize(total_npanes);
  maybe_shared.resize(total_npanes);
  lid2ind.resize(total_npanes);
  ind2lid.resize(total_npanes);
  adj_pane_set.resize(total_npanes);
  adj_pane_id.resize(total_npanes);
  adj_pane_recv.resize(total_npanes);
  send_sizes.resize(total_npanes);

  // Register an attribute for our fake pconn.
  COM::Attribute * false_pconn = wrk_window->new_attribute("false_pconn", 'p', COM_INT, 1, "");

  // Register an attribute for sending and receiving sizes
  // also used for sending buffer
  COM::Attribute *com_buff = wrk_window->new_attribute("com_buff", 'p', COM_INT, 1,"");
  int w_com_buff_id = com_buff->id();

  // Loop through panes
  for(int i=0; i<total_npanes; ++i){
    int size_of_real_nodes = allpanes[i]->size_of_real_nodes();
    is_border.clear();
    is_shared.clear();
    is_isolated.clear();
    is_border.resize(size_of_real_nodes,0);
    is_shared.resize(size_of_real_nodes,0);
    is_isolated.resize(size_of_real_nodes,0);
    send_sizes[i] = 0;
    
    //   1 Get the pane boundary nodes and faces using Rocmap::Pane_boundary
    // Determine the border nodes.
    MAP::Pane_boundary pb(allpanes[i]);
    pb.determine_border_nodes(is_border, is_isolated, &border_facet_list);

    // put the list of border facets into a set for faster searching
    for(int j =0, nj=border_facet_list.size(); j<nj; ++j)
      border_facets[i].insert(border_facet_list[j]);
    
    //   2 Identify adjacent panes, and shared nodes.  Also, create an id mappings
    //     from a node's id to its index in the shared node array and vice versa.

    // Obtain the pane connectivity of the local pane
    const COM::Attribute *pconn = allpanes[i]->attribute(COM::COM_PCONN);
    const int *vs = (const int*)pconn->pointer() + pconn_offset;
    int vs_size = pconn->size_of_real_items() - pconn_offset;

    // Determine the number of communicating panes for shared nodes.
    int count=0;
    for (int j=0, nj=vs_size; j<nj; j+=vs[j+1]+2) {
      if (wrk_window->owner_rank( vs[j]) >=0) ++count;
    }

    // Resize per-adjacent-pane data structures
    //false_pconn.resize(count);
    maybe_shared[i].resize(count);
    lid2ind[i].resize(count);
    ind2lid[i].resize(count);
    adj_pane_set[i].resize(count);
    adj_pane_id[i].resize(count);
    adj_pane_recv[i].resize(count);
    
    int index = 0;
    // Loop through communicating panes for shared nodes.
    for ( int j=0; j<count; ++j, index+=vs[index+1]+2) {
      // We skip the panes that are not in the current window 
      while ( wrk_window->owner_rank(vs[index])<0) {
        index+=vs[index+1]+2;
        COM_assertion_msg( index<=vs_size, "Invalid communication map");
      } 
      adj_pane_id[i][j] = vs[index];
      // Add this pane to current shared node's list, mark node as shared and
      // keep track of the mapping between array index and local node id
      for(int k=0; k<vs[index+1]; ++k){
        is_shared[vs[index+2+k]-1]=1;
        adj_pane_set[i][j].insert(vs[index+2+k]);
        lid2ind[i][j].insert(std::make_pair(vs[index+2+k],k));
        ind2lid[i][j].insert(std::make_pair(k,vs[index+2+k]));
      }
    } 
#if 0

    std::cout << "Pane " << allpanes[i]->id() << " lid2ind = \n";
    for(int j=0; j < count; ++j){
      std::cout << " to Pane " << adj_pane_id[i][j] << " =\n";
      for(int k =0; k < ind2lid[i][j].size();++k){
        std::cout << "(" << k << "," << ind2lid[i][j][k] << ") ";
      }
      std::cout << "\n";
    }
    std::cout << "Pane " << i << "\n"
              << " shared nodes = \n";

    for(int j=0; j<is_shared.size(); ++j){
      if(is_shared[j])
        std::cout << j+1 << " ";
    }
    std::cout << "\n\n nodes shared with othe panes";

    aps_it = adj_pane_set[i][0].begin();
    for(; aps_it != adj_pane_set[i][0].end(); ++aps_it){
      std::cout << *aps_it << " ";
    }

#endif

    //   3 Identify "possibly shared faces", those faces whose nodes are all
    //     shared with the same adjacent pane. Make a list of these faces.

    bf_it2 = border_facets[i].end();
    bf_it = border_facets[i].begin();
    //for(int j=0, nj=border_facets[i].size(); j<nj; ++j){
    for(; bf_it != bf_it2; ++bf_it){
      Element_node_enumerator ene(allpanes[i], (*bf_it).eid());
      Facet_node_enumerator fne (&ene, (*bf_it).lid());
      //if all the nodes are shared
      if( is_shared[fne[0]-1] &&
          is_shared[fne[1]-1] &&
          is_shared[fne[2]-1] &&
          fne.size_of_edges()>3?is_shared[fne[3]-1]:1
          ){
        // then see if they are all shared by the same panes
        Four_tuple ns( fne[0], fne[1], fne[2], fne.size_of_edges()>3?fne[3]:-1);
        for(int k=0; k<count; ++k){
          aps_it2 = adj_pane_set[i][k].end();
          if(aps_it2 != adj_pane_set[i][k].find(fne[0]) &&
             aps_it2 != adj_pane_set[i][k].find(fne[1]) &&
             aps_it2 != adj_pane_set[i][k].find(fne[2]) &&
             (fne.size_of_edges()>3?(adj_pane_set[i][k].find(fne[3]) != aps_it2):1)){
            // and if they are
            // then this facet is possibly shared with the adjacent node
            Four_tuple ns( fne[0], fne[1], fne[2], fne.size_of_edges()>3?fne[3]:-1);
            std::sort(&ns[0],&ns[4]);
            maybe_shared[i][k].insert(std::make_pair(ns,(*bf_it)));
          }
        }
      }
    }
    
#if 0
    // print out the list of possibly shared nodes:
    std::cout << "Pane " << allpanes[i]->id() << " possibly shared panes\n";
    for(int j = 0; j<count; ++j){
      std::cout << "  faces possibly shared with Pane " << adj_pane_id[i][j] << "\n";
      ms_it2 = maybe_shared[i][j].end();
      for(ms_it = maybe_shared[i][j].begin(); ms_it != ms_it2; ++ms_it){
        std::cout << "(" << (ms_it->first)[0] << " " 
                  << (ms_it->first)[1] << " "
                  << (ms_it->first)[2] << " "
                  << (ms_it->first)[3] << ")  ";          
      }
      std::cout << "\n";
    }
#endif

    wrk_window->set_size("com_buff", allpanes[i]->id(), count, 0);
    int* my_buff;
    wrk_window->resize_array("com_buff", (const int)allpanes[i]->id(),
                             reinterpret_cast<void**>(&my_buff));

    // find the size of the pconn to create
    // 3*(pconn_offset)           // number of communicating blocks
    // + 3                        // shared node block
    // + 6*count                  // (pane id,node count,array-id for each adj. pane)*2blocks

    int my_size = 3*pconn_offset + 3 + 6*count;
    int* my_pconn;
    wrk_window->set_size("false_pconn", allpanes[i]->id(), my_size, 
                         my_size-(3+pconn_offset));
    wrk_window->resize_array("false_pconn", (const int)allpanes[i]->id(),
                             reinterpret_cast<void**>(&my_pconn));
    // Add shared node data
    if(pconn_offset){my_pconn[0]=1;++my_pconn;}
    my_pconn[0]=1; my_pconn[1]=1; my_pconn[2]=1;
    my_pconn += 3;

    // filling in send buffer and real nodes to send info.
    if(pconn_offset){my_pconn[0]=count;++my_pconn;}
    for(int j=0; j<count; ++j){
      // print communicating pane id
      my_pconn[0] = adj_pane_id[i][j];
      // print number of ints to communicate
      my_pconn[1] = 1;
      // print index-id of information
      my_pconn[2] = j+1;
      my_pconn+=3;
      // store size of information to send
      my_buff[j] = 4*maybe_shared[i][j].size();
      send_sizes[i] += my_buff[j];
    }

    // add ghost nodes to receive info.
    if(pconn_offset){my_pconn[0]=count;++my_pconn;}
    for(int j=0; j<count; ++j){
      my_pconn[0] = adj_pane_id[i][j];
      my_pconn[1] = 1;
      my_pconn[2] = j+1;
      my_pconn+=3;
    }
    
#if 0
    for ( int j=0; j<count; ++j, index+=vs[index+1]+2) {
      // We skip the panes that are not in the current window 
      while ( wrk_window->owner_rank(vs[index])<0)
        index+=vs[index+1]+2;
      // List the communicating pane, 1 item to send, and the index of that item in a 
      // local panel array
      my_pconn[3*j] = vs[index];
      my_pconn[3*j+1] = 1;
      my_pconn[3*j+2] = j+1;

      my_pconn[3*(count+j)+pconn_offset] = vs[index];
      my_pconn[3*(count+j)+1+pconn_offset] = 1;
      my_pconn[3*(count+j)+2+pconn_offset] = j+1;

      my_pconn[3*(2*count+j)+2*pconn_offset] = vs[index];
      my_pconn[3*(2*count+j)+1+2*pconn_offset] = 1;
      my_pconn[3*(2*count+j)+2+2*pconn_offset] = j+1;


      // store the size of the buffer to send across
      my_buff[j] = 4*maybe_shared[i][j].size();
      send_sizes[i] += my_buff[j];
    }
#endif
#if 0
    std::cout << "pconn_offset = " << pconn_offset << "\n";
    std::cout << "total size = " << my_size << "\n";
    std::cout << "ghost sized = " << my_size - pconn_offset << "\n\n";
#endif
  }
  wrk_window->init_done();

#if 0
  //print out pconn, for debugging.
  for(int i=0; i<total_npanes; ++i){
    std::cout << "FALSE pconn for pane " << allpanes[i]->id() << std::endl;
    COM::Attribute *my_pconn = allpanes[i]->attribute(w_false_pconn_id);
    int* ptr = (int*)my_pconn->pointer();
    for(int k=0,nk=my_pconn->size_of_items();k<nk;++k){
      std::cout << ptr[k] << " ";
    }
    std::cout << "\n\n";

    // print out size of buffer being sent across
    std::cout << "size of face list\n";
    COM::Attribute *my_bsize = allpanes[i]->attribute(w_com_buff_id);
    ptr = (int*)my_bsize->pointer();
    for(int j =0; j<my_bsize->size_of_items(); ++j)
      std::cout << ptr[j] << " ";
    std::cout << "\n\n";
  }
#endif

  // 4 Communicate the size of the local face list to send to / recieve from
  //   adjacent panes.Pane_connectivity::
  Rocmap::update_ghosts(com_buff,
                        false_pconn);

  //print out pconn, for debugging.
#if 0
  for(int i=0; i<total_npanes; ++i){
    COM::Attribute *my_pconn = allpanes[i]->attribute(w_false_pconn_id);
    int* ptr = (int*)my_pconn->pointer();
    for(int k=0,nk=my_pconn->size_of_items();k<nk;++k){
      std::cout << ptr[k] << " ";
    }
    std::cout << "\n\n";
    // print out possibly shared faces
    std::cout << "possibly shared faces\n";
    for(ms_it = maybe_shared[i][0].begin(), ms_it2 = maybe_shared[i][0].end(); 
        ms_it != ms_it2; ++ms_it){
      std::cout << "(" 
                << (ms_it->first)[0] << " "
                << (ms_it->first)[1] << " "
                << (ms_it->first)[2] << " "
                << (ms_it->first)[3] << ") ";
      }
    std::cout << "\n\n";
    // print out size of buffer being sent across
    std::cout << "size of communicated face list\n";
    COM::Attribute *my_bsize = allpanes[i]->attribute(w_com_buff_id);
    ptr = (int*)my_bsize->pointer();
    for(int j =0; j <adj_pane_id[i].size(); ++j)
      std::cout << ptr[j] << " ";
    std::cout << "\n\n";
  }
#endif

  // 5 Create buffers for the incoming face list and communicate faces lists.
  // loop through panes
  for(int i=0; i<total_npanes; ++i){

    // get pane level attributes and pointers
    COM::Attribute *p_com_buff = allpanes[i]->attribute(w_com_buff_id);
    int *com_buff_ptr = (int*)p_com_buff->pointer();

    // find the number of adjacent panes
    int count = adj_pane_id[i].size();

    // find the size of data to be received and sent.
    int recv_size = 0;
    // add up the sizes of all the buffers being sent
    for(int j=0; j< count; ++j){
      recv_size += com_buff_ptr[j];
      adj_pane_recv[i][j] = com_buff_ptr[j];
    }
    int send_size = send_sizes[i];

    // find the size of the pconn to create
    // 3*(pconn_offset)           // number of communicating blocks
    // + 3                        // shared node block
    // + 4*count                  // (pane id and node count for each adj. pane)*2blocks
    // + recv_size+send_size      // number of nodes in facet lists
    int pconn_size = 3*(pconn_offset)+3+4*count+recv_size+send_size;

#if 0

    std::cout << " SIZES\nrecv_size = " << recv_size << " send_size = " << send_size
              << " pconn size = " << pconn_size << "\n"
              << " 3*pconn_offset = " << 3*pconn_offset
              << " 4*count = " << 4*count << "\n\n";

#endif

    wrk_window->set_size("com_buff", allpanes[i]->id(),std::max(send_size,recv_size),0);
    wrk_window->resize_array("com_buff", (const int)allpanes[i]->id(),
                             reinterpret_cast<void**>(&com_buff_ptr));    
    int* pconn_ptr;
    wrk_window->set_size("false_pconn", allpanes[i]->id(), pconn_size, 
                         pconn_size - (3+pconn_offset));
    wrk_window->resize_array("false_pconn", (const int)allpanes[i]->id(),
                             reinterpret_cast<void**>(&pconn_ptr));

    // add shared node data
    if(pconn_offset){pconn_ptr[0]=1;++pconn_ptr;}
    pconn_ptr[0]=1; pconn_ptr[1]=1; pconn_ptr[2]=1;
    pconn_ptr += 3;
    
    // now fill in the send buffer and real nodes to send info. in the pconn
    if(pconn_offset){pconn_ptr[0]=count;++pconn_ptr;}
    int com_ind = 0; // index of next place in com_buff
    for(int j =0; j<count; ++j){
      pconn_ptr[0]=adj_pane_id[i][j];
      pconn_ptr[1]=4*maybe_shared[i][j].size();
      pconn_ptr+=2;
      ms_it2 = maybe_shared[i][j].end();
      for(ms_it = maybe_shared[i][j].begin(); ms_it != ms_it2; 
          ++ms_it, com_ind+=4, pconn_ptr+=4){
        // Convert nodal id's to pconn based index id's and place in the comm buffer
        const Four_tuple* ft = &(ms_it->first); 
        if ((*ft)[0]==-1)
          com_buff_ptr[com_ind] = -1;
        else
          idm_it = lid2ind[i][j].find((*ft)[0]);
        idm_it = lid2ind[i][j].find((*ft)[1]);
        com_buff_ptr[com_ind+1] = idm_it->second;
        idm_it = lid2ind[i][j].find((*ft)[2]);
        com_buff_ptr[com_ind+2] = idm_it->second;
        idm_it = lid2ind[i][j].find((*ft)[3]);
        com_buff_ptr[com_ind+3] = idm_it->second;
        std::sort(&com_buff_ptr[com_ind],&com_buff_ptr[com_ind+4]);
        // Also fill in the pconn
        pconn_ptr[0] = com_ind+1;
        pconn_ptr[1] = com_ind+2;
        pconn_ptr[2] = com_ind+3;
        pconn_ptr[3] = com_ind+4;
      }
    }

    // now fill in the ghost nodes to receive info. in the pconn
    if(pconn_offset){pconn_ptr[0]=count;++pconn_ptr;}
    com_ind = 1;
    for(int j=0; j<count; ++j){
      pconn_ptr[0]=adj_pane_id[i][j];
      pconn_ptr[1]=adj_pane_recv[i][j];
      for(int k=0; k<adj_pane_recv[i][j]; ++k){
        pconn_ptr[k+2] = k+com_ind;     
      }   
      com_ind += adj_pane_recv[i][j];
      pconn_ptr += 2+adj_pane_recv[i][j];
    }
  }

#if 0
  //print out pconn, for debugging.
  for(int i=0; i<total_npanes; ++i){
    std::cout << "Facet list FALSE pconn for pane " << allpanes[i]->id() << std::endl;
    COM::Attribute *my_pconn = allpanes[i]->attribute(w_false_pconn_id);
    int* ptr = (int*)my_pconn->pointer();
    for(int k=0,nk=my_pconn->size_of_items();k<nk;++k){
      std::cout << ptr[k] << " ";
    }
    std::cout << "\n\n";

    std::cout << "Facet list buffer for pane " << allpanes[i]->id() << std::endl;
    COM::Attribute *my_buff = allpanes[i]->attribute(w_com_buff_id);
    ptr = (int*)my_buff->pointer();
    for(int k=0,nk=my_buff->size_of_items();k<nk;++k){
      std::cout << ptr[k] << " ";
    }
    std::cout << "\n\n";
  }
#endif

  // update
  Rocmap::update_ghosts(com_buff,
                        false_pconn);

#if 0
  //print out pconn, for debugging.
  for(int i=0; i<total_npanes; ++i){
    std::cout << "Updated Facet list buffer for pane " << allpanes[i]->id() << std::endl;
    COM::Attribute *my_buff = allpanes[i]->attribute(w_com_buff_id);
    int* ptr = (int*)my_buff->pointer();
    for(int k=0,nk=my_buff->size_of_items();k<nk;++k){
      std::cout << ptr[k] << " ";
    }
    std::cout << "\n\n";
  }
#endif
  // 6 For every face received, remove any matching faces in the list of 
  //   boundary faces.
  // loop through all panes
  for(int i =0; i<total_npanes; ++i){
    COM::Attribute *my_buff = allpanes[i]->attribute(w_com_buff_id);
    int* buff_ptr = (int*)my_buff->pointer();

    int count = adj_pane_recv[i].size();
    // loop through adjacent panes
    //std::cout << "Pane " << allpanes[i]->id() << "\n";
    for(int j=0;j<count; buff_ptr += adj_pane_recv[i][j], ++j){

      //std::cout << "Converting node array\n";
      // convert all node array indices to local id
      //std::cout << adj_pane_recv[i][j] << " items being converted\n";
      for(int k=0, nk =adj_pane_recv[i][j]; k<nk; ++k){
        if(buff_ptr[k]!=-1){
          buff_ptr[k] = ind2lid[i][j][buff_ptr[k]];
        }
      }
      //      std::cout << "Checking for matching facets with pane" << adj_pane_id[i][j] << "\n";
      // check for matching facets
      ms_it2 = maybe_shared[i][j].end();
      bf_it2 = border_facets[i].end();
      for(int k=0, nk=adj_pane_recv[i][j]; k<nk; k+=4){
        Four_tuple ns(buff_ptr[k],buff_ptr[k+1],buff_ptr[k+2],buff_ptr[k+3]);
        std::sort(&ns[0],&ns[4]);
        //std::cout << " (" << ns[0] << " " << ns[1] << " " << ns[2] << " " << ns[3] << ") ";
        ms_it = maybe_shared[i][j].find(ns);
        if(ms_it != ms_it2){
          //std::cout << "found ";
          Element_node_enumerator ene(allpanes[i], (ms_it->second).eid());
          Facet_node_enumerator fne (&ene, (ms_it->second).lid());
          //std::cout << "(" << fne[0] << " " << fne[1]
          //        << " " << fne[2];
          //      if(fne.size_of_edges()>3)
          //  std::cout << " " << fne[3];
          //std::cout << ")\n";
          MAP::Facet_ID fid = ms_it->second;
          bf_it = border_facets[i].find(fid);
          if(bf_it != bf_it2)
            border_facets[i].erase(fid);
        }
        else
          ;
          //std::cout << "not found\n";
      }
    }
  }

  //  std::cout << "On to step 7 \n";

  // 7 Mark all nodes still contained in boundary faces as physical boundary
  //   nodes.  
  for(int i =0; i < total_npanes; ++i){
    //    std::cout << "Pane " << i << " detected surface faces\n";
    bf_it2 = border_facets[i].end();
    for(bf_it = border_facets[i].begin(); bf_it!= bf_it2; ++bf_it){
      Element_node_enumerator ene(allpanes[i], (*bf_it).eid());
      Facet_node_enumerator fne (&ene, (*bf_it).lid());
      
      //std::cout << "(" << fne[0] << " " << fne[1]
      //                << " " << fne[2];
      //  if(fne.size_of_edges()>3)
      //        std::cout << " " << fne[3];
      //std::cout << ")";
    }
    //std::cout << "\n\n";
  }

  for(int i = 0; i < total_npanes; ++i){
    COM::Attribute *p_is_surface = allpanes[i]->attribute(w_is_surface_id);
    int* surf_ptr = (int*)p_is_surface->pointer();
    // initialize surface to 0s
    for(int j=0, nj= p_is_surface->size_of_items();j<nj; ++j){
      surf_ptr[j] = 0;
    }
    bf_it2 = border_facets[i].end();
    //std::cout << "Pane " << allpanes[i]->id() << " marking surface nodes\n";
    for(bf_it = border_facets[i].begin(); bf_it!= bf_it2; ++bf_it){
      Element_node_enumerator ene(allpanes[i], (*bf_it).eid());
      Facet_node_enumerator fne (&ene, (*bf_it).lid());
      surf_ptr[fne[0]-1] = 1;
      surf_ptr[fne[1]-1] = 1;
      surf_ptr[fne[2]-1] = 1;
      //std::cout << fne[0] << " " << fne[1] << " " << fne[2] << "\n";
      if(fne.size_of_edges()>3)
        surf_ptr[fne[3]-1] = 1;
    }
    //    std::cout << "\n\n";
  }

#if 0
  for(int i = 0; i < total_npanes; ++i){
    COM::Attribute *p_pconn = allpanes[i]->attribute(COM::COM_PCONN);
    //    std::cout << "PANE " << allpanes[i]->id() << "'s real pconn size = "
    //        << p_pconn->size_of_real_items() << " and ghost size = "
    //        << p_pconn->size_of_items()-p_pconn->size_of_real_items() << "\n\n";

    //int* ptr = (int*)p_pconn->pointer();
    //for(int j=0, nj=p_pconn->size_of_items(); j<nj; ++j){
      //      std::cout << ptr[j] << " ";
    //}
    //    std::cout << "\n\n";
  }
#endif

  // update
  Rocmap::reduce_maxabs_on_shared_nodes(w_is_surface);
  Rocmap::update_ghosts(w_is_surface);

  // delete buffer and false pconn attributes
  wrk_window->delete_attribute(com_buff->name());
  wrk_window->delete_attribute(false_pconn->name());  

}

Here is the call graph for this function:

void determine_physical_border

(

)

Determine which nodes and elements are on the physical border.

void determine_physical_border

(

)

Determine which nodes and elements are on the physical border.

Definition at line 677 of file Rocmop.C.

References COM_assertion_msg, i, and j.

                                      {
  if(_verb) 
    std::cout << "Entering Rocmop::determine_physical_border()" << std::endl;

  const std::string surf_attr("is_surface");
  COM::Attribute* w_is_surface = _wrk_window->attribute(surf_attr);
  COM_assertion_msg( w_is_surface, "Unexpected NULL pointer");
  int is_surface_id = w_is_surface->id();

  std::vector<const COM::Pane*> allpanes;
  _wrk_window->panes(allpanes);
  int local_npanes = (int)allpanes.size();
    
  _is_pane_bnd_node.resize(local_npanes);

  for(int i=0; i < local_npanes; ++i){
    _is_pane_bnd_node[i].resize(allpanes[i]->size_of_real_nodes());

    // get pane level pointer to physical border property.
    const COM::Attribute *p_is_surface = allpanes[i]->attribute(is_surface_id);
    int *is_surface_ptr = (int*)p_is_surface->pointer();
    
    // loop through real nodes
    for(int j=0; j< allpanes[i]->size_of_real_nodes(); ++j){
      if (is_surface_ptr[j])
        _is_pane_bnd_node[i][j] = true;
    }
  }

  mark_elems_from_nodes(_is_pane_bnd_node,_is_pane_bnd_elem);

  if(_verb > 1)
    std::cout << "Exiting Rocmop::determine_physical_border()" << std::endl;
}

void determine_physical_border

(

COM::Attribute *

w_is_surface

)

Determine which nodes are on the physical border.

void determine_physical_border

(

)

Determine which nodes and elements are on the physical border.

void determine_shared_border

(

)

Determine which nodes and elements are shared.

void determine_shared_border

(

)

Determine which nodes and elements are shared.

Definition at line 631 of file Rocmop.C.

References COM_assertion_msg, COM_PCONN, i, j, k, nj, and Pane_connectivity::pconn_offset().

                                    {
  if(_verb) 
    std::cout << "Entering Rocmop::determine_shared_nodes" << std::endl;

  std::vector<const COM::Pane*> allpanes;
  _wrk_window->panes(allpanes);
  int local_npanes = (int)allpanes.size();
    
  _is_shared_node.resize(local_npanes);
  
  //First, get the list of shared nodes.
  for(int i=0; i < (int)(local_npanes); ++i){
    // Obtain the pane connectivity of the local pane.
    const COM::Attribute *pconn = allpanes[i]->attribute(COM::COM_PCONN);
    // Use the pconn offset
    const int *vs = (const int*)pconn->pointer()+Pane_connectivity::pconn_offset();
    int vs_size=pconn->size_of_real_items()-Pane_connectivity::pconn_offset();    
    _is_shared_node[i].resize(allpanes[i]->size_of_real_nodes(),0);
    
    // Determine the number of communicating panes for shared nodes.
    int count=0;
    for (int j=0, nj=vs_size; j<nj; j+=vs[j+1]+2) {
      if (_wrk_window->owner_rank( vs[j]) >=0) ++count;
    }
    
    int index = 0;
    // Loop through communicating panes for shared nodes.
    for ( int j=0; j<count; ++j, index+=vs[index+1]+2) {
      // We skip the panes that are not in the current window 
      while ( _wrk_window->owner_rank(vs[index])<0) {
        index+=vs[index+1]+2;
        COM_assertion_msg( index<=vs_size, "Invalid communication map");
      } 
      // Add shared nodes to the list
      for(int k=0; k<vs[index+1]; ++k){
        _is_shared_node[i][vs[index+2+k]-1] = 1;
      }
    }
  }

  mark_elems_from_nodes(_is_shared_node,_is_shared_elem);

  if(_verb > 1) 
    std::cout << "Exiting Rocmop::determine_shared_nodes" << std::endl;
}

Here is the call graph for this function:

void determine_shared_border

(

)

Determine which nodes and elements are shared.

int eigen_analyze_vertex ( Vector_3< double >  A_io[3], Vector_3< double > &  b_io, Vector_3< double > *  nrm_nz, double  angle_defect  )

protected

Definition at line 1081 of file smooth_medial.C.

References _dir_thres, _eig_thres, _reorthog, _saliency_crn, NTS::abs(), cimg_library::atan(), compute_eigenvectors(), k, max(), Vector_3< Type >::normalize(), pi, sqrt(), and x.

Referenced by compute_medial_quadric().

                                                           {

  // A_io will be replaced by its eigenvectors.
  Vector_3<double> *es = A_io;

  // Create a backup of b_io, as b_io will be replaced by eigenvalues.
  const Vector_3<double> b = b_io;
  Vector_3<double> &lambdas = b_io;

  // Compute the eigenvalues and eigenvectors of the matrix. 
  // Note that lambdas will be in decreasing order!
  compute_eigenvectors( es, lambdas);

  int nrank;
  double eps = lambdas[0]*1.e-8;
  double gs[] = { b*es[0]/lambdas[0],
                  b*es[1]/std::max(eps,lambdas[1]), 
                  b*es[2]/std::max(eps,lambdas[2]) };
  double gs_abs[] = {std::abs(gs[0]), std::abs(gs[1]), std::abs(gs[2])};
  
  // Classify offset intersection based on acos(-b.norm()/rv) and lambdas
  if ( lambdas[2] >= _saliency_crn *
       std::max(lambdas[0]-lambdas[1], lambdas[1]-lambdas[2]) || ad >= .5*pi)
    nrank = 3;
  else if ( gs_abs[1] >= gs_abs[0] || lambdas[0]*_dir_thres<lambdas[1]) {
    if ( lambdas[1] < lambdas[0]*_eig_thres) {
      lambdas[1] = lambdas[0]*_eig_thres;
      std::cerr << "Rocprop Warning: Mesh contains near cusp." 
                << "lambdas[1]=" << lambdas[1] << " and lambdas[0]="
                << lambdas[0] << std::endl;
    }
    nrank = 2; // ridge vertex
#if PRINT_FEATURE
    if ( gs_abs[1] < gs_abs[0]) {
      double theta = (360/pi)*std::atan(std::sqrt(lambdas[1]/lambdas[0]));
      if ( theta<fangle_min_eigen) fangle_min_eigen=theta;
      if ( theta>fangle_max_eigen) fangle_max_eigen=theta;
    }
#endif
  }
  else
    nrank = 1; // smooth vertex

  if ( !_reorthog || nrank==1)
    *nrm_nz=( es[0]*b>0)?-es[0]:es[0];
  else {
    // Solve for x within primary space
    Vector_3<double> x(0,0,0);
    for ( int k=0; k<nrank; ++k) x -= gs[k]*es[k];
    *nrm_nz = x.normalize();
  }
  return nrank;
}

Here is the call graph for this function:

Here is the caller graph for this function:

int eigen_analyze_vertex ( Vector_3< double >  A_io[3], Vector_3< double > &  b_io, Vector_3< double > *  nrm_nz, double  angle_defect  )

protected

int eigen_analyze_vertex ( Vector_3< double >  A_io[3], Vector_3< double > &  b_io, Vector_3< double > *  nrm_nz, double  angle_defect  )

protected

void evaluate_face_normals ( )

protected

Evaluate face normals (copied from FaceOffset_3.[hC].

Definition at line 903 of file smooth_medial.C.

References _verb, _wrk_window, COM_assertion_msg, COM_NC, Vector_3< Type >::cross_product(), i, j, Element_node_enumerator::next(), nj, Vector_3< Type >::normalize(), Element_node_vectors_k_const< Value >::set(), Element_node_enumerator::size_of_edges(), and Element_node_enumerator::size_of_nodes().

Referenced by smooth_surf_medial().

                                   {

  if(_verb > 1) 
    std::cout << "Entering Rocmop::evaluate_face_normals" << std::endl;

  COM_assertion_msg(_wrk_window, "Unexpected NULL pointer encountered.");

  if(_verb > 3) 
    std::cout << "Getting attribute ids" << std::endl;

  const std::string att1("facenormals");
  const std::string att2("facecenters");
  int facenormals_id = _wrk_window->attribute(att1)->id();
  int facecenters_id = _wrk_window->attribute(att2)->id();

  // Loop through the panes and its real faces
  std::vector< COM::Pane*> allpanes;
  _wrk_window->panes(allpanes);

  std::vector< COM::Pane*>::iterator it = allpanes.begin();
  for ( int i=0, local_npanes = allpanes.size(); i<local_npanes; ++i, ++it) { 
    COM::Pane *pane = *it;

    // Obtain nodal coordinates of current pane, assuming contiguous layout
    COM_assertion_msg( pane->attribute( COM_NC)->stride()==3 ||
                       pane->size_of_real_nodes()==0, 
                       "Coordinates must be stored in contiguous layout.");

    const Vector_3<double> *pnts = reinterpret_cast< Vector_3<double>* >
      (pane->attribute(COM_NC)->pointer());
    Vector_3<double> *nrms = reinterpret_cast< Vector_3<double>* >
      ( pane->attribute(facenormals_id)->pointer());
    Vector_3<double> *cnts = reinterpret_cast< Vector_3<double>* >
      ( pane->attribute(facecenters_id)->pointer());

    COM_assertion_msg(pnts, "NULL pointer to nodal coords.");
    COM_assertion_msg(nrms, "NULL pointer to face normals.");
    COM_assertion_msg(cnts, "NULL pointer to face centers.");    

    // 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()) {

      // Obtain current face normal
      Element_node_vectors_k_const<Vector_3<double> > ps; ps.set( pnts, ene, 1);
      SURF::Generic_element_2 e(ene.size_of_edges(), ene.size_of_nodes());

      Vector_2<double> nc(0.5,0.5);
      Vector_3<double> J[2];

      e.Jacobian( ps, nc, J);

      nrms[j] = Vector_3<double>::cross_product( J[0], J[1]);

      nrms[j].normalize();

      e.interpolate_to_center( ps, &cnts[j]);

    }
  }

  if(_verb > 1) 
    std::cout << "Exiting Rocmop::evaluate_face_normals" << std::endl;
}

Here is the call graph for this function:

Here is the caller graph for this function:

void evaluate_face_normals ( )

protected

Evaluate face normals (copied from FaceOffset_3.[hC].

void evaluate_face_normals ( )

protected

Evaluate face normals (copied from FaceOffset_3.[hC].

void get_constraint_directions ( int  type, const Vector_3< double > &  dir, int &  ndirs, Vector_3< double >  dirs[2]  )

staticprotected

Get orthonormals of the constraints.

If nnrms is 0, then not constrained. If nnrms is 3, then the point is fixed. Otherwise, the point is constrained either in a plane or a line.

Definition at line 1245 of file smooth_medial.C.

References NTS::abs(), COM_assertion_msg, s, sqrt(), and Vector_3< Type >::squared_norm().

                                                                 {
  if ( type == 1 || type == -1)
    COM_assertion_msg( &dir && std::abs(dir.squared_norm()-1)<=1.e-6, 
                       "Constrained direction must be normalized");

  switch ( type) {
  case 0: 
    ndirs = 0;
    break;
  case 2: 
    ndirs = 3;
    break;
  case 1:
    ndirs = 1; dirs[0] = dir;
    break;
  case -1: {
    ndirs = 2;

    double sqnrm;
    if ( std::abs( dir[2])>0.7) {
      sqnrm = dir[1]*dir[1]+dir[2]*dir[2];
      dirs[0] = Vector_3<double>( 0, -dir[2], dir[1]);
      dirs[1] = Vector_3<double>( sqnrm, -dir[0]*dir[1], -dir[0]*dir[2]);
    }
    else {
      sqnrm = dir[0]*dir[0]+dir[1]*dir[1];
      dirs[0] = Vector_3<double>( dir[1], -dir[0], 0);
      dirs[1] = Vector_3<double>( -dir[0]*dir[2], -dir[1]*dir[2], sqnrm);
    }

    double s = 1./std::sqrt( sqnrm);
    dirs[0] *= s; dirs[1] *= s;
    break;
  }
  case 'x':  
    ndirs = 1; 
    dirs[0] = Vector_3<double>(1,0,0); 
    break;
  case -'x': 
    ndirs = 2; 
    dirs[0] = Vector_3<double>(0,1,0); 
    dirs[1] = Vector_3<double>(0,0,1); 
    break;
  case 'y': 
    ndirs = 1; 
    dirs[0] = Vector_3<double>(0,1,0); 
    break;
  case -'y':
    ndirs = 2;
    dirs[0] = Vector_3<double>(1,0,0); 
    dirs[1] = Vector_3<double>(0,0,1); 
    break;
  case 'z': 
    ndirs = 1; 
    dirs[0] = Vector_3<double>(0,0,1); 
    break;
  case -'z':
    ndirs = 2;
    dirs[0] = Vector_3<double>(1,0,0); 
    dirs[1] = Vector_3<double>(0,1,0); 
    break;
  default: COM_assertion_msg( false, "Unknown type of constraint");
  }
}

Here is the call graph for this function:

static void get_constraint_directions ( int  type, const Vector_3< double > &  dir, int &  ndirs, Vector_3< double >  dirs[2]  )

staticprotected

Get orthonormals of the constraints.

If nnrms is 0, then not constrained. If nnrms is 3, then the point is fixed. Otherwise, the point is constrained either in a plane or a line.

static void get_constraint_directions ( int  type, const Vector_3< double > &  dir, int &  ndirs, Vector_3< double >  dirs[2]  )

staticprotected

Get orthonormals of the constraints.

If nnrms is 0, then not constrained. If nnrms is 3, then the point is fixed. Otherwise, the point is constrained either in a plane or a line.

void get_redist_safe_factor ( COM::Attribute *  c_attr, COM::Attribute *  w_attr, int  rank  )

protected

Definition at line 825 of file smooth_medial.C.

References _wrk_window, COM_NC, Rocblas::copy_scalar(), Vector_3< Type >::cross_product(), i, j, k, max(), Element_node_enumerator::next(), nj, Vector_3< Type >::normalize(), Window_manifold_2::OP_MAX, Element_node_enumerator::size_of_edges(), and solve().

Referenced by redistribute_vertices_ridge(), and redistribute_vertices_smooth().

                                  {
  double zero=0;
  Rocblas::copy_scalar( &zero, w_attr);

  const std::string att1("tangranks"), att2("weights2"),
    att3("barycrds"), att4("PNelemids");
  int tangranks_id = _wrk_window->attribute(att1)->id(),
    weights2_id = _wrk_window->attribute(att2)->id(),
    barycrds_id = _wrk_window->attribute(att3)->id(),
    PNelemid_id = _wrk_window->attribute(att4)->id();

  std::vector< COM::Pane*> allpanes;
  _wrk_window->panes(allpanes);
  std::vector< COM::Pane*>::iterator it = allpanes.begin();
  for (int i=0, local_npanes = allpanes.size(); i<local_npanes; ++i, ++it) { 
    COM::Pane *pane = *it;
    
    const Point_3<double> *pnts = reinterpret_cast<Point_3<double>*>
      (pane->attribute(COM_NC)->pointer());
    Vector_3<double> *cs = reinterpret_cast<Vector_3<double>*>
      ( pane->attribute(c_attr->id())->pointer());
    double   *ws = reinterpret_cast<double*>
      ( pane->attribute(w_attr->id())->pointer());
    int *tranks = reinterpret_cast<int*>
      ( pane->attribute(tangranks_id)->pointer());
    double *ws2 = reinterpret_cast<double*>
      ( pane->attribute(weights2_id)->pointer());
    double *bcs = reinterpret_cast<double*>
      ( pane->attribute(barycrds_id)->pointer());
    int *PNids = reinterpret_cast<int*>
      ( pane->attribute(PNelemid_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;
      Vector_3<double> nrm_nz;

      // Loop through all vertices of current face.
      for ( int k=0; k<ne; ++k) {
        int windex = ene[(k+1==ne)?0:k+1]-1;

        // Only redistribute correct classification of vertices 
        if ( tranks[vindex]==rank) {
          Vector_3<double> tngs[2] = {pnts[uindex]-pnts[vindex], 
                              pnts[windex]-pnts[vindex]};

          // Compute normal of face offset
          if ( k==0 || ne>3) {
            nrm_nz = Vector_3<double>::cross_product( tngs[1], tngs[0]);
            nrm_nz.normalize();
          }
          
          // Solve for vector x, to determine maximum time step
          Vector_3<double> S[3] = { tngs[1], tngs[0], nrm_nz};
          Vector_3<double> xs; solve( S, cs[windex], xs);

          double temp = ws[windex];
          ws[windex] = std::max( (6-ne)*std::max( xs[0], xs[1]), ws[windex]);
          if(temp != ws[windex]){
            ws2[windex] = ws[windex];
            bcs[windex*2]= xs[0];
            bcs[windex*2+1] = xs[1];
            PNids[windex] = j+1;
          }
        }

        uindex=vindex; vindex=windex;
      }
    }
  }

  _wm->reduce_on_shared_nodes( w_attr, Window_manifold_2::OP_MAX);
}

Here is the call graph for this function:

Here is the caller graph for this function:

void get_redist_safe_factor ( COM::Attribute *  c_attr, COM::Attribute *  w_attr, int  rank  )

protected

void get_redist_safe_factor ( COM::Attribute *  c_attr, COM::Attribute *  w_attr, int  rank  )

protected

void get_usr_disp ( const COM::Attribute *  pmesh, COM::Attribute *  disp, double *  timestep  )

protected

Get displacement in _usr_window based on nc in _buf_window.

Use nodal coordinates in _buf_window as target and nodal coordinates in _usr_window as initial position. Determine displacement in _usr_window.

Definition at line 1001 of file Rocmop_2.C.

References COM_assertion_msg, COM_NC, Rocblas::div_scalar(), i, j, ni, nj, Rocblas::sub(), and Rocmap::update_ghosts().

{
  /*
   *    Disp array contains all nodal displacements. 
   *    For the real part, all displacements must not exceed 
   *    a given constraint. For the ghost part, they can since 
   *    they need to make up for incorrect initial positions 
   *    given by Rocflu (Rocflu never updates ghost nodes).
   *    
   *    The following steps are done.
   *
   * 1. Store displacements of REAL nodes in disp
   *
   * 2. Check if any displacements exceed the constraint
   *    2.1 If so, scale down all REAL displacements
   *
   * 3. Put updated NC in disp
   *
   * 4. Update ghost NC using Rocmap::update_ghosts(disp)
   *
   * 5. Disp now can be correctly calculated using 
   *    Rocblas::sub(disp, orig_nc, disp)
   *
   *
   *    Pornput Suriyamongkol 
   *         02/20/2007
   */
  
  std::vector<const Pane*> allusrpanes;
  std::vector<const Pane*> allbufpanes;
  _usr_window->panes(allusrpanes);
  _buf_window->panes(allbufpanes);
  COM_assertion_msg(allbufpanes.size() == allusrpanes.size(),
                    "Different number of panes on buffer and user windows.");
  
  // 1.  Store REAL displacment in disp 
  
  int disp_id = disp->id();
  for(int i = 0, ni = allusrpanes.size(); i < ni; ++i){
    
    const COM::Attribute *old_nc = allusrpanes[i]->attribute(COM::COM_NC);
    const COM::Attribute *new_nc = allbufpanes[i]->attribute(COM::COM_NC);
    COM::Attribute *p_disp       = const_cast<COM::Attribute*>
      (allusrpanes[i]->attribute(disp_id));
    
    COM_assertion_msg((p_disp->size_of_real_items() == new_nc->size_of_real_items()) &&
                      (p_disp->size_of_real_items() == old_nc->size_of_real_items()),
                      "Number of real items differs between buffer and user windows");
    
    double* old_nc_ptr = (double*)old_nc->pointer();
    double* new_nc_ptr = (double*)new_nc->pointer();
    double* p_disp_ptr = (double*)p_disp->pointer();
    
    int old_nc_stride = old_nc->stride();
    int new_nc_stride = new_nc->stride();
    int p_disp_stride = p_disp->stride();
    
    // Set distances to the next component for each data structure
    int old_nc_next_comp;
    int new_nc_next_comp;
    int p_disp_next_comp;
    
    if (old_nc_stride == 1)
      old_nc_next_comp = old_nc->size_of_items();
    else
      old_nc_next_comp = 1;
    
    if (new_nc_stride == 1)
      new_nc_next_comp = new_nc->size_of_items();
    else
      new_nc_next_comp = 1;
    
    if (p_disp_stride == 1)
      p_disp_next_comp = p_disp->size_of_items();
    else
      p_disp_next_comp = 1;
    
    // Put displacements for real nodes in p_disp
    for(int j = 0, nj = new_nc->size_of_real_items(); j < nj; ++j){
      
      *p_disp_ptr                      = *new_nc_ptr                      - *old_nc_ptr;
      *(p_disp_ptr+1*p_disp_next_comp) = *(new_nc_ptr+1*new_nc_next_comp) - *(old_nc_ptr+1*old_nc_next_comp);
      *(p_disp_ptr+2*p_disp_next_comp) = *(new_nc_ptr+2*new_nc_next_comp) - *(old_nc_ptr+2*old_nc_next_comp);
      
      // go to next item
      old_nc_ptr += old_nc_stride;
      new_nc_ptr += new_nc_stride;
      p_disp_ptr += p_disp_stride;
    }
  }
  
  // 2. Constrain displacements
  
  constrain_displacements(disp);
  
  // 3. Update real NC in disp
  
  for(int i = 0, ni = allusrpanes.size(); i < ni; ++i){
    
    const COM::Attribute *old_nc = allusrpanes[i]->attribute(COM::COM_NC);
    COM::Attribute *p_disp       = const_cast<COM::Attribute*>
      (allusrpanes[i]->attribute(disp_id));
    
    double* old_nc_ptr = (double*)old_nc->pointer();
    double* p_disp_ptr = (double*)p_disp->pointer();
    
    int old_nc_stride = old_nc->stride();
    int p_disp_stride = p_disp->stride();
    
    // Set distances to the next component for each data structure
    
    int old_nc_next_comp;
    int p_disp_next_comp;
    
    if (old_nc_stride == 1)
      old_nc_next_comp = old_nc->size_of_items();
    else
      old_nc_next_comp = 1;
    
    if (p_disp_stride == 1)
      p_disp_next_comp = p_disp->size_of_items();
    else
      p_disp_next_comp = 1;
    
    // Update real part of NC in disp
    for(int j = 0, nj = p_disp->size_of_real_items(); j < nj; ++j){
      
      *p_disp_ptr                      = *old_nc_ptr                      + *p_disp_ptr;
      *(p_disp_ptr+1*p_disp_next_comp) = *(old_nc_ptr+1*old_nc_next_comp) + *(p_disp_ptr+1*p_disp_next_comp);
      *(p_disp_ptr+2*p_disp_next_comp) = *(old_nc_ptr+2*old_nc_next_comp) + *(p_disp_ptr+2*p_disp_next_comp);
      
      // go to next item
      old_nc_ptr += old_nc_stride;
      p_disp_ptr += p_disp_stride;
    }
  }
  
  // 4. Update ghost part of disp 
  
  Rocmap::update_ghosts(disp);  
  
  // 5. Get displacment
  
  const COM::Attribute *orig_nc = pmesh->window()->attribute(COM::COM_NC);
  Rocblas::sub(disp, orig_nc, disp);
  
  if(timestep)
    Rocblas::div_scalar (disp,(const void*)timestep, disp);
}

Here is the call graph for this function:

void identify_ridge_edges ( )

protected

Identify ridge edges.

Definition at line 969 of file smooth_medial.C.

References _edges, _wrk_window, COM_DOUBLE, i, j, Element_node_enumerator::next(), nj, Vector_3< Type >::normalize(), Window_manifold_2::OP_MAX, Window_manifold_2::OP_MIN, and REAL_PANE.

{
  // Allocate buffer space for maximum and minedgev
  COM::Attribute *maxedgev = 
    _wrk_window->new_attribute( "maxedgev", 'n', COM_DOUBLE, 1, "");
  _wrk_window->resize_array( maxedgev, 0);

  COM::Attribute *minedgev = 
    _wrk_window->new_attribute( "minedgev", 'n', COM_DOUBLE, 1, "");
  _wrk_window->resize_array( minedgev, 0);
  _wrk_window->init_done();

  // Loop through the panes and its real faces
  std::vector< COM::Pane*> allpanes;
  _wrk_window->panes(allpanes);
  std::vector< COM::Pane*>::iterator it = allpanes.begin();
  Window_manifold_2::PM_iterator pm_it=_wm->pm_begin();

  const std::string att1("eigvecs"), att2("tangranks");

  int eigvecs_id = _wrk_window->attribute(att1)->id();
  int tangranks_id = _wrk_window->attribute(att2)->id();

  std::vector<std::map<Edge_ID, double> > edge_maps( allpanes.size());
  for ( int i=0, local_npanes = allpanes.size(); 
        i<local_npanes; ++i, ++it, ++pm_it) { 
    COM::Pane *pane = *it;
    std::map< Edge_ID, double> &edges_pn = edge_maps[i];

    // Obtain pointers to A_attr, bm_attr, and r_attr 
    Vector_3<double> *es = reinterpret_cast<Vector_3<double>*>
      ( pane->attribute(eigvecs_id)->pointer());
    int *tranks = reinterpret_cast<int*>
      ( pane->attribute(tangranks_id)->pointer());
    double *minvs = reinterpret_cast<double*>
      ( pane->attribute(minedgev->id())->pointer());
    double *maxvs = reinterpret_cast<double*>
      ( pane->attribute(maxedgev->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()) {
      Halfedge h( &*pm_it, Edge_ID(j+1,0), SURF::REAL_PANE), h0=h;
      
      int vindex = ene[ h.id().lid()]-1;
      do {
        int windex = ene[ h.next().id().lid()]-1;
        
        // If vertex is on a ridge
        if ( tranks[ vindex]==1) {
          int m = 1 + (tranks[windex]<2);
          Vector_3<double> tng = h.tangent(); tng.normalize();
          double feg = m * (es[3*vindex+2]* tng);

          edges_pn[ h.id()] = feg;
          if ( feg > maxvs[vindex]) maxvs[vindex] = feg;
          if ( feg < minvs[vindex]) minvs[vindex] = feg;
        }

        // Process border edge
        if ( h.is_border() && tranks[windex]==1) {
          int m = 1 + (tranks[vindex]<2);
          Vector_3<double> tng = -h.tangent(); tng.normalize();
          double feg = m * (es[3*windex+2]* tng);

          edges_pn[ h.opposite().id()] = feg;
          if ( feg > maxvs[windex]) maxvs[windex] = feg;
          if ( feg < minvs[windex]) minvs[windex] = feg;
        }

        // Update vindex
        vindex = windex;
      } while ( (h=h.next()) != h0);
    }
  }

  // Reduce on shared nodes for minedgev and maxedgev
  _wm->reduce_on_shared_nodes( minedgev, Window_manifold_2::OP_MIN);
  _wm->reduce_on_shared_nodes( maxedgev, Window_manifold_2::OP_MAX);

  it = allpanes.begin(); 
  // Insert ridge edges onto _edges
  _edges.clear(); _edges.resize( allpanes.size());
  for ( int i=0, local_npanes = allpanes.size(); i<local_npanes; ++i, ++it) { 
    COM::Pane *pane = *it;
    std::map< Edge_ID, double> &emap_pn = edge_maps[i];
    std::set< Edge_ID> &eset_pn = _edges[i];

    const double *minvs = reinterpret_cast<double*>
      ( pane->attribute(minedgev->id())->pointer());
    const double *maxvs = reinterpret_cast<double*>
      ( pane->attribute(maxedgev->id())->pointer());

    // Loop throug ebuf and put ridges edges onto edges_pn
    for ( std::map< Edge_ID, double>::const_iterator eit=emap_pn.begin(), 
            eend=emap_pn.end(); eit!=eend; ++eit) {
      Element_node_enumerator ene( pane, eit->first.eid());
      int vindex = ene[eit->first.lid()]-1;

      // Push edge onto edge
      if ( eit->second == minvs[vindex] || eit->second == maxvs[vindex])
        eset_pn.insert( eit->first);
    }
  }

  // Deallocate minedgev and maxedgev
  _wrk_window->delete_attribute( minedgev->name());
  _wrk_window->delete_attribute( maxedgev->name());
  _wrk_window->init_done();
}

Here is the call graph for this function:

void identify_ridge_edges ( )

protected

Identify ridge edges.

void identify_ridge_edges ( )

protected

Identify ridge edges.

void invert_elements ( int  conn_type )

protected

Repair inverted tets or hexes.

Definition at line 1713 of file Rocmop_2.C.

References i, MesqPane::invert(), and ni.

                                         {
#ifdef ROCPROF
  Profile_begin("Rocmop::invert_elements");
#endif
  print_legible(1,"Entering Rocmop::invert_elements");
  std::vector<Pane*> allpanes;
  _buf_window->panes(allpanes);
  for(int i=0, ni = allpanes.size(); i<ni; ++i){
    MesqPane* mp = new MesqPane(allpanes[i]);
    mp->invert();
    //mp->invert(conn_type);
    if(mp)
      delete mp;
    mp = NULL;
  }  
  print_legible(1,"Exiting Rocmop::invert_elements");
#ifdef ROCPROF
  Profile_end("Rocmop::invert_tets");
#endif
}

Here is the call graph for this function:

void invert_tets ( )

protected

Repair inverted tets.

void invert_tets ( )

protected

Repair inverted tets.

Definition at line 748 of file Rocmop.C.

References i, MesqPane::invert(), and ni.

                        {
  std::vector<Pane*> allpanes;
  _wrk_window->panes(allpanes);
  for(int i=0, ni = allpanes.size(); i<ni; ++i){
    MesqPane* mp = new MesqPane(allpanes[i]);
    mp->invert();
    if(mp)
      delete mp;
    mp = NULL;
  }  
}

Here is the call graph for this function:

void load ( const std::string &  mname )

static

Loads Rocmop onto Roccom with a given module name.

Definition at line 88 of file Rocmop.C.

References COM_METADATA, COM_new_attribute(), COM_new_window(), COM_RAWDATA, COM_set_member_function(), COM_set_object(), COM_STRING, COM_VOID, COM_window_init_done(), set_value(), smooth(), and smooth_in_place().

Referenced by COM_F_FUNC2(), and Rocmop_load_module().

                                         {

  Rocmop *mop = new Rocmop();

  COM_new_window( mname.c_str());

  std::string glb=mname+".global";

  COM_new_attribute( glb.c_str(), 'w', COM_VOID, 1, "");
  COM_set_object( glb.c_str(), 0, mop);

  COM_Type types[3];
  types[0] = COM_RAWDATA; types[1] = COM_METADATA;
  types[2] = COM_METADATA;

  COM_set_member_function( (mname+".smooth").c_str(), 
                           (Member_func_ptr)(&Rocmop::smooth),
                           glb.c_str(), "biB", types);

  COM_set_member_function( (mname+".smooth_in_place").c_str(), 
                           (Member_func_ptr)(&Rocmop::smooth_in_place),
                           glb.c_str(), "bb", types);

  types[1] = COM_STRING; types[2] = COM_VOID;
  COM_set_member_function( (mname+".set_value").c_str(), 
                           (Member_func_ptr)(&Rocmop::set_value),
                           glb.c_str(), "bii", types);    
  
  COM_window_init_done( mname.c_str());

}

Here is the call graph for this function:

Here is the caller graph for this function:

static void load ( const std::string &  mname )

static

Loads Rocmop onto Roccom with a given module name.

static void load ( const std::string &  mname )

static

Loads Rocmop onto Roccom with a given module name.

virtual SURF::Window_manifold_2* manifold ( )

inlineprotectedvirtual

Obtain a reference to the manifold.

Definition at line 201 of file Rocmop_1.h.

  { return _wm;}

virtual SURF::Window_manifold_2* manifold ( )

inlineprotectedvirtual

Obtain a reference to the manifold.

Definition at line 202 of file Rocmop.h.

  { return _wm;}

virtual SURF::Window_manifold_2* manifold ( )

inlineprotectedvirtual

Obtain a reference to the manifold.

Definition at line 223 of file Rocmop_2.h.

  { return _wm;}

void mark_elems_from_nodes	(	std::vector< std::vector< bool > > &	marked_nodes,
		std::vector< std::vector< bool > > &	marked_elems
	)

Mark the nodes which contain marked elems.

void mark_elems_from_nodes	(	std::vector< std::vector< bool > > &	marked_nodes,
		std::vector< std::vector< bool > > &	marked_elems
	)

Mark the nodes which contain marked elems.

Definition at line 712 of file Rocmop.C.

References i, j, k, and nk.

                                                                             {
  if(_verb) 
    std::cout << "Entering Rocmop::mark_elems_from_nodes" << std::endl;

  std::vector<const COM::Pane*> allpanes;
  _wrk_window->panes(allpanes);
  int local_npanes = (int)allpanes.size();

  marked_elems.clear();
  marked_elems.resize(local_npanes);
  
  //Loop through panes
  for(int i=0; i < (int)(local_npanes); ++i){

    marked_elems[i].clear();
    marked_elems[i].resize(allpanes[i]->size_of_real_elements(),false);

    // Loop through real elements.
    // Mark for quality check if they contain shared nodes.
    int s_real_elems = allpanes[i]->size_of_real_elements();
    std::vector<int> nodes;
    for(int j=1; j<= s_real_elems; ++j){
      Element_node_enumerator ene(allpanes[i],j);
      ene.get_nodes(nodes);
      for(int k=0, nk=nodes.size(); k<nk; ++k){
        if (marked_nodes[i][nodes[k]-1])
          marked_elems[i][j-1] = true;
      }
    }
  }

  if(_verb > 1) 
    std::cout << "Exiting Rocmop::mark_elems_from_nodes" << std::endl;  
}

void mark_elems_from_nodes	(	std::vector< std::vector< bool > > &	marked_nodes,
		std::vector< std::vector< bool > > &	marked_elems
	)

Mark the nodes which contain marked elems.

void obtain_extremal_dihedrals ( const COM::Attribute *  att, double *  min, double *  max  )

protected

Obtain the min and max dihedral angles.

Definition at line 1991 of file Rocmop_2.C.

References COMMPI_Initialized(), Angle_Metric_3::compute(), i, Geo_Metric_Base_3::initialize(), k, MPI_MAX, MPI_MIN, ni, nk, and rank.

                                                                                                   {

  const COM::Window *win = att->window();

  bool cominit = COMMPI_Initialized();

  int rank =0;
  if(cominit){
    int ierr = MPI_Comm_rank( win->get_communicator(),
                              &rank); assert( ierr == 0);
  }

  std::vector<const Pane*> allpanes;
  win->panes(allpanes);

  max_angle[0] = 0.0;
  min_angle[0] = 180.0;
  double angles[] = {0.0, 0.0};

  for(int i=0,ni = allpanes.size(); i<ni; ++i){
    for(int k =0,nk=allpanes[i]->size_of_real_elements(); k<nk; ++k){
      Element_node_enumerator ene(allpanes[i],k+1);
      Angle_Metric_3 am;
      am.initialize(ene);
      am.compute(angles);

      if(angles[1]>max_angle[0])
        max_angle[0] = angles[1];
      if(angles[0]<min_angle[0])
        min_angle[0] = angles[0];
    }
  }

  if(cominit){
    double temp = max_angle[0];
    MPI_Allreduce(&max_angle[0], &temp, 1, MPI_DOUBLE, MPI_MAX,
                  win->get_communicator());
    max_angle[0] = temp;

    temp = min_angle[0];
    MPI_Allreduce(&min_angle[0], &temp, 1, MPI_DOUBLE, MPI_MIN,
                  win->get_communicator());
    min_angle[0] = temp;
  }
}

Here is the call graph for this function:

void perform_iterative_smoothing ( )

protected

Perform iterative smoothing.

Run an iterative smoothing method if needed until either the maximum number of iterations is reached or the required quality is achieved.

void perform_iterative_smoothing ( )

protected

Perform iterative smoothing.

Run an iterative smoothing method if needed until either the maximum number of iterations is reached or the required quality is achieved.

Definition at line 798 of file Rocmop_1.C.

References COM_assertion_msg.

                                        {
  // All internal iterative behavior has been removed.
#ifdef ROCPROF
  Profile_begin("Rocmop::perform_itersmooth");
#endif

  print_legible(1,"    Entering Rocmop::perform_iterative_smoothing");

  COM_assertion_msg(_buf_window, "Unexpected NULL pointer encountered.");

  std::vector<const Pane*> allpanes;
  _buf_window->panes(allpanes);
  
#ifdef MESQUITE

  if(_method==SMOOTH_VOL_MESQ_WG)
    smooth_vol_mesq_wg();

  else if(_method==SMOOTH_VOL_MESQ_NG)
    smooth_vol_mesq_ng();
  
#else

  if((_method==SMOOTH_VOL_MESQ_WG) || (_method==SMOOTH_VOL_MESQ_NG))
    COM_assertion_msg(0,"Rocmop not compiled with MESQUITE");

#endif
  
  else if(_method==SMOOTH_SURF_MEDIAL)
    smooth_surf_medial();

  else COM_assertion_msg(0, "No valid iterative smoothing method selected");  
  
  print_legible(1,"    Exiting Rocmop::perform_iterative_smoothing");
#ifdef ROCPROF
  Profile_end("Rocmop::perform_itersmooth");
#endif
}

void perform_iterative_smoothing ( double  pre_quality )

protected

Perform iterative smoothing.

Run an iterative smoothing method if needed until either the maximum number of iterations is reached or the required quality is achieved.

Parameters

pre_quality

see perform_smoothing()

Definition at line 290 of file Rocmop.C.

References COM_assertion_msg, i, and MPI_MIN.

                                                          {
  if(_verb) 
    std::cout << "    Entering Rocmop::perform_iterative_smoothing" << std::endl;

  COM_assertion_msg(_wrk_window, "Unexpected NULL pointer encountered.");

  // bool_iter is true until the max number of iterations (_niter)
  // is reached or the minimum quality (_tol) is reached.
  int to_iter = true; 
  std::vector<const Pane*> allpanes;
  _wrk_window->panes(allpanes);
  double cur_qual = pre_quality;
  
  for(int i = 0; 
      ( (i < _niter) && (to_iter) ); 
      ++i, agree_int(to_iter, MPI_MIN)){

    if(_verb > 2)
      std::cout << "      Smoothing iteration " << i << std::endl;

#ifdef MESQUITE
    if(_method==SMOOTH_VOL_MESQ_WG){
      //std::string msg("\nQuality prior to smoothing = ");
      //print_quality(msg);
      smooth_vol_mesq_wg(cur_qual);
    }
    //    else if(_method == SMOOTH_LAPLACE)
    //  smooth_laplace();
    else if(_method == SMOOTH_VOL_MESQ_NG){
      //std::string msg("Quality prior to smoothing = ");
      //print_quality(msg);
      smooth_vol_mesq_ng(cur_qual);
    }
#else
    if((_method==SMOOTH_VOL_MESQ_WG) || (_method==SMOOTH_VOL_MESQ_NG))
      COM_assertion_msg(0,"Rocmop not compiled with MESQUITE");
#endif
    else if(_method==SMOOTH_SURF_MEDIAL)
      smooth_surf_medial();
    else COM_assertion_msg(0, "No valid iterative smoothing method selected");  

    // If a non zero quality tolerance is set, then determine if quality
    // is low enough to warrant smoothing.
    if( _tol != 0.0 ){
      cur_qual = check_all_elem_quality(allpanes);
      if(cur_qual >= _tol)
        to_iter = false;
    }
  }

  if(_verb > 1) 
    std::cout << "    Exiting Rocmop::perform_iterative_smoothing" << std::endl;
}

void perform_noniterative_smoothing ( )

protected

Perform noniterative smoothing.

void perform_noniterative_smoothing ( )

protected

Perform noniterative smoothing.

void perform_noniterative_smoothing ( )

protected

Perform noniterative smoothing.

Definition at line 344 of file Rocmop.C.

References COM_assertion_msg.

                                           {
  if(_verb) 
    std::cout << "    Entering Rocmop::perform_noniterative_smoothing" << std::endl;

  if(_niter != 1){
    std::cerr << "Although the maximum number of iterations > 1 is selected,\n"
              << "the smoothing method is noniterative, and will run once.\n";
  }

  if (0)
    ;// Currently, no noniterative methods exist.
  else COM_assertion_msg(0, "No valid noniterative smoothing method selected");

  if(_verb > 1) 
    std::cout << "    Exiting Rocmop::perform_noniterative_smoothing" << std::endl;
}

void perform_smoothing ( )

protected

Perform smoothing on _buf_window.

Calls smoother_specific_init(), and then selects between perform_iterative_smoothing() and perform_noniterative_smoothing() as required by the currently selected smoothing method.

void perform_smoothing ( double  pre_quality )

protected

Perform smoothing on _wrk_window.

Calls smoother_specific_init(), and then selects between perform_iterative_smoothing() and perform_noniterative_smoothing() as required by the currently selected smoothing method.

Parameters

pre_quality

pre-smoothed mesh quality or the worst quality possible if not in lazy or monotone mode.

Definition at line 269 of file Rocmop.C.

References COM_assertion_msg.

                                                {
  if(_verb) 
    std::cout << "  Entering Rocmop::perform_smoothing" << std::endl;

  COM_assertion_msg(_wrk_window, "Unexpected NULL pointer encountered.");

  // Perform initialization specific to the smoother.
  smoother_specific_init();

  // Select iterative or noniterative option depending on the smoother.
  if (0) // No noniterative smoothers currently
    perform_noniterative_smoothing();
  else if( _method < SMOOTH_NONE)
    perform_iterative_smoothing(pre_quality);
  else
    COM_assertion_msg(0, "No valid smoothing method selected");

  if(_verb) 
    std::cout << "  Exiting Rocmop::perform_smoothing" << std::endl;
}

void perform_smoothing ( )

protected

Perform smoothing on _buf_window.

Calls smoother_specific_init(), and then selects between perform_iterative_smoothing() and perform_noniterative_smoothing() as required by the currently selected smoothing method.

Definition at line 776 of file Rocmop_1.C.

References COM_assertion_msg.

                              {
#ifdef ROCPROF
  Profile_begin("Rocmop::perform_smoothing");
#endif
  print_legible(1,"  Entering Rocmop::perform_smoothing");

  COM_assertion_msg(_buf_window, "Unexpected NULL pointer encountered.");

  // Select iterative or noniterative option depending on the smoother.
  if (0) // No noniterative smoothers currently
    perform_noniterative_smoothing();
  else if( _method < SMOOTH_NONE)
    perform_iterative_smoothing();
  else
    COM_assertion_msg(0, "No valid smoothing method selected");

  print_legible(1,"  Exiting Rocmop::perform_smoothing");
#ifdef ROCPROF
  Profile_end("Rocmop::perform_smoothing");
#endif
}

void perturb_stationary ( )

protected

Randomly perturn stationary nodes on pane boundary, not on phys. surface.

Definition at line 1689 of file Rocmop_1.C.

References COM_assertion_msg, COM_NC, COMMPI_Initialized(), Angle_Metric_3::compute(), i, Geo_Metric_Base_3::initialize(), j, k, Mesquite::length(), MPI_MAX, ni, nj, nk, and cimg_library::cimg::rand().

                               {
#ifdef ROCPROF
  Profile_begin("Rocmop::perturb_stat");
#endif
  // first calculate current quality of all nodes.
  std::vector<const Pane*> allpanes;
  std::vector<const Pane*> allpanes_usr;
  _buf_window->panes(allpanes);
  _usr_window->panes(allpanes_usr);

  COM::Attribute *w_fgpn_bnd = 
    _buf_window->attribute("is_fgpn_bnd");
  COM::Attribute *w_qual_b_perturb =
    _buf_window->attribute("qual_b_perturb");
  COM::Attribute *w_qual_a_perturb =
    _buf_window->attribute("qual_a_perturb");

  double angles[] = {0.0, 0.0};
  
  // Get worst quality on local panes
  for(int i=0,ni=allpanes.size(); i<ni; ++i){
  
    COM::Attribute *p_fgpn_bnd = 
      const_cast<COM::Attribute*>(allpanes[i]->attribute(w_fgpn_bnd->id()));
    COM::Attribute *p_qual_b_perturb = 
      const_cast<COM::Attribute*>(allpanes[i]->attribute(w_qual_b_perturb->id()));
 
    int* fgpn_bnd_ptr = reinterpret_cast<int*>(p_fgpn_bnd->pointer());
    double* qual_b_perturb_ptr = 
      reinterpret_cast<double*>(p_qual_b_perturb->pointer());

    std::vector<int> elist;
    // only worry about real nodes, ghost nodes should be fixed      
    for(int j=0, nj= allpanes[i]->size_of_real_nodes(); j<nj; ++j){
      if(fgpn_bnd_ptr[j]){
        _dcs[i]->incident_elements(j+1,elist);
        qual_b_perturb_ptr[j] = 0.0;
        for(uint k=0, nk=elist.size(); k<nk; ++k){
          Element_node_enumerator ene(allpanes[i],k+1);
          Angle_Metric_3 am;
          am.initialize(ene);
          am.compute(angles);
          if(angles[1]>qual_b_perturb_ptr[j])
            qual_b_perturb_ptr[j] = angles[1];  
        }
      }
    }
  }

  // Get the worst quality across all panes
  if(COMMPI_Initialized()){
    MAP::Pane_communicator pc(_buf_window, _buf_window->get_communicator());
    pc.init(w_qual_b_perturb);
    pc.begin_update_shared_nodes();
    pc.reduce_on_shared_nodes(MPI_MAX);
    pc.end_update_shared_nodes();
  }    

  // Now generate random perturbations
  for(int i=0,ni=allpanes.size(); i<ni; ++i){
  
    COM::Attribute *p_fgpn_bnd = 
      const_cast<COM::Attribute*>(allpanes[i]->attribute(w_fgpn_bnd->id()));
    COM::Attribute *p_nc = 
      const_cast<COM::Attribute*>(allpanes[i]->attribute(COM::COM_NC));
 
    int* fgpn_bnd_ptr = 
      reinterpret_cast<int*>(p_fgpn_bnd->pointer());
    Vector_3<double>* nc_ptr = 
      reinterpret_cast<Vector_3<double>*>(p_nc->pointer());

    std::vector<int> elist;
    // only worry about real nodes, ghost nodes should be fixed      
    for(int j=0, nj= (int)allpanes[i]->size_of_real_nodes(); j<nj; ++j){

      if(fgpn_bnd_ptr[j]){
        _dcs[i]->incident_elements(j+1,elist);
        
        if(_verb > 1)
          std::cout << "Rocmop: Perturbing node " << j+1 << std::endl;
        
        //select a random element
        int rand_el = (std::rand()%elist.size());
        Element_node_enumerator ene(allpanes[i],elist[rand_el]);
        std::vector<int> enodes;
        ene.get_nodes(enodes);
        
        // get length of a randomly selected incident edge of the element
        int rand_ed = (std::rand()%3)+1;
        int nindex = 0;
        for(int nk= (int)enodes.size(); nindex<nk; ++nindex){
          if(enodes[nindex]==j+1)        
            break;
        }      
        
        COM_assertion_msg((nindex>=0 && nindex<= (int)enodes.size()),
                          "Node not found in supposedly adjacent element\n");
        double length = (nc_ptr[j]-nc_ptr[enodes[(nindex+rand_ed)%4]-1]).norm();
        
        Vector_3<double> perturb(length,length,length);
        for(int k=0; k<3; ++k){
          perturb[0] *= (std::rand()%2 == 0) ?
            (double)((std::rand()%1000)+1)/1000.0 :
            -1.0*((double)((std::rand()%1000)+1)/1000.0);         
        }
        nc_ptr[j] += perturb;
      }
    }
  }

    // Get perturbed worst quality on local panes
  for(int i=0,ni=(int)allpanes.size(); i<ni; ++i){
  
    COM::Attribute *p_fgpn_bnd = 
      const_cast<COM::Attribute*>(allpanes[i]->attribute(w_fgpn_bnd->id()));
    COM::Attribute *p_qual_a_perturb = 
      const_cast<COM::Attribute*>(allpanes[i]->attribute(w_qual_a_perturb->id()));
 
    int* fgpn_bnd_ptr = 
      reinterpret_cast<int*>(p_fgpn_bnd->pointer());
    double* qual_a_perturb_ptr = 
      reinterpret_cast<double*>(p_qual_a_perturb->pointer());

    std::vector<int> elist;
    // only worry about real nodes, ghost nodes should be fixed      
    for(int j=0, nj= (int)allpanes[i]->size_of_real_nodes(); j<nj; ++j){
      if(fgpn_bnd_ptr[j]){
        _dcs[i]->incident_elements(j+1,elist);
        qual_a_perturb_ptr[j] = 0.0;
        for(int k=0, nk=(int)elist.size(); k<nk; ++k){
          Element_node_enumerator ene(allpanes[i],k+1);
          Angle_Metric_3 am;
          am.initialize(ene);
          am.compute(angles);
          if(angles[1]>qual_a_perturb_ptr[j])
            qual_a_perturb_ptr[j] = angles[1];  
        }
      }
    }
  }

  // Get the perturbed worst quality across all panes
  if(COMMPI_Initialized()){
    MAP::Pane_communicator pc(_buf_window, _buf_window->get_communicator());
    pc.init(w_qual_a_perturb);
    pc.begin_update_shared_nodes();
    pc.reduce_on_shared_nodes(MPI_MAX);
    pc.end_update_shared_nodes();
  }    

  // Reset positions of any nodes whose worst adj. element quality
  // has gotten worse.
  for(int i=0,ni=allpanes.size(); i<ni; ++i){
    
    COM::Attribute *p_fgpn_bnd = 
      const_cast<COM::Attribute*>(allpanes[i]->attribute(w_fgpn_bnd->id()));
    COM::Attribute *p_nc = 
      const_cast<COM::Attribute*>(allpanes[i]->attribute(COM::COM_NC));
    COM::Attribute *p_nc_usr = 
      const_cast<COM::Attribute*>(allpanes_usr[i]->attribute(COM::COM_NC));
      
    COM::Attribute *p_qual_b_perturb = 
      const_cast<COM::Attribute*>(allpanes[i]->attribute(w_qual_b_perturb->id()));
    COM::Attribute *p_qual_a_perturb = 
      const_cast<COM::Attribute*>(allpanes[i]->attribute(w_qual_a_perturb->id()));
 
    int* fgpn_bnd_ptr = 
      reinterpret_cast<int*>(p_fgpn_bnd->pointer());
    Vector_3<double>* nc_ptr = 
      reinterpret_cast<Vector_3<double>*>(p_nc->pointer());
    Vector_3<double>* nc_ptr_usr = 
      reinterpret_cast<Vector_3<double>*>(p_nc_usr->pointer());
    double* qual_b_perturb_ptr = 
      reinterpret_cast<double*>(p_qual_b_perturb->pointer());
    double* qual_a_perturb_ptr = 
      reinterpret_cast<double*>(p_qual_a_perturb->pointer());

    // only worry about real nodes, ghost nodes should be fixed      
for(int j=0, nj= (int)allpanes[i]->size_of_real_nodes(); j<nj; ++j){
      if(fgpn_bnd_ptr[j] &&
         (qual_a_perturb_ptr[j] > qual_b_perturb_ptr[j]))
        nc_ptr[j] = nc_ptr_usr[j];
    }
  }  
#ifdef ROCPROF
  Profile_end("Rocmop::perturb_stat");
#endif
}

Here is the call graph for this function:

void perturb_stationary ( )

protected

Randomly perturn stationary nodes on pane boundary, not on phys. surface.

void print_extremal_dihedrals ( COM::Window *  window )

protected

Print the min and max dihedral angles along with their locations.

Definition at line 1463 of file Rocmop_1.C.

References COMMPI_Initialized(), Angle_Metric_3::compute(), i, Geo_Metric_Base_3::initialize(), k, ni, nk, and rank.

                                                       {
#ifdef ROCPROF
  Profile_begin("Rocmop::ext_dihedrals");
#endif

  // Find the rank of this processor, making sure to only make
  // MPI calls on parallel runs
  int myrank =0;
  if(COMMPI_Initialized()){ // true if in parallel
    int ierr = MPI_Comm_rank( window->get_communicator(),
                              &myrank); 
    assert( ierr == 0);
  }

  // Calculate extremal angles on local panes, and
  // record their locations by pane id and element id

  double max_angle = 0.0;
  double min_angle = 180.0; 
  int max_pane = -1; // id of pane w/ the largest dihedral angle
  int min_pane = -1; // id of pane w/ the smallest dihedral angle
  int max_elem = -1; // id of element w/ the largest dihedral angle
  int min_elem = -1; // id of element w/ the smallest dihedral angle
  double angles[] = {0.0, 0.0};
    
  std::vector<Pane*> allpanes;
  window->panes(allpanes);
  
  for(int i=0,ni = allpanes.size(); i<ni; ++i){
    for(int k =0,nk=allpanes[i]->size_of_real_elements(); k<nk; ++k){
      Element_node_enumerator ene(allpanes[i],k+1);
      Angle_Metric_3 am;
      am.initialize(ene);
      am.compute(angles);
      if(angles[1]>max_angle){
        max_angle = angles[1];
        max_pane = allpanes[i]->id();
        max_elem = k+1;
      }
      if(angles[0]<min_angle){
        min_angle = angles[0];
        min_pane = allpanes[i]->id();
        min_elem = k+1;
      }
    }
  }

  // Find minimum angle across all panes.

  // create MPI_DOUBLE_INT data types for send/recv buffs
  struct {
    double value;
    int rank;
  } send, recv;
  send.value = min_angle;
  send.rank = myrank;

  if(COMMPI_Initialized()){
    MPI_Allreduce(&send,&recv,1,MPI_DOUBLE_INT,MPI_MINLOC,
               window->get_communicator());
  }
    if(recv.rank == myrank && _verb > 1){
      std::cout << "Rocmop:" << std::endl << "Rocmop: " 
                << std::setw(10) << min_angle << " on element "
                << min_elem << " of pane " << min_pane 
                << "." << std::endl;
    }
    
  // Sync. output.
  if(COMMPI_Initialized()){
    MPI_Barrier(window->get_communicator());
  }

  // Find maximum angle across all panes.

  send.value = max_angle;

  if(COMMPI_Initialized()){
    MPI_Allreduce(&send,&recv,1,MPI_DOUBLE_INT,MPI_MAXLOC,
               window->get_communicator());
  }
  if(recv.rank == myrank && _verb > 1){
    std::cout << "Rocmop:" << std::endl << "Rocmop: " 
              << std::setw(10) << max_angle << " on element "
              << max_elem << " of pane " << max_pane << std::endl;
  }
  // Sync. output.
  if(COMMPI_Initialized()){
    MPI_Barrier(window->get_communicator());
  }
#ifdef ROCPROF
  Profile_end("Rocmop::ext_dihedrals");
#endif

}

Here is the call graph for this function:

void print_extremal_dihedrals ( COM::Window *  window )

protected

Print the min and max dihedral angles along with their locations.

void print_legible ( int  verb, const char *  msg  )

protected

Single process print message if given verbosity level is exceeded.

Definition at line 1408 of file Rocmop_1.C.

References COMMPI_Initialized(), and rank.

                                                   {
  if(_verb > verb){

  int rank =0;

  if(COMMPI_Initialized()){
    int ierr = MPI_Comm_rank( _buf_window->get_communicator(),
                              &rank); 
    assert( ierr == 0);
  }

  if (rank==0)
    std::cout << "Rocmop: " << msg << std::endl;
  }
}

Here is the call graph for this function:

void print_legible ( int  verb, const char *  msg  )

protected

Single process print message if given verbosity level is exceeded.

void print_mquality ( std::string &  s, std::vector< std::vector< bool > > &  to_check  )

protected

Print the quality range of marked elements, for debugging.

Definition at line 883 of file Rocmop.C.

References COM_assertion_msg, COMMPI_Initialized(), Angle_Metric_3::compute(), i, Geo_Metric_Base_3::initialize(), k, MPI_MAX, MPI_MIN, ni, nk, and rank.

                                                                  {

  std::vector<std::vector<bool> > elem_to_check;
  mark_elems_from_nodes(to_check, elem_to_check);

  std::vector<Pane*> allpanes;
  _wrk_window->panes(allpanes);
  
  double max_angle = 0.0;
  double min_angle = 180.0;
  double angles[] = {0.0, 0.0};
  int id = -1;
  for(int i=0,ni = allpanes.size(); i<ni; ++i){
    for(int k =0,nk = elem_to_check[i].size(); k<nk; ++k){
      if(elem_to_check[i][k]){
        Element_node_enumerator ene(allpanes[i],k+1);
        Angle_Metric_3 am;
        am.initialize(ene);
        am.compute(angles);
        if(angles[1]>max_angle)
          max_angle = angles[1];
        if(angles[0]<min_angle){
          id = k;
          min_angle = angles[0];
        }
      }
    }
  }
  int rank =0;
  double temp = min_angle;
  if(COMMPI_Initialized()){
    agree_double(max_angle,MPI_MAX);
    agree_double(min_angle,MPI_MIN);
    int ierr = MPI_Comm_rank( _wrk_window->get_communicator()
                              , &rank); 
    assert( ierr == 0);
  }  

  if (rank==0){
    string outstr("angles.txt");
    ofstream file (outstr.c_str(), std::ios::app);
    COM_assertion_msg( file, "File failed to open\n");

    file << std::left << std::setw(30) << s << std::setw(0)
         << "(" << min_angle << " , " << max_angle << ")";

    file.close();
  }

  if(COMMPI_Initialized())
    int ierr = MPI_Barrier(_wrk_window->get_communicator());

  if(min_angle == temp){
    
    string outstr("angles.txt");
    ofstream file (outstr.c_str(), std::ios::app);
    COM_assertion_msg( file, "File failed to open\n");
    
    file << "  worst = (" << rank << " , " << id << ")\n";
    file.close();
  }

  if(COMMPI_Initialized())
    int ierr = MPI_Barrier(_wrk_window->get_communicator());
}

Here is the call graph for this function:

void print_mquality ( std::string &  s, std::vector< std::vector< bool > > &  to_check  )

protected

Print the quality range of marked elements, for debugging.

void print_mquality ( std::string &  s, std::vector< std::vector< bool > > &  to_check  )

protected

Print the quality range of marked elements, for debugging.

void print_quality ( std::string &  s )

protected

Print the quality range of all elements, for debugging.

Definition at line 846 of file Rocmop.C.

References COM_assertion_msg, COMMPI_Initialized(), Angle_Metric_3::compute(), i, Geo_Metric_Base_3::initialize(), k, MPI_MAX, MPI_MIN, ni, nk, and rank.

                                      {

  string outstr("angles.txt");
  ofstream file (outstr.c_str(), std::ios::app);
  COM_assertion_msg( file, "File failed to open\n");
  
  std::vector<Pane*> allpanes;
  _wrk_window->panes(allpanes);

  double max_angle = 0.0;
  double min_angle = 180.0;
  double angles[] = {0.0, 0.0};
  for(int i=0,ni = allpanes.size(); i<ni; ++i){
    for(int k =0,nk=allpanes[i]->size_of_real_elements(); k<nk; ++k){
      Element_node_enumerator ene(allpanes[i],k+1);
      Angle_Metric_3 am;
      am.initialize(ene);
      am.compute(angles);
      if(angles[1]>max_angle)
        max_angle = angles[1];
      if(angles[0]<min_angle)
        min_angle = angles[0];
    }
  }
  int rank =0;
  if(COMMPI_Initialized()){
    agree_double(max_angle,MPI_MAX);
    agree_double(min_angle,MPI_MIN);
    int ierr = MPI_Comm_rank( _wrk_window->get_communicator(),
                              &rank); assert( ierr == 0);
  }
  if (rank==0){
    file << std::left << std::setw(30) << s << std::setw(0)
         << "(" << min_angle << " , " << max_angle << ")\n";
  }
}

Here is the call graph for this function:

void print_quality ( std::string &  s, const std::string &  s2  )

protected

Print the quality range of all elements, for debugging.

Definition at line 1559 of file Rocmop_1.C.

References COM_assertion_msg, COMMPI_Initialized(), Angle_Metric_3::compute(), i, Geo_Metric_Base_3::initialize(), k, MPI_MAX, MPI_MIN, ni, nk, and rank.

                                                          {

  string outstr("angles.txt");


  ofstream file;
  ofstream file2;

  bool cominit = COMMPI_Initialized();

  int rank =0;
  if(cominit){
    int ierr = MPI_Comm_rank( _buf_window->get_communicator(),
                              &rank); assert( ierr == 0);
  }

  std::vector<Pane*> allpanes;
  _buf_window->panes(allpanes);

  double max_angle = 0.0;
  double min_angle = 180.0;
  double angles[] = {0.0, 0.0};

  file2.open(s2.c_str());
  for(int i=0,ni = allpanes.size(); i<ni; ++i){
    for(int k =0,nk=allpanes[i]->size_of_real_elements(); k<nk; ++k){
      Element_node_enumerator ene(allpanes[i],k+1);
      Angle_Metric_3 am;
      am.initialize(ene);
      am.compute(angles);
      file2 << angles[0] << " " << angles[1] << std::endl;
      if(angles[1]>max_angle)
        max_angle = angles[1];
      if(angles[0]<min_angle)
        min_angle = angles[0];
    }
  }
  file2.close();

  agree_double(max_angle,MPI_MAX);
  agree_double(min_angle,MPI_MIN);

  if (rank==0){
    file.open(outstr.c_str(), std::ios::app);
    COM_assertion_msg( file, "File failed to open\n");
    file << std::left << std::setw(30) << s << std::setw(0)
         << "(" << min_angle << " , " << max_angle << ")\n";
    file.close();
  }
}

Here is the call graph for this function:

void print_quality ( std::string &  s, const std::string &  s2  )

protected

Print the quality range of all elements, for debugging.

void read_config_file

(

const std::string &

cfname

)

Definition at line 122 of file Rocmop_1.C.

References COMMPI_Initialized(), get_next_line(), MPI_COMM_WORLD, and rank.

{
  std::ifstream Inf;
  Inf.open(cfname.c_str());

  // We can't use the usual print_legible function here because the
  // buffer winodw may not be created yet.
  int rank =0;
  if(COMMPI_Initialized()){
    int ierr = MPI_Comm_rank(MPI_COMM_WORLD,&rank); 
    assert( ierr == 0);
  }

  // If the file doesn't exist, fire off a warning and continue
  if(!Inf){
    if(rank == 0)
      std::cerr << "Rocmop: Warning: Could not open configfile, " 
                << cfname << "." << std::endl;
    return;
  }

  std::string line;
  std::istringstream Istr;

  // Read and set the verbosity
  int verbosity;
  get_next_line(Inf,line,'#');
  //     At this point, "line" should have data in it.  By using the 
  //     istringstream object this way, we eliminate all whitespace
  //     and extract only the valid part of the line.  If the file 
  //     has already closed, line will be empty, and verbosity will
  //     usually take on some nonsense value.  
  Istr.str(line);
  Istr >> verbosity;
  //     If the line had only the one chunk of data in it, Istr will 
  //     have a bit set that needs to be cleared with a call to "clear()"
  //     before populating it with another string.
  Istr.clear();
  if(verbosity < 0 || verbosity > 10){
    std::cerr << "Rocmop: Warning: Invalid verbosity from configfile. " 
              << "Giving up on " << cfname << "." << std::endl;
    Inf.close();
    return;
  }
  _verb = verbosity;
  // Speak up if verbosity is turned on!
  if (rank==0 && _verb){
    std::cout << "Rocmop: Found configfile, " << cfname << "." << std::endl
              << "Rocmop: Setting verbosity to " << _verb << "." << std::endl;
  }  

  // Load the next line from configfile and set smoothing method
  get_next_line(Inf,line,'#');
  Istr.str(line);
  Istr >> _method;
  Istr.clear();
  if(_method < 0 || _method > SMOOTH_NONE){
    std::cerr << "Rocmop: Warning: Invalid smoothing method from configfile. " 
              << "Giving up on " << cfname << "." << std::endl;
    Inf.close();
    return;
  }
  // I love these ternary operators
  if (rank==0 && _verb){
    std::cout << "Rocmop: Setting method to "
              << (_method == SMOOTH_VOL_MESQ_WG ? "SMOOTH_VOL_MESQ_WG" :
                  (_method == SMOOTH_VOL_MESQ_NG ? "SMOOTH_VOL_MESQ_NG" :
                   (_method == SMOOTH_SURF_MEDIAL ? "SMOOTH_SURF_MEDIAL" :
                    "SMOOTH_NONE"))) << "." << std::endl;
  }

  // Load the next line and set laziness
  int lazy;
  get_next_line(Inf,line,'#');
  Istr.str(line);
  Istr >> lazy;
  Istr.clear();
  if(lazy > 0) 
    _lazy = 1;
  if (rank==0 && _verb){
    std::cout << "Rocmop: Setting lazy to " << _lazy << "." << std::endl;
  }  

  // Load the next line and set tolerance for lazy threshold
  float tol;
  get_next_line(Inf,line,'#');
#if ! NO_STRINGSTREAM
  Istr.str(line);
  Istr >> tol;
  Istr.clear();
#else
  sscanf(line.c_str(), "%f", &tol);
#endif
  if(tol < 0. || tol > 180.){
    std::cerr << "Rocmop: Warning: Invalid dihedral angle tolerance"
              << " from configfile. " 
              << "Giving up on " << cfname << "." << std::endl;
    Inf.close();
    return;
  }
  _tol = tol;
  if (rank==0 && _verb){
    std::cout << "Rocmop: Setting tolerance to " << _tol << "." << std::endl;
  }  

  // Load the next line and set node displacement constraint for Mesquite smoothing 
  float max_disp;
  get_next_line(Inf,line,'#');
#if ! NO_STRINGSTREAM
  Istr.str(line);
  Istr >> max_disp;
  Istr.clear();
#else
  sscanf(line.c_str(), "%f", &max_disp);
#endif
  if(max_disp < 0. || max_disp > 10.0){
    std::cerr << "Rocmop: Warning: Invalid displacement constraint from configfile. " 
              << "Giving up on " << cfname << "." << std::endl;
    Inf.close();
    return;
  }
  _maxdisp = max_disp;
  if (rank==0 && _verb){
    std::cout << "Rocmop: Setting displacement constraint to " 
              << _maxdisp << "." << std::endl;
  }  

  // Load the next line and set N, which forces Rocmop to smooth every Nth step
  get_next_line(Inf,line,'#');
  Istr.str(line);
  Istr >> _smoothfreq;
  Istr.clear();
  if(_smoothfreq <= 0 || _method == SMOOTH_NONE)
    _smoothfreq = 0;
  if (rank==0 && _verb){
    if(_method == SMOOTH_NONE){
      std::cout << "Rocmop: No method selected, setting N to 0"
                << ", disabling smoothing." << std::endl;
    }
    else{
      std::cout << "Rocmop: Setting N to " << _smoothfreq 
                << (_smoothfreq==0 ? ", disabling smoothing." : ".") 
                << std::endl;
    }
  }  
  
  // Displacement threshold.  When the max displacement among all processors 
  // exceeds this value, smooth.
  get_next_line(Inf,line,'#');
#if ! NO_STRINGSTREAM
  Istr.str(line);
  Istr >> _disp_thresh;
  Istr.clear();
#else
  sscanf(line.c_str(), "%f", &_disp_thresh);
#endif
  if(_disp_thresh < 0.0)
    _disp_thresh = 0.0;
  else if ((_smoothfreq > 1) && (_disp_thresh > 0.0) ){
    _smoothfreq = 1;
    if(rank==0 && _verb)
      std::cout << "Rocmop: WARNING: N reset to 1 to enable displacement thresholding."
                << std::endl;
  }
  if (rank==0 && _verb){
    std::cout << "Rocmop: Setting displacement threshold to " 
              << _disp_thresh << "." << std::endl;
  }  
  
  // Close up the file, we're done.
  Inf.close();
  
  N = _smoothfreq-1;          // initial value used in smooth() function
}

Here is the call graph for this function:

void read_config_file

(

const std::string &

)

void redistribute_vertices_ridge ( )

protected

Redistribute ridge vertices within their tangent spaces.

Definition at line 547 of file smooth_medial.C.

References _edges, _wrk_window, COM_NC, Rocblas::copy_scalar(), get_redist_safe_factor(), i, min(), Window_manifold_2::OP_SUM, Element_node_enumerator::size_of_edges(), and swap().

Referenced by smooth_surf_medial().

                              {
  // Use the following buffer spaces: _eigvalues for vector c (only first one)

  const std::string att1("disps"), att2("lambda"), att3("weights"),
    att4("cntnvecs"), att5("cntnranks");
  int disps_id = _wrk_window->attribute(att1)->id();
  // c_attr are the edge centers
  // w_attr is the nodal safe factor 
  COM::Attribute *c_attr = _wrk_window->attribute(att2);
  COM::Attribute *w_attr = _wrk_window->attribute(att3);
  int cntnvecs_id = _wrk_window->attribute(att4)->id();
  int cntnranks_id = _wrk_window->attribute(att5)->id();

  double zero = 0.;
  Rocblas::copy_scalar( &zero, c_attr);

  // Loop through the panes and its real faces
  std::vector< COM::Pane*> allpanes;
  _wrk_window->panes(allpanes);
  std::vector< COM::Pane*>::iterator it = allpanes.begin();
  Window_manifold_2::PM_iterator pm_it=_wm->pm_begin();
  for ( int i=0, local_npanes = allpanes.size(); 
        i<local_npanes; ++i, ++it, ++pm_it) { 
    COM::Pane *pane = *it;
    
    Vector_3<double> *es = reinterpret_cast<Vector_3<double>*>
      ( pane->attribute(cntnvecs_id)->pointer());
    const Point_3<double> *pnts = reinterpret_cast<Point_3<double>*>
      (pane->attribute(COM_NC)->pointer());
    Vector_3<double> *cs = reinterpret_cast<Vector_3<double>*>
      ( pane->attribute(c_attr->id())->pointer());
    int *cranks = reinterpret_cast<int*>
      ( pane->attribute(cntnranks_id)->pointer());

    std::set< Edge_ID> &eset_pn = _edges[i];
    for ( std::set< Edge_ID>::const_iterator eit=eset_pn.begin(),
            eend=eset_pn.end(); eit!=eend; ++eit) {
      Edge_ID eid = *eit;
      bool    is_border = eid.is_border();
      if ( is_border) eid = pm_it->get_opposite_real_edge( eid);

      // Loop through real elements of the current pane
      Element_node_enumerator ene( pane, eid.eid()); 
      int ne = ene.size_of_edges();
      
      int vindex=ene[eid.lid()]-1, windex=ene[(eid.lid()+1)%ne]-1;
      if ( is_border) std::swap( vindex, windex);

      if ( cranks[vindex]>0) {
        Vector_3<double> diff_wv = pnts[windex]-pnts[vindex];
        // OLD:: Vector_3<double> diff_wv = pnts[windex]-pnts[vindex] + ds[windex] - ds[vindex];

        // get average of edge centers
        cs[vindex][0] += 0.5*(es[2*vindex]*diff_wv);
      }
    }
  }

  _wm->reduce_on_shared_nodes( c_attr, Window_manifold_2::OP_SUM);

  // Loop through all vertices to computer center
  it = allpanes.begin();  pm_it=_wm->pm_begin();
  for ( int i=0, local_npanes = allpanes.size(); 
        i<local_npanes; ++i, ++it, ++pm_it) { 
    COM::Pane *pane = *it;

    Vector_3<double> *es = reinterpret_cast<Vector_3<double>*>
      ( pane->attribute(cntnvecs_id)->pointer());
    Vector_3<double> *cs = reinterpret_cast<Vector_3<double>*>
      ( pane->attribute( c_attr->id())->pointer());
    int *cranks = reinterpret_cast<int*>
      ( pane->attribute(cntnranks_id)->pointer());

    // Loop through all real nodes of the pane
    std::set< Edge_ID> &eset_pn = _edges[i];
    for ( std::set< Edge_ID>::const_iterator eit=eset_pn.begin(),
            eend=eset_pn.end(); eit!=eend; ++eit) {
      Edge_ID eid = *eit;
      bool    is_border = eid.is_border();
      if ( is_border) eid = pm_it->get_opposite_real_edge( eid);

      // Loop through real elements of the current pane
      Element_node_enumerator ene( pane, eid.eid()); 
      int ne = ene.size_of_edges();
      
      int vindex=ene[eid.lid()]-1, windex=ene[(eid.lid()+1)%ne]-1;
      if ( is_border) std::swap( vindex, windex);

      if ( cranks[vindex]>0)
        cs[vindex] = 0.5*cs[vindex][0]*es[2*vindex];
    }
  }

  // Loop through the panes and its real faces to compute safe factors.
  get_redist_safe_factor( c_attr, w_attr, 1);

  // Loop through all vertices to update displacement vector 
  it = allpanes.begin();  pm_it=_wm->pm_begin();
  for ( int i=0, local_npanes = allpanes.size(); 
        i<local_npanes; ++i, ++it, ++pm_it) { 
    COM::Pane *pane = *it;

    Vector_3<double> *cs = reinterpret_cast<Vector_3<double>*>
      ( pane->attribute( c_attr->id())->pointer());
    Vector_3<double> *ds = reinterpret_cast<Vector_3<double>*>
      ( pane->attribute( disps_id)->pointer());
    double   *ws = reinterpret_cast<double*>
      ( pane->attribute(w_attr->id())->pointer());

    // Loop through all real nodes of the pane
    std::set< Edge_ID> &eset_pn = _edges[i];
    for ( std::set< Edge_ID>::const_iterator eit=eset_pn.begin(),
            eend=eset_pn.end(); eit!=eend; ++eit) {
      Edge_ID eid = *eit;
      bool    is_border = eid.is_border();
      if ( is_border) eid = pm_it->get_opposite_real_edge( eid);

      // Loop through real elements of the current pane
      Element_node_enumerator ene( pane, eid.eid()); 
      int ne = ene.size_of_edges();
      
      int vindex=ene[eid.lid()]-1, windex=ene[(eid.lid()+1)%ne]-1;
      if ( is_border) std::swap( vindex, windex);

      double w=0.5;
      if ( ws[vindex]>0) w = std::min( w, 1./ws[vindex]);
      // displacement = w * average of two edge segments
      ds[vindex] += w*cs[vindex];
      // do PN triangle stuff here
    }
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

void redistribute_vertices_ridge ( )

protected

Redistribute ridge vertices within their tangent spaces.

void redistribute_vertices_ridge ( )

protected

Redistribute ridge vertices within their tangent spaces.

void redistribute_vertices_smooth ( )

protected

Redistribute smooth vertices within their tangent spaces.

Definition at line 682 of file smooth_medial.C.

References _wrk_window, cimg_library::acos(), COM_NC, Rocblas::copy_scalar(), get_redist_safe_factor(), half_pi, i, j, k, max(), min(), Element_node_enumerator::next(), nj, Window_manifold_2::OP_SUM, pi, s, Element_node_vectors_k_const< Value >::set(), Element_node_enumerator::size_of_edges(), and sqrt().

Referenced by smooth_surf_medial().

                               {

  // Use the following buffer spaces: _eigvalues for vector c 
  // (only lsast two components); _weights to store weight
  const std::string att1("disps"), att2("lambda"), att3("weights"),
    att4("cntnvecs"), att5("cntnranks"), att6("tangranks");
  int disps_id = _wrk_window->attribute(att1)->id();
  COM::Attribute *c_attr = _wrk_window->attribute(att2);
  COM::Attribute *w_attr = _wrk_window->attribute(att3);
  int cntnvecs_id = _wrk_window->attribute(att4)->id();
  int cntnranks_id = _wrk_window->attribute(att5)->id();
  int tangranks_id = _wrk_window->attribute("tangranks")->id();

  double zero = 0.;
  Rocblas::copy_scalar( &zero, c_attr);
  Rocblas::copy_scalar( &zero, w_attr);

  std::vector< COM::Pane*> allpanes;
  _wrk_window->panes(allpanes);
  // Loop through the panes and its real faces
  std::vector< COM::Pane*>::iterator it = allpanes.begin();
  for (int i=0, local_npanes = allpanes.size(); i<local_npanes; ++i, ++it) { 
    COM::Pane *pane = *it;
    
    Vector_3<double> *es = reinterpret_cast<Vector_3<double>*>
      ( pane->attribute(cntnvecs_id)->pointer());
    const Point_3<double> *pnts = reinterpret_cast<Point_3<double>*>
      (pane->attribute(COM_NC)->pointer());
    Vector_3<double> *cs = reinterpret_cast<Vector_3<double>*>
      ( pane->attribute(c_attr->id())->pointer());
    double   *ws = reinterpret_cast<double*>
      ( pane->attribute(w_attr->id())->pointer());
    Vector_3<double> *ds = reinterpret_cast<Vector_3<double>*>
      ( pane->attribute(disps_id)->pointer());
    int *tranks = reinterpret_cast<int*>
      ( pane->attribute(tangranks_id)->pointer());
    int *cranks = reinterpret_cast<int*>
      ( pane->attribute(cntnranks_id)->pointer());

    // Loop through real elements of the current pane
    Element_node_enumerator ene( pane, 1); 
    Element_node_vectors_k_const<Point_3<double> > ps; ps.set( pnts, ene, 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;

      // Loop through all vertices of current face.
      for ( int k=0; k<ne; ++k) {
        int windex = ene[(k+1==ne)?0:k+1]-1;

        // Only redistribute smooth vertices 
        if ( tranks[vindex]==2) {
          std::cout << "cranks[vindex]-1 = " << cranks[vindex]-1 << std::endl;
          Vector_3<double> tngs[2] = {pnts[uindex]-pnts[vindex], 
                              pnts[windex]-pnts[vindex]};
          Vector_3<double> disp_v = ds[vindex];
          Vector_3<double> diff_uv = ds[uindex]-disp_v+tngs[0];
          Vector_3<double> diff_wv = ds[windex]-disp_v+tngs[1];

          // Compute angle at vertex
          double theta = std::acos
            ( diff_uv*diff_wv/std::sqrt((diff_wv*diff_wv)*(diff_uv*diff_uv)));
          std::cout << "theta = " << theta << std::endl;

          if ( theta>half_pi) theta = std::max(0.,pi-theta);
          std::cout << "theta 2 = " << theta << std::endl;
  
          for ( int k=cranks[vindex]-1; k>=0; --k) {
            cs[vindex][k] += theta*(es[2*vindex+k]*(diff_uv+diff_wv));
            //      std::cout << "cs[" << vindex << "][" << k << "] += " 
            //      << theta << "*(" << es[2*vindex+k] << " * (" 
            //      << diff_uv << " + " << diff_wv << "))\n";
          }
          ws[vindex] += theta;
          //      std::cout << "ws[" << vindex << "] = " << ws[vindex]  << std::endl;
        }

        uindex=vindex; vindex=windex;
      }
    }
  }
  _wm->reduce_on_shared_nodes( c_attr, Window_manifold_2::OP_SUM);
  _wm->reduce_on_shared_nodes( w_attr, Window_manifold_2::OP_SUM);

  // Loop through all vertices to compute center of stars.
  it = allpanes.begin();
  for (int i=0, local_npanes = allpanes.size(); i<local_npanes; ++i, ++it) { 
    COM::Pane *pane = *it;

    Vector_3<double> *cs = reinterpret_cast<Vector_3<double>*>
      ( pane->attribute( c_attr->id())->pointer());
    double   *ws = reinterpret_cast<double*>
      ( pane->attribute( w_attr->id())->pointer());
    int *tranks = reinterpret_cast<int*>
      ( pane->attribute(tangranks_id)->pointer());
    Vector_3<double> *es = reinterpret_cast<Vector_3<double>*>
      ( pane->attribute(cntnvecs_id)->pointer());
    int *cranks = reinterpret_cast<int*>
      ( pane->attribute(cntnranks_id)->pointer());

    // Loop through all real nodes of the pane
    for ( int j=0, nj=pane->size_of_real_nodes(); j<nj; ++j) {
      if ( tranks[j]==2) {
        Vector_3<double> s(0,0,0);
        for ( int k=cranks[j]-1; k>=0; --k) {
          s += (cs[j][k]*es[2*j+k])/ws[j];
        }
        cs[j] = s;
      }
    }
  }

  // Loop through the panes and its real faces to compute safe factors.
  get_redist_safe_factor( c_attr, w_attr, 2);

  // Loop through all vertices to update displacement vector 
  it = allpanes.begin();
  for (int i=0, local_npanes = allpanes.size(); i<local_npanes; ++i, ++it) { 
    COM::Pane *pane = *it;

    Vector_3<double> *cs = reinterpret_cast<Vector_3<double>*>
      ( pane->attribute( c_attr->id())->pointer());
    double   *ws = reinterpret_cast<double*>
      ( pane->attribute( w_attr->id())->pointer());
    Vector_3<double> *ds = reinterpret_cast<Vector_3<double>*>
      ( pane->attribute( disps_id)->pointer());
    int *tranks = reinterpret_cast<int*>
      ( pane->attribute( tangranks_id)->pointer());

    // Loop through all real nodes of the pane
    for ( int j=0, nj=pane->size_of_real_nodes(); j<nj; ++j) {
      if ( tranks[j]==2) {
        double w=0.5;
        if ( ws[j]>0) w = std::min( w, 1./ws[j]);
        ds[j] += w*cs[j];
        std::cout << "ds[" << j << "] += " << w << " * " << cs[j] << std::endl;
      }
    }
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

void redistribute_vertices_smooth ( )

protected

Redistribute smooth vertices within their tangent spaces.

void redistribute_vertices_smooth ( )

protected

Redistribute smooth vertices within their tangent spaces.

static void reduce_sum_on_shared_nodes ( COM::Attribute *  att )

staticprotected

Perform a sum-reduction on the shared nodes for the given attribute.

void reduce_sum_on_shared_nodes ( COM::Attribute *  att )

staticprotected

Perform a sum-reduction on the shared nodes for the given attribute.

Definition at line 595 of file Rocmop.C.

References MPI_SUM.

                                                        {
  Pane_communicator pc(att->window(), att->window()->get_communicator());
  pc.init(att);
  pc.begin_update_shared_nodes();
  pc.reduce_on_shared_nodes(MPI_SUM);
  pc.end_update_shared_nodes();
}

static void reduce_sum_on_shared_nodes ( COM::Attribute *  att )

staticprotected

Perform a sum-reduction on the shared nodes for the given attribute.

void set_value	(	const char *	opt,
		const void *	val
	)

Set a Rocomp option.

Parameters

opt	the option to be set.
val	the new value for the option.

Option List:

• name (data type, default value)
• method (int, SMOOTH_NONE) smoothing method to use
  • SMOOTH_VOL_MESQ_WG (Volume smoothing via mesquite with ghost)
  • SMOOTH_VOL_MESQ_NG (Volume smoothing via mesquite no ghosts)
  • SMOOTH_SURF_MEDIAL (Surface smoothing via medial quadric)
  • SMOOTH_NONE (None, default).
• verbose (int, 0) verbose level, an integer greater than -1
• lazy (int, 0) check quality before smoothing? 0 or 1
• tol (float, 0.0) looping quality tolerance [0,180]
• niter (int, 1) max # iterations for smoothing
• ctol (float, .99) convergence subcycle tolerance [0,1]
• ncycle (int, 1) max # subcycles for convergence
• inverted (int, 0) tets are inverted? 0 or 1
• monotone (int, 0) apply a non-decreasing invariant?

void set_value	(	const char *	opt,
		const void *	val
	)

Set a Rocomp option.

Parameters

opt	the option to be set.
val	the new value for the option.

Option List:

• name (data type, default value)
• method (int, SMOOTH_NONE) smoothing method to use
  • SMOOTH_VOL_MESQ_WG (Volume smoothing via mesquite with ghost)
  • SMOOTH_VOL_MESQ_NG (Volume smoothing via mesquite no ghosts)
  • SMOOTH_SURF_MEDIAL (Surface smoothing via medial quadric)
  • SMOOTH_NONE (None, default).
• verbose (int, 0) verbose level, an integer greater than -1
• lazy (int, 0) check quality before smoothing? 0 or 1
• tol (float, 0.0) looping quality tolerance [0,180]
• niter (int, 1) max # iterations for smoothing
• ctol (float, .99) convergence subcycle tolerance [0,1]
• ncycle (int, 1) max # subcycles for convergence
• inverted (int, 0) tets are inverted? 0 or 1
• monotone (int, 0) apply a non-decreasing invariant?

Definition at line 529 of file Rocmop.C.

References COM_assertion_msg, and cimg_library::cimg::option().

Referenced by load().

{
  if(_verb) 
    std::cout << "Entering Rocmop::set_value" << std::endl;

  COM_assertion_msg( validate_object()==0, "Invalid object");

  COM_assertion_msg( opt && value,
                      "Rocmop::set_value does not except NULL parameters");
  std::string option;
  if(opt) option = opt;
  if ( option == "method") {
    COM_assertion_msg( *((int*)value) <= SMOOTH_NONE && *((int*)value)>=0 
                      ,"Illegal value for 'method' option");
    _method = *((int*)value);
  }
  else if ( option == "verbose"){ 
    _verb = *((int*)value); }
  else if ( option == "lazy"){ 
    _lazy = *((int*)value); }
  else if ( option == "tol"){ 
    COM_assertion_msg( *((float*)value) <= 180. && *((float*)value)>=0. 
                      ,"Illegal value for 'method' option");
    _tol = *((float*)value); }
  else if ( option == "niter"){ 
    _niter = *((int*)value); }
  else if ( option == "ctol"){ 
    COM_assertion_msg( *((float*)value) <= 1. && *((float*)value)>=0. 
                      ,"Illegal value for 'method' option");
    _ctol = *((float*)value); }
  else if ( option == "ncycle"){
    _ncycle = *((int*)value); }
  else if ( option == "inverted"){
    _invert = *((int*)value); 
  }
  else COM_assertion_msg( false, "Unknown option");

  if(_verb > 1) 
    std::cout << "Exiting Rocmop::set_value" << std::endl;
}

Here is the call graph for this function:

Here is the caller graph for this function:

void set_value	(	const char *	opt,
		const void *	val
	)

Set a Rocomp option.

Parameters

opt	the option to be set.
val	the new value for the option.

Option List:

• name (data type, default value)
• method (int, SMOOTH_NONE) smoothing method to use
  • SMOOTH_VOL_MESQ_WG (Volume smoothing via mesquite with ghost)
  • SMOOTH_VOL_MESQ_NG (Volume smoothing via mesquite no ghosts)
  • SMOOTH_SURF_MEDIAL (Surface smoothing via medial quadric)
  • SMOOTH_NONE (None, default).
• verbose (int, 0) verbose level, an integer greater than -1
• lazy (int, 0) check quality before smoothing? 0 or 1
• tol (float, 0.0) looping quality tolerance [0,180]
• niter (int, 1) max # iterations for smoothing
• ctol (float, .99) convergence subcycle tolerance [0,1]
• ncycle (int, 1) max # subcycles for convergence
• inverted (int, 0) tets are inverted? 0 or 1
• monotone (int, 0) apply a non-decreasing invariant?

void smooth	(	const COM::Attribute *	pmesh,
		COM::Attribute *	disp
	)

Smooth the mesh in a Rocmop managed buffer.

Creates a copy of the mesh, which is then optimized in place via perform_smoothing().

Parameters

pmesh	The pmesh to be smoothed.
disp	The attribute where displacements to new positions should be stored.

Definition at line 132 of file Rocmop.C.

References COM_assertion_msg, COM_MESH, COM_NC, COM_PMESH, MPI_MAX, and Rocblas::sub().

Referenced by load().

                                       {
  COM_assertion_msg( validate_object()==0, "Invalid object");
  if(_verb)
    cout << "Entering Rocmop::smooth\n";

  int pmesh_id = pmesh->id();
  COM_assertion_msg( pmesh &&
                     (pmesh_id==COM::COM_PMESH || pmesh_id==COM::COM_MESH),
                     "Input to Rocmop::smooth must be a pmesh.");
  _is_pmesh = (pmesh_id == COM_PMESH) ? true : false;

  if(_wrk_window) delete _wrk_window;
  _usr_window = pmesh->window();

  // Create a buffer window by inheriting(clone) from the user's mesh.
  if (_verb >3) 
    std::cout << "  Creating(clone) buffer window." << std::endl;
  std::string buf_name(_usr_window->name()+"-Rocmopbuf");
  _wrk_window = new COM::Window(buf_name, 
                                _usr_window->get_communicator());
  _wrk_window->inherit( const_cast<COM::Attribute*>(pmesh), "", 
                        COM::Pane::INHERIT_CLONE, true, NULL, 0);
  _wrk_window->init_done();

  std::vector<const Pane*> allpanes;
  _wrk_window->panes(allpanes);

  double pre_worst = 180.0, post_worst = 0.0;

  // Check initial quality if we are in lazy or monotone mode
  if(_monotone || _lazy){
    pre_worst = check_all_elem_quality(allpanes);
    // get the worst pre smoothing quality across all panes
    agree_double(pre_worst, MPI_MAX);
  }

  bool to_smooth = (!_lazy || (pre_worst > 180.*_tol));
    
  int degraded = 0;
  if (to_smooth){
    perform_smoothing(pre_worst);
    
    // Check final mesh quality if a tolerance is set, or we're in monotonic mode
    if(_monotone || (_tol!=0.0)){
      post_worst = check_all_elem_quality(allpanes);
      agree_double(post_worst, MPI_MAX);
    }
    
    // If in monotone mode, see if quality degraded on any pane
    if(_monotone && (pre_worst > post_worst)){
      cerr << "Warning, smoothing would have degraded mesh quality \n" 
           << "so the mesh has not been modified.\n";
      degraded = 1;
    }
    
    // Warn the user if could not reach a specified minimum mesh quality.
    if(_tol && (post_worst < _tol)){
      cerr << "Warning, post-smoothing mesh quality is below specified minimum " 
           << "quality \n";
    }
  }
  
  const COM::Attribute *old_nc = pmesh->window()->attribute( COM::COM_NC);
  const COM::Attribute *new_nc = (degraded || !to_smooth) ? 
    old_nc :_wrk_window->attribute( COM::COM_NC);
  Rocblas::sub (new_nc, old_nc, disp);

  if (_verb > 1) 
    std::cout << "Exiting rocmop::smooth" << std::endl;
}

Here is the call graph for this function:

Here is the caller graph for this function:

void smooth	(	const COM::Attribute *	pmesh,
		COM::Attribute *	disp
	)

Smooth the mesh in a Rocmop managed buffer.

Creates a copy of the mesh, which is then optimized in place via perform_smoothing().

Parameters

pmesh	The pmesh to be smoothed.
disp	The attribute where displacements to new positions should be stored.

void smooth	(	const COM::Attribute *	pmesh,
		COM::Attribute *	disp,
		double *	timestep = NULL
	)

Smooth the mesh in a Rocmop managed buffer.

Creates a copy of the mesh, which is then optimized in place via perform_smoothing().

Parameters

pmesh	The pmesh to be smoothed.
disp	The attribute where displacements to new positions should be stored.
timestep	If a timestep is supplied, then disp will contain velocity vectors instead of displacement vectors

Definition at line 700 of file Rocmop_2.C.

References COM_assertion_msg, COM_MESH, COM_PMESH, and MPI_MAX.

                                     {

  // The static var N will statically store the current value of 
  // N.  When N == 0, we will smooth and set N back to (_smoothfreq - 1).
  // What this amounts to is that we smooth only every _smoothfreqth call
  // to Rocmop.  _smoothfreq is specified as "N" in the config file.
  static int N = (_smoothfreq-1);
  
  if(N == 0){ // Go ahead and actually smooth if this is true
    N = (_smoothfreq-1); // Reset N 
    
    // We want to start profiling from here as it's smoothing we
    // want to time, not the NOOP calls.
#ifdef ROCPROF
    Profile_begin("Rocmop::smooth");
#endif
    
    COM_assertion_msg( validate_object()==0, "Invalid object");
    
    // Verify that we are working with a mesh or pmesh
    int pmesh_id = pmesh->id();
    COM_assertion_msg( pmesh &&
                       (pmesh_id==COM::COM_PMESH || pmesh_id==COM::COM_MESH),
                       "Input to Rocmop::smooth must be a mesh or pmesh attribute.");
    _is_pmesh = (pmesh_id == COM_PMESH) ? true : false;
    
    // If buffer window exists and was cloned from the current user
    // window, then just update the coordinates.
   
    if(_buf_window && (_usr_window == pmesh->window()))
    {
      // update nodel coordinates (real part) from _usr_window
      update_buf_real_nc();    
    }
    
    // Else, inherit(clone) a buffer window from the user's mesh.
    else
    {
      // Eliminate any old buffers
      if(_buf_window)
        delete _buf_window;
      
      _usr_window = pmesh->window();

      // First, check if pconn of _usr_window is complete
      int pconnErr = 0;
      pconnErr = check_input_pconn();
      agree_int(pconnErr, MPI_MAX);
      if (pconnErr == 1)
        COM_assertion_msg(pconnErr == 0, "*** Rocmop:ERROR: PCONN incomplete ***");

      // Create a buffer window by cloning from the user window
      std::string buf_name(_usr_window->name()+"-Rocmopbuf");
      _buf_window = new COM::Window(buf_name, 
                                    _usr_window->get_communicator());
      _buf_window->inherit( const_cast<COM::Attribute*>(_usr_window->attribute(COM::COM_MESH)), "", 
                            COM::Pane::INHERIT_CLONE, false, NULL, 0);
      _buf_window->init_done();

      // Add the pconn
      if(_buf_window->size_of_panes_global() > 1){
        Pane_ghost_connectivity pgc(_buf_window);
        pgc.build_pconn();
      }
     
      // Perform initialization specific to the smoother.
      smoother_specific_init();
    }
   
    // Max possible worst element quality is 180.0, so this
    // default should result in mesh smoothing if _lazy==0
    double mesh_qual = 190.0;
    
    // Check initial quality if we are in lazy mode
    if(_lazy){
      
      // Fill a vector with pointers to the local panes
      std::vector<const Pane*> allpanes;
      _buf_window->panes(allpanes);
      mesh_qual = check_all_elem_quality(allpanes);
    }
    
    // Mechanism to trigger smoothing on tallied displacements
    bool exceeded = true;
    if(_disp_thresh > 0.0)
      exceeded = check_displacements(disp);
    
    if(exceeded || ((mesh_qual > _tol) && _lazy)){     
      print_legible(0,"Smoothing...");
      perform_smoothing();    
      print_legible(0,"Smoothing complete.");

      // At this point NC in _buf_window has been smoothed.
      // Get displacement.
      get_usr_disp(pmesh, disp, timestep);
    }
    else
    {
      _usr_window = pmesh->window();
      zero_displacements(disp);
    }
    
#ifdef ROCPROF
    Profile_end("Rocmop::smooth");
#endif
  } // if (N says it's smoothing time
  else{
    if(_smoothfreq > 0)
      N--;
    _usr_window = pmesh->window();
    zero_displacements(disp);
  }
}

void smooth_boeing ( COM::Attribute *  att, int *  niter  )

protected

Definition at line 1471 of file Rocmop_2.C.

                                                       {
  /*
  std::vector<Pane*> allpanes;
  att->window()->panes(allpanes);

  COM_assertion_msg(allpanes.size() == 1,
                    "This function only supports winows with a single pane\n");

  MesqPane mp(allpanes[0],false);

  if(_wrapper == WRAPPER_SHAPE){
    
    Mesquite::MsqError err;
    ShapeImprovementWrapper mesh_quality_algorithm(err);
    MSQ_ERRRTN(err);
    for(int i=0; i< *niter; ++i){
      mesh_quality_algorithm.run_instructions(&mp,err);    
      MSQ_ERRRTN(err);
    }
  }
  else if(_wrapper == WRAPPER_BOEING){
    Mesquite::MsqError err;
    CGWrapper mesh_quality_algorithm(err);
    MSQ_ERRRTN(err);
    for(int i=0; i< *niter; ++i){
      mesh_quality_algorithm.run_instructions(&mp,err);    
      MSQ_ERRRTN(err);
    }
  }
  */
}

void smooth_in_place

(

COM::Attribute *

pmesh

)

Smooth a mesh in place..

The mesh is optimized in place via perform_smoothing().

Parameters

pmesh

The pmesh to be smoothed.

Definition at line 204 of file Rocmop.C.

References COM_assertion_msg, COM_MESH, COM_PMESH, and MPI_MIN.

Referenced by load().

                                               {
  COM_assertion_msg( validate_object()==0, "Invalid object");
  if (_verb)
    std::cout << "Entering rocmop::smooth_in_place" << std::endl;

  int pmesh_id = pmesh->id();

  COM_assertion_msg( pmesh && 
                     (pmesh_id==COM::COM_PMESH || pmesh_id==COM::COM_MESH) , 
                     "Input to Rocmop::smooth_in_place must be a mesh or pmesh");

  if(_wrk_window) 
    delete _wrk_window;
  _usr_window = pmesh->window();

  // FIXME
  // Create a buffer window by inheriting(use) from the user's mesh.
  // can't use the _usr_window directly because of problems with
  // Element_node_enumerator
  _wrk_window = new COM::Window(_usr_window->name()+"-Rocmopbuf", 
                                _usr_window->get_communicator());
  _wrk_window->inherit( const_cast<COM::Attribute*>(pmesh), "", 
                        false, true, NULL, 0);
  _wrk_window->init_done();

  std::vector<const Pane*> allpanes;
  _wrk_window->panes(allpanes);

  double pre_worst = 180.0, post_worst = 0.0;

  // Check initial quality if we are in lazy or monotone mode
  if(_monotone || _lazy) {
    pre_worst = check_all_elem_quality(allpanes);
    // get the worst pre smoothing quality on all panes
    agree_double(pre_worst, MPI_MIN);
  }

  bool to_smooth = (!_lazy || (pre_worst > 180.*_tol));

  if (to_smooth){
    perform_smoothing(pre_worst);
    
    // Check final mesh quality if a tolerance is set, or we're in monotonic mode
    if(_monotone || (_tol!=0.0)){
      post_worst = check_all_elem_quality(allpanes);
      agree_double(post_worst, MPI_MIN);
    }
    
    // If in monotone mode, see if quality degraded on any pane
    if(_monotone && (pre_worst > post_worst))
      cerr << "Warning, mesh quality degraded during smoothing in place." << endl;

    // If a quality tolerance is set, and we didn't meet it, then warn
    if(_tol && (post_worst < _tol)){
      cerr << "Warning, post-smoothing mesh quality is below specified minimum " 
           << "quality \n";
    }
  }

  if (_wrk_window){ delete _wrk_window;_wrk_window = NULL;}

  if (_verb > 1)
    std::cout << "Exiting rocmop::smooth_in_place" << std::endl;
}

Here is the caller graph for this function:

void smooth_in_place

(

COM::Attribute *

pmesh

)

Smooth a mesh in place..

The mesh is optimized in place via perform_smoothing().

Parameters

pmesh

The pmesh to be smoothed.

void smooth_mesquite ( std::vector< COM::Pane * > &  allpanes, int  ghost_level = 0  )

protected

Smooth the panes of the working window using MESQUITE.

This method creates MesqPane attributes for each of the local panes, which are then smoothed serially.

Parameters

ghost_level

whether to use ghost cells and nodes.

Definition at line 570 of file Rocmop.C.

References MeshSet::add_mesh(), j, MSQ_CHKERR, ShapeImprovementWrapper::run_instructions(), MesqPane::set_verb(), and total_npanes.

                                             {

  Mesquite::MsqError err;
  MesqPane *mp;
  ShapeImprovementWrapper mesh_quality_algorithm;

  int total_npanes = (int)allpanes.size();
  bool wg = (ghost_level == 0) ? false : true;
  for(int j=0; j < total_npanes; ++j){
    MeshSet* mesh_set1 = new MeshSet;
    mp = new MesqPane(allpanes[j], wg);
    if(_verb > 4) 
      mp->set_verb(_verb - 4);
    mesh_set1->add_mesh(mp, err); 
    MSQ_CHKERR(err);
    mesh_quality_algorithm.run_instructions(*mesh_set1, err); 
    MSQ_CHKERR(err);
    delete mesh_set1;
    if(mp)
      delete mp;
    mp = NULL;
  }
}

Here is the call graph for this function:

void smooth_mesquite ( std::vector< COM::Pane * > &  allpanes, int  ghost_level = 0  )

protected

Smooth the panes of the buffer window using MESQUITE.

This method creates MesqPane attributes for each of the local panes, which are then smoothed serially.

Parameters

ghost_level

whether to use ghost cells and nodes.

void smooth_mesquite ( std::vector< COM::Pane * > &  allpanes, int  ghost_level = 0  )

protected

Smooth the panes of the buffer window using MESQUITE.

This method creates MesqPane attributes for each of the local panes, which are then smoothed serially.

Parameters

ghost_level

whether to use ghost cells and nodes.

void smooth_surf_medial ( )

protected

Smooths a surface using the medial quadric.

Definition at line 230 of file smooth_medial.C.

References _rediter, _usr_window, _verb, _wrk_window, COM_assertion_msg, COM_NC, compute_medial_quadric(), Rocblas::copy(), evaluate_face_normals(), i, redistribute_vertices_ridge(), and redistribute_vertices_smooth().

                               {
  if(_verb > 1) 
    std::cout << "Entering Rocmop::smooth_medial" << std::endl;

  COM_assertion_msg(_wrk_window, "Unexpected NULL pointer encountered.");

  // Calculate the face normals
  Rocmop::evaluate_face_normals();

  // Compute the medial quadric
  Rocmop::compute_medial_quadric();

#if EXPORT_DEBUG_INFO
 COM::Attribute *vnormals = const_cast<COM::Attribute*>(_usr_window->attribute("vnormals"));
 COM::Attribute *vnrms = _wrk_window->attribute("vnormals");
 COM::Attribute *awnormals = const_cast<COM::Attribute*>(_usr_window->attribute("angleweighted"));
 COM::Attribute *awnrms = _wrk_window->attribute("awnormals");
 COM::Attribute *uwnormals = const_cast<COM::Attribute*>(_usr_window->attribute("unitweighted"));
 COM::Attribute *uwnrms = _wrk_window->attribute("uwnormals");
 COM::Attribute *tangranks = const_cast<COM::Attribute*>(_usr_window->attribute("tangranks"));
 COM::Attribute *tranks = _wrk_window->attribute("tangranks");
 COM::Attribute *evals = _wrk_window->attribute("lambda");
 COM::Attribute *eigenvalues = const_cast<COM::Attribute*>(_usr_window->attribute("eigenvalues"));
 COM::Attribute *ws = _wrk_window->attribute("weights");
 COM::Attribute *weights = const_cast<COM::Attribute*>(_usr_window->attribute("adefects"));

 if ( vnormals)  Rocblas::copy( vnrms, vnormals);
 if ( awnormals)  Rocblas::copy( awnrms, awnormals);
 if ( uwnormals)  Rocblas::copy( uwnrms, uwnormals);
 if ( tangranks) Rocblas::copy( tranks, tangranks);
 if ( eigenvalues) Rocblas::copy(evals, eigenvalues);
 if ( weights) Rocblas::copy(ws,weights);
#endif
 
 // Identify ridge edges
 //identify_ridge_edges();
 
 // Redistribute vertices within their tangent spaces
 // FIXME.... MODIFY THESE TO USE PN TRIANGLES
 for ( int i=0; i<_rediter; ++i) {
   redistribute_vertices_ridge();
   redistribute_vertices_smooth();
 }
 
 // Add the displacements onto the nodal coords
 // FIXME
 const std::string att1("disps");
 COM::Attribute *disps = _wrk_window->attribute(att1);
 COM::Attribute *nc = _wrk_window->attribute( COM::COM_NC);
 // Rocblas::add(disps,nc,nc);
 
 if(_verb > 2) 
   std::cout << "Exiting Rocmop::smooth_medial" << std::endl;
}

Here is the call graph for this function:

void smooth_surf_medial ( )

protected

Smooths a surface using the medial quadric.

void smooth_surf_medial ( )

protected

Smooths a surface using the medial quadric.

void smooth_vol_mesq_ng ( )

inlineprotected

Smooths a volume using Mesquite with only shared node information.

This method operates in two steps. First the shared nodes which are not on the physical boundary of the mesh are redistributed via element based laplacian smoothing. Then, the real nodes of the panes are smoothed individually with pane boundaries fixed.

Parameters

pre_quality

Mesh quality prior to smoothing.

Definition at line 276 of file Rocmop_1.h.

                           {
    ;
  }

void smooth_vol_mesq_ng ( double  pre_quality )

protected

Smooths a volume using Mesquite with only shared node information.

This method operates in two steps. First the shared nodes which are not on the physical boundary of the mesh are redistributed via element based laplacian smoothing. Then, the real nodes of the panes are smoothed individually with pane boundaries fixed.

Parameters

pre_quality

Mesh quality prior to smoothing.

Definition at line 64 of file smooth_mesquite_ng.C.

References COM_assertion_msg, COM_DOUBLE, COM_INT, COM_NC, COM_PCONN, COMMPI_Initialized(), Angle_Metric_3::compute(), Rocblas::copy(), Rocblas::copy_scalar(), dist(), Rocblas::div(), i, Geo_Metric_Base_3::initialize(), j, k, MPI_MAX, MPI_MIN, ni, nj, nk, Vector_3< Type >::norm(), nvc::norm(), Vector_3< Type >::normalize(), rank, Rocmap::reduce_maxabs_on_shared_nodes(), and total_npanes.

                                                     {
  if(_verb) 
    std::cout << "      Entering Rocmop::smooth_mesquite_ng" << std::endl;

  const std::string surf_attr("is_surface");
  COM::Attribute* w_is_surface = _wrk_window->attribute(surf_attr);
  COM_assertion_msg( w_is_surface, "Unexpected NULL pointer");

  // First Perform Element-based Laplacian smoothing on pane boundary volume nodes
  // implemented as follows:

  // * this part should go into smoother_specific_init
  // 0 Declare and resize variables and buffers
  // 1 Loop through panes
  //   a. Find the set of shared nodes.
  //   b. Remove surface nodes from above set, "leaving submerged boundary nodes"
  //   c. Loop through all elements
  //      if(element contains submersed boundary nodes)
  //          increment submersed_boundary nodes' adj. element count
  //          add element center to each submersed_boundary nodes' new position


  // 2 Sum reduction on nodal positions and adj. element counts
  // 3 Divide new nodal positions by element counts

  // 0 Declare and resize data structures.
  if(_verb > 2) std::cout << "  Declaring Variables\n";

  // Allocate buffer space for new nodal positions
  COM::Attribute *w_new_coords = 
    _wrk_window->new_attribute("new_coords", 'n', COM_DOUBLE, 3, "");
  _wrk_window->resize_array(w_new_coords, 0);

  // Allocate buffer space for adjacent element counts
  COM::Attribute *w_adj_elem_cnts =
    _wrk_window->new_attribute("adj_elem_cnts", 'n', COM_DOUBLE, 1, "");
  _wrk_window->resize_array(w_adj_elem_cnts, 0);
  _wrk_window->init_done();

  // Allocate buffer space for safe distance
  COM::Attribute *w_safe_dist =
    _wrk_window->new_attribute("safe_dist", 'n', COM_DOUBLE, 1, "");
  _wrk_window->resize_array(w_safe_dist, 0);
  _wrk_window->init_done();

  // Allocate buffer space for backup coords
  COM::Attribute *w_backup =
    _wrk_window->new_attribute("backup", 'n', COM_DOUBLE, 3, "");
  _wrk_window->resize_array(w_backup, 0);
  _wrk_window->init_done();

  // Allocate buffer space for reset flags
  COM::Attribute *w_reset =
    _wrk_window->new_attribute("reset", 'n', COM_INT, 1, "");
  _wrk_window->resize_array(w_reset, 0);
  _wrk_window->init_done();

  // Initialize buffer spaces to zeroes.
  double dzero = 0.;
  double dlarge =99999999;
  int izero = 0;
  Rocblas::copy_scalar( &dzero, w_new_coords);
  Rocblas::copy_scalar( &dzero, w_adj_elem_cnts);
  Rocblas::copy_scalar( &dlarge, w_safe_dist);
  Rocblas::copy_scalar( &izero, w_reset);

  //backup nodal coords.
  Rocblas::copy( _wrk_window->attribute(COM::COM_NC), w_backup);

  // Get worst qualities, pre smoothing.

  // Get a list of local panes
  std::vector<COM::Pane*> allpanes;
  _wrk_window->panes(allpanes);
  int total_npanes = (int)allpanes.size();

  // Allocate space for the set of submerged_boundary.
  std::vector<std::vector<bool> > is_sub_bnd; // Is submerged boundary node?
  std::vector<std::vector<bool> > is_pan_bnd; // Is pane boundary node?
  is_sub_bnd.resize(total_npanes);
  is_pan_bnd.resize(total_npanes);

  // Get attribute ids
  int new_coords_id = w_new_coords->id();
  int adj_elem_cnts_id = w_adj_elem_cnts->id();
  int is_surface_id = w_is_surface->id();
  int safe_dist_id = w_safe_dist->id();
  int reset_id = w_reset->id();
  int backup_id = w_backup->id();

  // get the pconn offset
  int pconn_offset = MAP::Pane_connectivity::pconn_offset();

  // 1 Loop through panes

  if(_verb > 2) std::cout << "  Step 1\n";
  for(int i=0; i < total_npanes; ++i){
    int size_of_real_nodes = allpanes[i]->size_of_real_nodes();
    int size_of_real_elems = allpanes[i]->size_of_real_elements();
    is_sub_bnd[i].resize(size_of_real_nodes,0);
    is_pan_bnd[i].resize(size_of_real_nodes,0);

    std::vector<bool> is_isolated; // is a node isolated?
    MAP::Pane_boundary pb (allpanes[i]);
    pb.determine_border_nodes(is_pan_bnd[i], is_isolated);

    //   a. Find the set of shared nodes.

    // get pane level pointers
    COM::Attribute *p_is_surface = allpanes[i]->attribute(is_surface_id);
    int *is_surface_ptr = (int*)p_is_surface->pointer();

    double * adj_elem_cnts_ptr = reinterpret_cast<double*>
      (allpanes[i]->attribute(adj_elem_cnts_id)->pointer());

    double * safe_dist_ptr = 
      (double*)(allpanes[i]->attribute(safe_dist_id)->pointer());

    const Vector_3<double> *pnts = reinterpret_cast<Vector_3<double>*>
      (allpanes[i]->attribute(COM_NC)->pointer());

    Vector_3<double> *new_coords_ptr = reinterpret_cast<Vector_3<double>*>
      (allpanes[i]->attribute(new_coords_id)->pointer());
    
    // Obtain the pane connectivity of the local pane
    const COM::Attribute *pconn = allpanes[i]->attribute(COM::COM_PCONN);
    const int *vs = (const int*)pconn->pointer() + pconn_offset;
    int vs_size = pconn->size_of_real_items() - pconn_offset;

    // Determine the number of communicating panes for shared nodes.
    int count=0;
    for (int j=0, nj=vs_size; j<nj; j+=vs[j+1]+2) {
      if (_wrk_window->owner_rank( vs[j]) >=0) ++count;
    }

    int index = 0;
    // Loop through communicating panes for shared nodes.
    for ( int j=0; j<count; ++j, index+=vs[index+1]+2) {
      // We skip the panes that are not in the current window 
      while ( _wrk_window->owner_rank(vs[index])<0) {
        index+=vs[index+1]+2;
        COM_assertion_msg( index<=vs_size, "Invalid communication map");
      } 
      // Mark node as shared
      for(int k=0; k<vs[index+1]; ++k){
        is_sub_bnd[i][vs[index+2+k]-1]=1;
      }
    }
    
    //   b. Remove surface nodes from above set, leaving submerged boundary nodes
    for (int j =0; j < size_of_real_nodes; ++j){
      if(is_surface_ptr[j])
        is_sub_bnd[i][j] = 0;
    }

    //   c. Loop through all elements
    //      if(element contains submersed boundary nodes)
    //          increment submersed_boundary nodes' adj. element count
    //          add element center to each submersed_boundary nodes' new position
    Element_node_enumerator ene(allpanes[i],1);
    for(int j =0; j < size_of_real_elems; ++j, ene.next()){
      // Collect element indices of submerged boundary nodes.
      int nn = ene.size_of_nodes(), flag = 0;
      for(int k =0; k < nn; ++k){
        if(is_sub_bnd[i][ene[k]-1]){
          flag = 1;
          adj_elem_cnts_ptr[ene[k]-1] += 1.0;
        }
      }
      // If the element contains submerged boundary nodes
      if(flag!=0){
        // Calculate the element's center
        Vector_3<double> elem_cntr_pos(0.0,0.0,0.0);
        for(int k =0; k < nn; ++k){
          elem_cntr_pos += pnts[ene[k]-1];
        }              
        elem_cntr_pos /= (double)nn;
        // Add the center position to the accumulating new positions
        // of the constituent submerged boundary nodes.
        // And update safe distance if needed   
        for(int k =0; k < nn; ++k){
          if(is_sub_bnd[i][ene[k]-1]){
            // initialize safe distance with distance to center
            double local_safe_dist = (pnts[ene[k]-1] - elem_cntr_pos).norm();
            int long_src = -1;
            double long_dist = -1.0;
            // replace safe distance with .5 shortest adj. edge length
            // if less than distance to center.
            int safe_source = -1;
            for(int ii=1; ii<nn; ++ii){
              double cur_edge_dist =
                .5*(pnts[ene[k]-1] - pnts[ene[(k+ii)%nn]-1]).norm();
              if(cur_edge_dist < local_safe_dist){
                safe_source = ii;
                local_safe_dist = cur_edge_dist;
              }
              if(cur_edge_dist > long_dist){
                long_dist = cur_edge_dist;
                long_src = ii;
              }
            }
            new_coords_ptr[ene[k]-1] += elem_cntr_pos;
#if 0
            if(safe_source == -1)
              new_coords_ptr[ene[k]-1] += elem_cntr_pos;
            else {
              Vector_3<double> direction = 
                (pnts[ene[(k+long_src)%nn]-1] - pnts[ene[k]-1]);
              direction.normalize();
              new_coords_ptr[ene[k]-1] += pnts[ene[k]-1];
              new_coords_ptr[ene[k]-1] += direction*local_safe_dist;          
              //new_coords_ptr[ene[k]-1] -= pnts[ene[(k+safe_source)%nn]-1];
            }
#endif
            if(local_safe_dist < safe_dist_ptr[ene[k]-1])
              safe_dist_ptr[ene[k]-1] = local_safe_dist;
          }
        }
      }
    }
  }

  // 2 Sum reduction on nodal positions and adj. element counts
  if(_verb > 1) std::cout << "  Step 2\n";
  reduce_sum_on_shared_nodes(w_new_coords);
  reduce_sum_on_shared_nodes(w_adj_elem_cnts);
  if(COMMPI_Initialized()){
    MAP::Pane_communicator pc(_wrk_window, _wrk_window->get_communicator());
    pc.init(w_safe_dist);
    pc.begin_update_shared_nodes();
    pc.reduce_on_shared_nodes(MPI_MIN);
    pc.end_update_shared_nodes();
  }

  // 3 Divide new nodal positions by element counts
  if(_verb > 1) std::cout << "  Step 3\n";
  Rocblas::div(w_new_coords, w_adj_elem_cnts, w_new_coords);

  //std::string msg("\n  I submerged quality = ");
  //print_mquality(msg, is_sub_bnd);
  //msg = "  I pane boundary quality = ";
  //print_mquality(msg, is_pan_bnd);

  std::vector<std::vector<bool> > elem_to_check;
  mark_elems_from_nodes(is_sub_bnd, elem_to_check);
  
  // get min,max angles for later use
  double max_angle = 0.0;
  double min_angle = 180.0;
  double angles[] = {0.0, 0.0};
  for(int i=0,ni = allpanes.size(); i<ni; ++i){
    for(int k =0,nk = elem_to_check[i].size(); k<nk; ++k){
      if(elem_to_check[i][k]){
        Element_node_enumerator ene(allpanes[i],k+1);
        Angle_Metric_3 am;
        am.initialize(ene);
        am.compute(angles);
        if(angles[1]>max_angle)
          max_angle = angles[1];
        if(angles[0]<min_angle)
          min_angle = angles[0];
      }
    }
  }
  int rank =0;
  double temp = min_angle;
  if(COMMPI_Initialized()){
    agree_double(max_angle,MPI_MAX);
    agree_double(min_angle,MPI_MIN);
    int ierr = MPI_Comm_rank( _wrk_window->get_communicator()
                              , &rank); 
    assert( ierr == 0);
  }  

  if(_verb > 1) std::cout << "  Step 4\n";
  // 4 Copy new nodal coords into mesh for submerged boundary nodes

  for(int i=0; i < total_npanes; ++i){

    // get pane level attributes
    Vector_3<double>* new_coords_ptr = 
      reinterpret_cast<Vector_3<double>* >(allpanes[i]->
                                           attribute(new_coords_id)->
                                           pointer());
    Vector_3<double>* coords_ptr = 
      reinterpret_cast<Vector_3<double>* >(allpanes[i]->
                                           attribute(COM::COM_NC)
                                           ->pointer());
    
    double* safe_dist_ptr = 
      (double*)(allpanes[i]->attribute(safe_dist_id)->pointer());

    for (int j = 0, nj = allpanes[i]->size_of_real_nodes(); j<nj; ++j){
      if(is_sub_bnd[i][j]){
        Vector_3<double> direction = 
          new_coords_ptr[j] - coords_ptr[j];
        double dist = direction.norm();
        if (dist > safe_dist_ptr[j]){
          double alpha = .9*safe_dist_ptr[j]/dist;
          coords_ptr[j] = 
            alpha*new_coords_ptr[j] +(1.0-alpha)*coords_ptr[j];
        }
        else{
          coords_ptr[j] = .5*new_coords_ptr[j] +.5*coords_ptr[j];
        }
      }
    }
  }
  //msg = ("\n  L submerged quality = ");
  //print_mquality(msg, is_sub_bnd);
  //msg = "  L pane boundary quality = ";
  //print_mquality(msg, is_pan_bnd);
  //msg = "  L overall quality = ";
  //print_quality(msg);

  // Decide which positions to reset
  for(int i=0; i < total_npanes; ++i){
    
    int* reset_ptr =
      (int*)(allpanes[i]->attribute(reset_id)->pointer());
    
    for(int k =0,nk = elem_to_check[i].size(); k<nk; ++k){
      if(elem_to_check[i][k]){
        double angles[] = {0.0, 0.0};
        
        Element_node_enumerator ene(allpanes[i],k+1);
        Angle_Metric_3 am;
        am.initialize(ene);
        am.compute(angles);
        if(angles[1]>max_angle ||
           angles[0]<min_angle){
          std::vector<int> nodes;
          ene.get_nodes(nodes);
          for(int j =0, nj=nodes.size(); j<nj; ++j){
            if (is_sub_bnd[i][nodes[j]-1])
              reset_ptr[nodes[j]-1] = 1;
          }
        }
      }
    }
  }

  // All nodes being moved are shared, so this should be valid
  Rocmap::reduce_maxabs_on_shared_nodes(w_reset);

  for(int i=0; i < total_npanes; ++i){
    
    // get pane level attributes
    Vector_3<double>* backup_ptr = 
      reinterpret_cast<Vector_3<double>* >(allpanes[i]->
                                           attribute(backup_id)->
                                           pointer());
    Vector_3<double>* coords_ptr = 
      reinterpret_cast<Vector_3<double>* >(allpanes[i]->
                                           attribute(COM::COM_NC)
                                           ->pointer());  
    
    int* reset_ptr =
      (int*)(allpanes[i]->attribute(reset_id)->pointer());
    
    for(int j=0, nj = allpanes[i]->size_of_real_nodes(); j<nj; ++j){
      if(reset_ptr[j])
        coords_ptr[j] = backup_ptr[j];
    }
  }

  //msg = ("\n  R submerged quality = ");
  //print_mquality(msg, is_sub_bnd);
  //msg = "  R pane boundary quality = ";
  //print_mquality(msg, is_pan_bnd);
  //msg = "  R overall quality = ";
  //print_quality(msg);

  // Smooth using mesquite.
  smooth_mesquite(allpanes,0);

  //msg = "\n  M submerged quality = ";
  //print_mquality(msg, is_sub_bnd);
  // msg = "  M pane boundary quality = ";
  //print_mquality(msg, is_pan_bnd);
  //msg = "  M overall quality = ";
  //print_quality(msg);
    
  //Deallocate buffer space
  _wrk_window->delete_attribute( w_adj_elem_cnts->name());
  _wrk_window->delete_attribute( w_new_coords->name());
  _wrk_window->init_done();

  if(_verb > 1) 
    std::cout << "      Exiting Rocmop::smooth_mesquite_ng" << std::endl;
}

Here is the call graph for this function:

void smooth_vol_mesq_ng ( )

inlineprotected

Smooths a volume using Mesquite with only shared node information.

This method operates in two steps. First the shared nodes which are not on the physical boundary of the mesh are redistributed via element based laplacian smoothing. Then, the real nodes of the panes are smoothed individually with pane boundaries fixed.

Parameters

pre_quality

Mesh quality prior to smoothing.

Definition at line 298 of file Rocmop_2.h.

                           {
    ;
  }

void smooth_vol_mesq_wg ( )

protected

Smooth a volume via Mesquite using ghost information.

This method operates in two steps. First the panes, including ghost nodes, are smoothed individually using Mesquite. Then, shared nodes are moved to the average of their positions across all panes. These steps two steps repeat until either ncycles sub-cycles are completed or the loss in quality between the two steps is less than ctol.

Parameters

pre_quality

Mesh quality prior to smoothing.

Definition at line 58 of file smooth_mesquite_1.C.

References COM_NC, COM_PCONN, Rocmap::reduce_average_on_shared_nodes(), and Rocmap::update_ghosts().

                               {
  
  print_legible(1,"      Entering Rocmop::smooth_vol_mesq_wg");

  print_legible(2,"        Updating ghost node positions.");
  Rocmap::update_ghosts(_buf_window->attribute(COM::COM_NC),
                        _buf_window->attribute(COM::COM_PCONN));

  std::vector<COM::Pane*> allpanes;
  _buf_window->panes(allpanes);
  
  // smooth the panes via MESQUITE with ghosts.
  smooth_mesquite(allpanes,1);
  
  print_legible(2,"        Updating shared and ghost node positions.");
  Rocmap::reduce_average_on_shared_nodes(_buf_window->attribute(COM::COM_NC));
  Rocmap::update_ghosts(_buf_window->attribute(COM::COM_NC),
                        _buf_window->attribute(COM::COM_PCONN));
  
  print_legible(1,"      Exiting Rocmop::smooth_vol_mesq_wg");
}

Here is the call graph for this function:

void smooth_vol_mesq_wg ( double  pre_quality )

protected

Smooth a volume via Mesquite using ghost information.

This method operates in two steps. First the panes, including ghost nodes, are smoothed individually using Mesquite. Then, shared nodes are moved to the average of their positions across all panes. These steps two steps repeat until either ncycles sub-cycles are completed or the loss in quality between the two steps is less than ctol.

Parameters

pre_quality

Mesh quality prior to smoothing.

Definition at line 59 of file smooth_mesquite.C.

References COM_NC, COM_PCONN, i, MPI_MAX, Rocmap::reduce_average_on_shared_nodes(), total_npanes, and Rocmap::update_ghosts().

                                                     {
   if(_verb) 
    std::cout << "      Entering Rocmop::smooth_vol_mesq_wg" << std::endl;

  std::vector<COM::Pane*> allpanes;
  _wrk_window->panes(allpanes);
  int total_npanes = (int)allpanes.size();
  
  // Get the worst dihedral angle of all of the elements
  double pre_worst = initial_quality, post_worst = 180.0;
  double pre_worst_all = 180.0;
  double post_worst_all = 0.0;
  int to_cycle = 1;
  
  //std::string msg("\n  Initial shared quality = ");
  //print_mquality(msg, _is_shared_node);
  
  for(int i=0; (i<_ncycle && to_cycle); ++i){
    if(_verb > 2)
      std::cout << "        Subcycle " << i << std::endl;

    // smooth the panes via MESQUITE with ghosts.
    smooth_mesquite(allpanes,1);

    //msg = "\n  Post-mesquite quality = ";
    //print_mquality(msg, _is_shared_node);

    //msg = "\n  Post-mesquite all quality = ";
    //print_quality(msg);

    //get the worst quality of all the shared elements
    pre_worst = check_marked_elem_quality(_is_shared_elem,allpanes);
    
    if(_verb > 2)
      std::cout << "        Updating shared and ghost node positions.\n";
    Rocmap::reduce_average_on_shared_nodes(_wrk_window->attribute(COM::COM_NC));
    Rocmap::update_ghosts(_wrk_window->attribute(COM::COM_NC),
                          _wrk_window->attribute(COM::COM_PCONN));

    //msg = "\n  Post-communicated quality = ";
    //print_mquality(msg, _is_shared_node);

    //msg = "\n  Post-communicated all quality = ";
    //print_quality(msg);

    if(_ctol != 0.0){
      post_worst = check_marked_elem_quality(_is_shared_elem,allpanes);
      if(pre_worst/post_worst > _ctol)
        to_cycle = 0;
    }
    agree_int(to_cycle, MPI_MAX);
  }
  
  if(_verb > 1) 
    std::cout << "      Exiting Rocmop::smooth_vol_mesq_wg" << std::endl;
}

Here is the call graph for this function:

void smooth_vol_mesq_wg ( )

protected

Smooth a volume via Mesquite using ghost information.

This method operates in two steps. First the panes, including ghost nodes, are smoothed individually using Mesquite. Then, shared nodes are moved to the average of their positions across all panes. These steps two steps repeat until either ncycles sub-cycles are completed or the loss in quality between the two steps is less than ctol.

Parameters

pre_quality

Mesh quality prior to smoothing.

void smoother_specific_init ( )

protected

Perform smoother specific initialization.

Perform smoother specific initializing, for example initializing the Window_manifold_2 for a surface mesh, adding smoother specific attributes to the buffer window, etc.

Definition at line 364 of file Rocmop.C.

References COM_assertion_msg, COM_compatible_types(), COM_DOUBLE, COM_INT, COM_MESH, COM_PMESH, and Rocsurf::initialize().

                                   {
  if(_verb) 
    std::cout << "    Entering Rocmop::smoother_specific_init" << std::endl;

  // get rid of any old data
  if(_wm){ delete _wm; _wm = NULL; }

  // perform initialization common to surface smoothers
  if(_method==SMOOTH_SURF_MEDIAL){
    // Initialize the Window_manifold
    if(_wm == NULL)
      Rocsurf::initialize(_wrk_window->attribute(_is_pmesh ? COM_PMESH : COM_MESH));    

  }
  
  // perform smoother specific initialization 
  switch (_method){
    
#ifdef MESQUITE
  case SMOOTH_VOL_MESQ_WG: {
    // Obtain a list of elements containing shared nodes for each pane.
    // Don't bother if _ctol == 0 or _ncycle ==1
    //    if(!(_ctol==0.) && !(_ncycle==1))
    determine_shared_border();
    if(_invert)
      invert_tets();
    break;
  }
  case SMOOTH_VOL_MESQ_NG: {

    // Check to see if the physical surface boundary exists
    const std::string surf_attr("is_surface");
    const COM::Attribute *w_is_surface = _usr_window->attribute(surf_attr);
    if(w_is_surface){
      
      COM_assertion_msg( COM_compatible_types( w_is_surface->data_type(), COM_INT),
                         "Surface-list must have integer type");
      COM_assertion_msg( w_is_surface->is_nodal() == 1,
                         "Surface-list must be nodal");
      COM_assertion_msg( w_is_surface->size_of_components() == 1,
                         "Surface-list must have a single component");
      COM_assertion_msg( w_is_surface->initialized() == 1,
                         "Surface-list must be initialized");      

      // Clone the attribute
      COM::Attribute * new_attr = 
        _wrk_window->inherit( const_cast<COM::Attribute *>(w_is_surface), 
                              surf_attr, COM::Pane::INHERIT_CLONE, true, NULL, 0);
    }
    // else, detect the physical boundary ourselves
    else{
      COM::Attribute* w_surf_attr =
        _wrk_window->new_attribute( "is_surface", 'n', COM_INT, 1, "");
      _wrk_window->resize_array( w_surf_attr, 0);
      
      determine_physical_border(w_surf_attr);
    }
    _wrk_window->init_done();
    
    if(_invert)
      invert_tets();
    break;
  }
#else
  case SMOOTH_VOL_MESQ_WG:
  case SMOOTH_VOL_MESQ_NG:
    COM_assertion_msg(0, "Not compiled with MESQUITE");
    break;
#endif

  case SMOOTH_SURF_MEDIAL: {

    // Extend buffer window
    COM::Attribute* w_disps =
      _wrk_window->new_attribute( "disps", 'n', COM_DOUBLE, 3, "");
    _wrk_window->resize_array( w_disps, 0);

    COM::Attribute* w_facenormals = 
      _wrk_window->new_attribute( "facenormals", 'e', COM_DOUBLE, 3, "");
    _wrk_window->resize_array( w_facenormals, 0);
    
    COM::Attribute* w_facecenters = 
      _wrk_window->new_attribute( "facecenters", 'e', COM_DOUBLE, 3, "");
    _wrk_window->resize_array( w_facecenters, 0);

    COM::Attribute* w_eigvalues = 
    _wrk_window->new_attribute( "lambda", 'n', COM_DOUBLE, 3, "");
    _wrk_window->resize_array( w_eigvalues, 0);

    COM::Attribute* w_vnormals = 
    _wrk_window->new_attribute( "vnormals", 'n', COM_DOUBLE, 3, "");
    _wrk_window->resize_array( w_vnormals, 0);

    COM::Attribute* w_awnormals = 
    _wrk_window->new_attribute( "awnormals", 'n', COM_DOUBLE, 3, "");
    _wrk_window->resize_array( w_awnormals, 0);

    COM::Attribute* w_uwnormals = 
    _wrk_window->new_attribute( "uwnormals", 'n', COM_DOUBLE, 3, "");
    _wrk_window->resize_array( w_uwnormals, 0);

    COM::Attribute* w_eigvecs = 
    _wrk_window->new_attribute( "eigvecs", 'n', COM_DOUBLE, 9, "");
    _wrk_window->resize_array( w_eigvecs, 0);

    COM::Attribute* w_tangranks = 
    _wrk_window->new_attribute( "tangranks", 'n', COM_INT, 1, "");
    _wrk_window->resize_array( w_tangranks, 0);

    COM::Attribute* w_cntnranks = 
    _wrk_window->new_attribute( "cntnranks", 'n', COM_INT, 1, "");
    _wrk_window->resize_array( w_cntnranks, 0);
    
    COM::Attribute* w_cntnvecs =
    _wrk_window->new_attribute( "cntnvecs", 'n', COM_DOUBLE, 6, "");
    _wrk_window->resize_array( w_cntnvecs, 0);

    COM::Attribute* w_scales = 
    _wrk_window->new_attribute( "scales", 'n', COM_DOUBLE, 1, "");
    _wrk_window->resize_array( w_scales, 0);
    
    COM::Attribute* w_weights = 
    _wrk_window->new_attribute( "weights", 'n', COM_DOUBLE, 1, "");
    _wrk_window->resize_array( w_weights, 0);

    COM::Attribute* w_weights2 = 
    _wrk_window->new_attribute( "weights2", 'n', COM_DOUBLE, 1, "");
    _wrk_window->resize_array( w_weights2, 0);

    COM::Attribute* w_barycrds = 
    _wrk_window->new_attribute( "barycrds", 'n', COM_DOUBLE, 2, "");
    _wrk_window->resize_array( w_barycrds, 0);

    COM::Attribute* w_PNelemids = 
    _wrk_window->new_attribute( "PNelemids", 'n', COM_INT, 1, "");
    _wrk_window->resize_array( w_PNelemids, 0);

    // Extend the buffer window to hold local contributions to new placement
    // and the number of contributing faces for Laplacian smoothing.
    COM::Attribute * w_pnt_contrib = 
      _wrk_window->new_attribute("pnt_contrib", 'n', COM_DOUBLE, 3, "");
    _wrk_window->resize_array(w_pnt_contrib, 0);

    COM::Attribute * w_disp_count =
      _wrk_window->new_attribute("disp_count", 'n', COM_DOUBLE, 1, "");
    _wrk_window->resize_array(w_disp_count, 0);

    _wrk_window->init_done();

    break;
  }
    //  case SMOOTH_LAPLACE : {
    //break;
    //}
  default:
    COM_assertion_msg(0, "Can't initialize for invalid smoother.");
    break;
  }

  COM_assertion_msg(_wrk_window, "Unexpected NULL pointer encountered.");

  if(_verb > 1) 
    std::cout << "    Exiting Rocmop::smoother_specific_init" << std::endl;
}

Here is the call graph for this function:

void smoother_specific_init ( )

protected

Perform smoother specific initialization.

Perform smoother specific initializing, for example initializing the Window_manifold_2 for a surface mesh, adding smoother specific attributes to the buffer window, etc.

void smoother_specific_init ( )

protected

Perform smoother specific initialization.

Perform smoother specific initializing, for example initializing the Window_manifold_2 for a surface mesh, adding smoother specific attributes to the buffer window, etc.

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 378 of file Rocmop.h.

References denom.

Referenced by get_redist_safe_factor(), and 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 391 of file Rocmop.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 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 403 of file Rocmop_1.h.

References denom.

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

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

inlinestaticprotected

Definition at line 407 of file Rocmop.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:

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 416 of file Rocmop_1.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 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 430 of file Rocmop_2.h.

References denom.

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

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

inlinestaticprotected

Definition at line 432 of file Rocmop_1.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:

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 443 of file Rocmop_2.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< double >  A[3], const Vector_3< double > &  q, Vector_3< double > &  x  )

inlinestaticprotected

Definition at line 459 of file Rocmop_2.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:

void unload ( const std::string &  mname )

static

Unloads Rocmop from Roccom.

Definition at line 120 of file Rocmop.C.

References COM_delete_window(), and COM_get_object().

Referenced by COM_F_FUNC2(), and Rocmop_unload_module().

                                           {

  Rocmop *mop;
  std::string glb=mname+".global";

  COM_get_object( glb.c_str(), 0, &mop);
  delete mop;

  COM_delete_window( mname.c_str());

}

Here is the call graph for this function:

Here is the caller graph for this function:

static void unload ( const std::string &  mname )

static

Unloads Rocmop from Roccom.

static void unload ( const std::string &  mname )

static

Unloads Rocmop from Roccom.

void update_buf_real_nc ( )

protected

Update real nodal coordinates of _buf_window from _usr_window.

Use real nodal coordinates of _usr_window to update real nodal coordinates of _buf_window

Definition at line 824 of file Rocmop_2.C.

References COM_assertion_msg, COM_NC, i, j, ni, and nj.

{
  std::vector<const Pane*> allbufpanes;
  _buf_window->panes(allbufpanes);   
  std::vector<const Pane*> allusrpanes;
  _usr_window->panes(allusrpanes);

  COM_assertion_msg(allbufpanes.size() == allusrpanes.size(),
                    "Different number of panes on buffer and user windows.");

  for(int i = 0, ni = allbufpanes.size(); i < ni; ++i)
  {
    // Obtain pane level attributes
    const COM::Attribute *usr_nc = allusrpanes[i]->attribute( COM::COM_NC);
    const COM::Attribute *buf_nc = allbufpanes[i]->attribute( COM::COM_NC);
    
    COM_assertion_msg(usr_nc->stride() != 0,
                      "Rocmop can not operate on meshes with stride == 0");

    COM_assertion_msg(usr_nc->size_of_real_items() == buf_nc->size_of_real_items(),
                      "Number of real items differs between buffer and user windows");

    double* usr_nc_ptr = (double*)usr_nc->pointer();
    double* buf_nc_ptr = (double*)buf_nc->pointer();

    // get stride and distance_to_next_component
    int usr_nc_stride = usr_nc->stride();
    int buf_nc_stride = buf_nc->stride();
    int usr_nc_next_comp;
    int buf_nc_next_comp;
          
    if (usr_nc_stride == 1)
      usr_nc_next_comp = usr_nc->size_of_items();
    else
      usr_nc_next_comp = 1;
          
    if (buf_nc_stride == 1)
      buf_nc_next_comp = buf_nc->size_of_items();
    else
      buf_nc_next_comp = 1;
        
    // update real node coordinates
    for(int j=0, nj = usr_nc->size_of_real_items(); j<nj; ++j){
      *buf_nc_ptr                       = *usr_nc_ptr;
      *(buf_nc_ptr+ 1*buf_nc_next_comp) = *(usr_nc_ptr + 1*usr_nc_next_comp);
      *(buf_nc_ptr+ 2*buf_nc_next_comp) = *(usr_nc_ptr + 2*usr_nc_next_comp);
      
      buf_nc_ptr += buf_nc_stride;
      usr_nc_ptr += usr_nc_stride;
    }
  }
}

int validate_object ( ) const

inlineprotected

Check that the object is valid.

Definition at line 205 of file Rocmop_1.h.

References _cookie, and MOP_COOKIE.

                              {
    if ( _cookie != MOP_COOKIE) return -1;
    else return COM_Object::validate_object();
  }

int validate_object ( ) const

inlineprotected

Check that the object is valid.

Definition at line 206 of file Rocmop.h.

References _cookie, and MOP_COOKIE.

                              {
    if ( _cookie != MOP_COOKIE) return -1;
    else return COM_Object::validate_object();
  }

int validate_object ( ) const

inlineprotected

Check that the object is valid.

Definition at line 227 of file Rocmop_2.h.

References _cookie, and MOP_COOKIE.

                              {
    if ( _cookie != MOP_COOKIE) return -1;
    else return COM_Object::validate_object();
  }

void zero_displacements ( COM::Attribute *  disp )

protected

Definition at line 426 of file Rocmop_1.C.

References i, j, ni, and nj.

{
#ifdef ROCPROF
  Profile_begin("Rocmop::zero_disp");
#endif
  
  int disp_id = w_disp->id();
  std::vector<Pane*> allpanes;
  const_cast<COM::Window*>(_usr_window)->panes(allpanes);
    
  for(int i=0,ni = allpanes.size(); i<ni; ++i){
    Vector_3<double> *ptr = NULL;
    COM::Attribute* ptr_att = allpanes[i]->attribute(disp_id);
    void* void_ptr = ptr_att->pointer();
    ptr = reinterpret_cast<Vector_3<double>*>(void_ptr);
    for(int j=0,nj = allpanes[i]->size_of_real_nodes(); j<nj; ++j){     
      // Zero to maintain original physics module behavior
      ptr[j][0] = 0.0;
      ptr[j][1] = 0.0;
      ptr[j][2] = 0.0;
    }
  }
#ifdef ROCPROF
  Profile_end("Rocmop::zero_disp");
#endif
}

void zero_displacements ( COM::Attribute *  disp )

protected

Member Data Documentation

COM::Window * _buf_window

protected

The buffer window.

Definition at line 464 of file Rocmop_1.h.

COM::Attribute * _cnstr_dirs

protected

Stores directions of nodal contraints.

Definition at line 425 of file Rocmop.h.

Referenced by compute_medial_quadric().

const COM::Attribute * _cnstr_types

protected

Stores types of nodal constraints.

Definition at line 424 of file Rocmop.h.

Referenced by compute_medial_quadric().

int _cookie

protected

For Roccom.

Definition at line 487 of file Rocmop.h.

Referenced by validate_object().

float _ctol

protected

Subcycling tolerance.

Definition at line 471 of file Rocmop.h.

std::vector< MAP::Pane_dual_connectivity * > _dcs

protected

Definition at line 466 of file Rocmop_1.h.

double _dir_thres

protected

Another threshold.

Definition at line 421 of file Rocmop.h.

Referenced by eigen_analyze_vertex().

float _disp_thresh

protected

originally a static in check_displacements

Smooth every _smoothfreq'th call.

Definition at line 520 of file Rocmop_1.h.

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

protected

ridge edges

Definition at line 426 of file Rocmop.h.

Referenced by identify_ridge_edges(), and redistribute_vertices_ridge().

double _eig_thres

protected

Eigenvalue thresholds.

Definition at line 420 of file Rocmop.h.

Referenced by eigen_analyze_vertex().

int _invert

protected

If true (default false), then invert the mesh.

Definition at line 489 of file Rocmop.h.

int _invert_hexes

protected

If true (default false), then invert hexes.

Definition at line 543 of file Rocmop_2.h.

int _invert_tets

protected

If true (default false), then invert tets.

Definition at line 541 of file Rocmop_2.h.

std::vector< std::vector< bool > > _is_pane_bnd_elem

protected

Does the element contain nodes on the pane boundary?

Definition at line 436 of file Rocmop.h.

std::vector< std::vector< bool > > _is_pane_bnd_node

protected

Is the node on the pane boundary?

Definition at line 435 of file Rocmop.h.

std::vector< std::vector< bool > > _is_phys_bnd_elem

protected

Does the element contain nodes on the phys. boundary?

Definition at line 434 of file Rocmop.h.

std::vector< std::vector< bool > > _is_phys_bnd_node

protected

Is the node on the physical boundary?

Definition at line 433 of file Rocmop.h.

bool _is_pmesh

protected

pmesh or mesh ?

Definition at line 441 of file Rocmop.h.

std::vector< std::vector< bool > > _is_shared_elem

protected

Does the element contain shared nodes?

Definition at line 432 of file Rocmop.h.

std::vector< std::vector< bool > > _is_shared_node

protected

Is the node shared?

Definition at line 431 of file Rocmop.h.

int _lazy

protected

Check quality before smoothing?

Definition at line 462 of file Rocmop.h.

float _maxdisp

protected

Maximum displacement allowed.

Definition at line 512 of file Rocmop_1.h.

int _method

protected

Choice of smoothing method.

Definition at line 447 of file Rocmop.h.

int _monotone

protected

Impose non-decreasing quality?

Definition at line 481 of file Rocmop.h.

int _ncycle

protected

Max number of subcycles for convergence.

Definition at line 479 of file Rocmop.h.

int _niter

protected

Maximum number of iterations for smoothing.

Definition at line 469 of file Rocmop.h.

int _rediter

protected

No.iterations for vertex redistribution.

Definition at line 423 of file Rocmop.h.

Referenced by smooth_surf_medial().

bool _reorthog

protected

Reorthogonalize?

Definition at line 418 of file Rocmop.h.

Referenced by eigen_analyze_vertex().

double _saliency_crn

protected

Definition at line 422 of file Rocmop.h.

Referenced by eigen_analyze_vertex().

int _smoothfreq

protected

Definition at line 514 of file Rocmop_1.h.

float _tol

protected

Smoother iterating tolerance.

Definition at line 464 of file Rocmop.h.

const COM::Window * _usr_window

protected

The user's window.

Definition at line 438 of file Rocmop.h.

Referenced by agree_double(), agree_int(), and smooth_surf_medial().

int _verb

protected

Verbose level.

Definition at line 453 of file Rocmop.h.

Referenced by compute_medial_quadric(), evaluate_face_normals(), and smooth_surf_medial().

char _wght_scheme

protected

Weighting scheme.

Definition at line 419 of file Rocmop.h.

Referenced by compute_medial_quadric().

int _wrapper

protected

Choice of Mesquite Smoothing Wrappers.

Definition at line 507 of file Rocmop_2.h.

COM::Window* _wrk_window

protected

The working window.

Definition at line 439 of file Rocmop.h.

Referenced by compute_medial_quadric(), evaluate_face_normals(), get_redist_safe_factor(), identify_ridge_edges(), redistribute_vertices_ridge(), redistribute_vertices_smooth(), and smooth_surf_medial().

double disp_tally

protected

originally a static in smooth()

Definition at line 518 of file Rocmop_1.h.

int N

protected

Smooth every _smoothfreq'th call.

Definition at line 516 of file Rocmop_1.h.

---

The documentation for this class was generated from the following files:

• Rocmop.h
• Rocmop_1.h
• Rocmop_2.h
• Rocmop.C
• smooth_medial.C
• smooth_mesquite.C
• smooth_mesquite_ng.C
• Rocmop_1.C
• Rocmop_2.C
• smooth_medial_1.C
• smooth_medial_2.C
• smooth_mesquite_1.C
• smooth_mesquite_2.C
• smooth_mesquite_ng_1.C
• smooth_mesquite_ng_2.C