Rocmop_2.h

Go to the documentation of this file.

/* *******************************************************************
 * Rocstar Simulation Suite                                          *
 * Copyright@2015, Illinois Rocstar LLC. All rights reserved.        *
 *                                                                   *
 * Illinois Rocstar LLC                                              *
 * Champaign, IL                                                     *
 * www.illinoisrocstar.com                                           *
 * sales@illinoisrocstar.com                                         *
 *                                                                   *
 * License: See LICENSE file in top level of distribution package or *
 * http://opensource.org/licenses/NCSA                               *
 *********************************************************************/
/* *******************************************************************
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,   *
 * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES   *
 * OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND          *
 * NONINFRINGEMENT.  IN NO EVENT SHALL THE CONTRIBUTORS OR           *
 * COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER       *
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,   *
 * ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE    *
 * USE OR OTHER DEALINGS WITH THE SOFTWARE.                          *
 *********************************************************************/
// $Id: Rocmop_2.h,v 1.4 2008/12/06 08:45:24 mtcampbe Exp $

#ifndef __ROCMOP_H_
#define __ROCMOP_H_

#include "mopbasic.h"
#include "roccom_devel.h"
#include "Pane_communicator.h"
#include "commpi.h"
#include "Manifold_2.h"
#include "Rocsurf.h"
#include "Dual_connectivity.h"

MOP_BEGIN_NAMESPACE

class Rocmop : public SURF::Rocsurf {
public:
  
  enum { SMOOTH_VOL_MESQ_WG = 0, SMOOTH_VOL_MESQ_NG, 
         SMOOTH_SURF_MEDIAL, SMOOTH_NONE};

  enum { WRAPPER_SHAPE = 0, WRAPPER_BOEING, WRAPPER_MAX};

public:

  Rocmop() : _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();
  //\}

public:

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

  static void unload( const std::string &mname);
  //\}

  void read_config_file(const std::string &);
public:


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

protected:


  void smoother_specific_init();


  void perform_smoothing();


  void perform_iterative_smoothing();

  void perform_noniterative_smoothing();


  void update_buf_real_nc();


  int check_input_pconn();


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

  //\} 

public:


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

  void determine_physical_border(COM::Attribute* w_is_surface);

  void determine_physical_border();

  void determine_shared_border();

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

protected:
  virtual SURF::Window_manifold_2 *manifold()
  { return _wm;}
  
  int validate_object() const {
    if ( _cookie != MOP_COOKIE) return -1;
    else return COM_Object::validate_object();
  }
  
  void invert_elements(int conn_type);
  
protected:
  static void reduce_sum_on_shared_nodes(COM::Attribute *att);


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


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

protected:


  void smooth_vol_mesq_wg();


  void smooth_vol_mesq_ng(){
    ;
  }

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


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

protected:


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


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


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

  void print_legible(int verb, const char* msg);

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

  void print_extremal_dihedrals(COM::Window * window);

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

  void print_mquality(std::string &s,std::vector<std::vector<bool> > &to_check);
  //\}

  void perturb_stationary();

  // Check and tally displacements for triggering smoothing
  bool check_displacements(COM::Attribute *disp);

  // Zero displacement array
  void zero_displacements(COM::Attribute *disp);

public:

  // Add aspect ratio measures to the users mesh.
  void add_aspect_ratios(COM::Attribute *usr_att, 
                         COM::Attribute *buf_att =NULL);

 protected:

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

   void evaluate_face_normals();

   void identify_ridge_edges();

  void redistribute_vertices_ridge();

  void redistribute_vertices_smooth();
  
  void compute_medial_quadric();
  
  // Solve the minimization problem (x^T)Ax+2b^Tx with eigen-decomposition
  // at a vertex, return the codimension of the vertex, and save mean
  // normal in x_nz.
  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);

  // Compute eigenvalues and eigenvectors of a squared matrix A of order 3. 
  // At output, the eigenvalues are saved in lambdas in decending order,
  // and A is replaced by the orthonormal eigenvectors.
  static void compute_eigenvectors( Vector_3<double> A[3], 
                                    Vector_3<double> &lambdas);
  
  // Compute eigenvalues and eigenvectors of a squared matrix A of order 2.
  // At output, the eigenvalues are saved in lambdas in decending order,
  // and A is replaced by the orthonormal eigenvectors.
  static void compute_eigenvectors( Vector_2<double> A[2], 
                                     Vector_2<double> &lambdas);

  // Helpers for linear solvers for 2x2 and 3x3 equations.
  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)
  {
    FT denom = a1*b2-b1*a2;
    
    x = - (b1*c2-c1*b2)/denom;
    
    y = (a1*c2-c1*a2)/denom;
  }

  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)
  {
    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;
  }

  // Wrapper for solving 3x3 equations
  static void solve ( const Vector_3<double> A[3],
                      const Vector_3<double> &q,
                      Vector_3<double> &x) {
    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]);
  }
  //\}
  
 protected: // Variables specific to medial quadric smoothing

  bool _reorthog;                           
  char _wght_scheme;                        
  double _eig_thres;                        
  double _dir_thres;                        
  double _saliency_crn;
  int _rediter;                             
   const COM::Attribute *_cnstr_types;      
   COM::Attribute *_cnstr_dirs;             
   std::vector<std::set< Edge_ID> > _edges; 

protected: // Variables not currently modifiable via set_option
   
   //std::vector<std::vector<bool> >  _elems_to_check; ///< Check quality of these.
   std::vector<std::vector<bool> > _is_shared_node;   
   std::vector<std::vector<bool> > _is_shared_elem;   
   std::vector<std::vector<bool> > _is_phys_bnd_node; 
   std::vector<std::vector<bool> > _is_phys_bnd_elem; 
   std::vector<std::vector<bool> > _is_pane_bnd_node; 
   std::vector<std::vector<bool> > _is_pane_bnd_elem; 

   const COM::Window           *_usr_window;         
   COM::Window                 *_buf_window;          

   std::vector<MAP::Pane_dual_connectivity*> _dcs;   // Pane dual connectivities 
   
  bool                         _is_pmesh;            

 protected: // Variables which can be modified via set option.

   enum { MOP_COOKIE=762667};

  int _method; 
  // 0 == mesquite
  // 1 == mesquite, no ghosts
  // 2 == medial_quadric
  // 3 == error (default)

  int _wrapper; 
  // 0 == ShapeImprovement Wrapper (Feasible Newton)
  // 1 == CGWrapper (Conjugate Gradient)

  int _verb; 
  // 0 = no debugging info
  // 1 = routine entries
  // 2 = routine entries and exits
  // 3 = highest level in routine info
  // 4 = low level in routine info
  // 5 = 4 + high level MesqPane info
  // 6 = 5 + low level MesqPane info
  
   int _lazy; 

   float _tol; 
   // Smoother loops until either 
   //   worst element quality > tol 
   //   OR _niter iterations reached

   int _niter; 

   float _ctol; 
   // Parallel mesquite subcycles for convergence until either 
   //      (post communication quality)/(pre communication quality) > _tol
   //      OR _ncycle cycles performed.
   // In this context, quality refers to worst quality of all real elements
   //      containing shared nodes.
   // _ctol should be set to a value between 0 and 1

   int _ncycle; 

   int _cookie; 

   int _invert_tets; 

   int _invert_hexes; 

   float _maxdisp; 

   int _smoothfreq;         

  float _disp_thresh; // Threshold for smoothing trigger
};

MOP_END_NAMESPACE

#endif