Manifold_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: Manifold_2.h,v 1.18 2008/12/06 08:43:23 mtcampbe Exp $

//=======================================================================
//  This file defines the data structures for accessing a surface mesh.
//  The data structures mimic the halfedge data structure, but with the
//  extension to support partitioned meshes. Typically, the user should 
//  use Node and Halfedge to access the entities of a surface mesh.
//
//  Author: Xiangmin Jiao
//  Last modified: Oct 11, 2004
//=======================================================================

#ifndef __MANIFOLD_2_H_
#define __MANIFOLD_2_H_

#include "surfbasic.h"
#include "../Rocmap/include/Simple_manifold_2.h"

namespace MAP { class Pane_communicator; };

SURF_BEGIN_NAMESPACE

using MAP::Simple_manifold_2;

class Window_manifold_2;
class Halfedge;
class Node;
typedef MAP::Facet_ID Edge_ID;

enum Access_Mode { REAL_PANE, WHOLE_PANE, ACROSS_PANE };

class Pane_manifold_2 : public Simple_manifold_2 {
  friend class Window_manifold_2;
  friend class Node;
  friend class Halfedge;
public:
  Pane_manifold_2();

  explicit Pane_manifold_2( const COM::Pane *p, 
                            const Pane_manifold_2 *parent=NULL);

  ~Pane_manifold_2();

  inline const Halfedge &get_counterpart( Edge_ID eID) const;

  inline const Vector_3<Real> &get_bd_normal( Edge_ID) const;

  inline int get_bd_flag( Edge_ID) const;

  inline bool get_bdedge_flag( Edge_ID) const;

  inline Halfedge get_opposite_edge( Edge_ID, Access_Mode mode) const;

  Halfedge get_prev_edge( Edge_ID, Access_Mode mode) const;

  Halfedge get_next_edge( Edge_ID, Access_Mode mode) const;

  inline Node get_primary( int vID, Access_Mode mode) const;

  inline bool is_primary( int vID, Access_Mode mode) const;

  inline bool is_border_edge( Edge_ID eID, Access_Mode mode) const;

  inline bool is_physical_border_edge( Edge_ID eID) const
  { return eID.is_border() && _cnt_pn[get_border_edgeID(eID)-1]==0; }

  int size_of_nodes( Access_Mode mode) const;

  int size_of_edges( Access_Mode mode) const;

  int size_of_faces( Access_Mode mode) const {
    switch ( mode) {
    case REAL_PANE:
    case ACROSS_PANE: return size_of_real_faces(); 
    default: return Simple_manifold_2::size_of_faces();
    }
  }

  int size_of_triangles( Access_Mode mode) const {
    switch ( mode) {
    case REAL_PANE:
    case ACROSS_PANE: return size_of_real_triangles(); 
    default: return Simple_manifold_2::size_of_triangles();
    }
  }

  int size_of_quadrilaterals( Access_Mode mode) const {
    switch ( mode) {
    case REAL_PANE:
    case ACROSS_PANE: return size_of_real_quadrilaterals(); 
    default: return Simple_manifold_2::size_of_quadrilaterals();
    }
  }

protected:
  void convert_pconn_edge2ghost( const COM::Attribute *pconn_e, 
                                 COM::Attribute *pconn_g);


  int get_bnode_index( int vID) const {
    Edge_ID eID = get_incident_real_edge( vID);

    if ( !eID.is_border()) return -1;

    int bid = eID.eid()-1;
    if ( bid>=int(_beIDs.size())) return bid-_size_ghost_borders;
    else return bid; 
  }

  /* The following data memters are for storing the counterpart of each
   * edge on pane boundaries, for the across-pane access mode. */
  std::vector<int>             _cnt_pn;

  //< Stores pane-id of the counterparts of border edges.
  //< If owner pane is remote, then save the negative ID. If there is no
  //< counterpart (i.e., edge is on physical boundary), then save 0.

  std::vector<Halfedge>        _bd_cnt;   //< Counterpart edges in ACROSS_PANE mode
  std::vector<Vector_3<Real> > _bd_nrms;  //< Normal of opposite face
  std::vector<int >            _bd_flgs;  //< Flag of opposite face
  std::vector<char >           _bd_bm;    //< Bitmap of opposite edge
  std::vector<Node>            _nd_prm;   //< Primary nodes in ACROSS_PANE mode.
};

class Window_manifold_2 {
public:
  typedef std::vector< Pane_manifold_2>::iterator PM_iterator;
  typedef std::vector< Pane_manifold_2>::const_iterator PM_const_iterator;

  Window_manifold_2() : _buf_window(NULL), _cc(NULL), _pconn_nb(0) {}

  explicit Window_manifold_2(const COM::Attribute *pmesh)
    : _buf_window(NULL), _cc(NULL), _pconn_nb(0) { init(pmesh); }

  ~Window_manifold_2();

  COM::Window *window() { return _buf_window; }
  const COM::Window *window() const { return _buf_window; }

  //\}

  Pane_manifold_2 *get_pane_manifold( int pid)
  { return &_pms[ _pi_map[pid]]; }

  const Pane_manifold_2 *get_pane_manifold( int pid) const
  { return &_pms[ _pi_map.find(pid)->second]; }

  int size_of_panes() const { return _pms.size(); }

  PM_iterator pm_begin() { return _pms.begin(); }
  PM_const_iterator pm_begin() const { return _pms.begin(); }

  PM_iterator pm_end()   { return _pms.end(); }
  PM_const_iterator pm_end() const { return _pms.end(); }

  //\}

  void init_communicator();

  void reduce_on_shared_nodes( COM::Attribute *attr, MPI_Op op);

  //\}

  void compute_normals( COM::Attribute *normals, 
                        int scheme = E2N_ANGLE,
                        bool to_normalize = true);
  
  void update_bd_normals( const COM::Attribute *normals,
                          bool to_normalize = false);

  void update_bd_flags( const COM::Attribute *flags);

  void update_bdedge_bitmap( const COM::Attribute *bitmap);

  void accumulate_bd_values( const COM::Attribute *vals);

  // Compute mean-curvature normals and their Laplace-Beltrami operator.
  void compute_mcn( COM::Attribute *mcn_in, COM::Attribute *lbmcn_in);

  void elements_to_nodes( const COM::Attribute *evals,
                          COM::Attribute *nvals,
                          const int scheme = E2N_ONE,
                          const COM::Attribute *ews = NULL,
                          COM::Attribute *nws = NULL,
                          const int tosum=false);

  void shortest_edge_length( COM::Attribute *lens);

  void perturb_mesh( double range);
  //\}

  int size_of_nodes( Access_Mode mode) const;

  int size_of_faces( Access_Mode mode) const;

  int size_of_triangles( Access_Mode mode) const;

  int size_of_quadrilaterals( Access_Mode mode) const;
  
  int size_of_edges( Access_Mode mode) const;

  int pconn_nblocks() const { return _pconn_nb; }

  void serialize_window( COM::Window *outwin) const;

protected:
  void init( const COM::Attribute *pmesh);

  void determine_counterparts( const COM::Attribute *pconn_e);
  
  void determine_primaries( const COM::Attribute *pconn_n);

  static void compute_shortest_edgelen_elements( COM::Attribute *lens);

  void compute_shortest_edgelen_nodes( COM::Attribute *lens);

  void compute_nodal_normals( COM::Attribute *nrm,
                              int scheme = E2N_ANGLE,
                              bool to_normalize = true,
                              const COM::Attribute *weights = NULL);

  // Assign nodal IDs. Called by serialize_window()
  void assign_global_nodeIDs( std::vector<std::vector<int> > &gids) const;

  //\}

public:
  static MPI_Op OP_MAX, OP_MIN, OP_SUM, OP_PROD, OP_MAXABS, OP_DIFF,
    OP_BOR, OP_BAND, OP_LOR, OP_LAND;

protected:
  // A buffer window created by the Manifold class by inheriting from
  // a user window when the init() function is called. in general, 
  // one must inherit a user attribute onto this window.
  COM::Window                     *_buf_window;
  std::vector< Pane_manifold_2>    _pms;
  // Create a mapping from pane ids to the corresponding indices in _pms
  std::map< int, int>              _pi_map;
  MAP::Pane_communicator          *_cc;       // Pane communicator.
  int                              _pconn_nb; // Number of blocks of pconn
  //\}
};

class Node {
public:
  Node() : _pm(NULL), _vID(0), _mode( ACROSS_PANE) {}

  Node( const Pane_manifold_2 *pm, int vid, Access_Mode mode) 
    : _pm( pm), _vID( vid), _mode( mode) {}

  inline Halfedge halfedge() const;

  const Pane_manifold_2 *pane_manifold() const { return _pm; }

  const COM::Pane *pane() const 
  { if ( _pm) return _pm->pane(); else return NULL; }

  int  id() const { return _vID; }

  const Point_3<Real> &point() const 
  { return (const Point_3<Real>&)_pm->pane()->coordinates()[ 3*_vID-3]; }
  
  const void *addr( const COM::Attribute *attr) const;

  const COM::Attribute *attr( const COM::Attribute *a) const { 
    COM_assertion_msg( a, "Unexpected NULL pointer");
    return a->window()->pane( _pm->pane()->id()).attribute(a->id()); 
  }

  inline bool is_border() const;

  bool is_isolated() const { return _pm->is_isolated_node( _vID); }

  bool is_primary() const { return _pm->is_primary( _vID, _mode); }

  Node get_primary() const { return _pm->get_primary( _vID, _mode); }

  bool coincide( const Node &n) const
  { return _pm->get_primary( _vID, _mode) == n._pm->get_primary( n._vID, n._mode); }

  bool operator==( const Node &n) const
  { return _pm == n._pm && _vID == n._vID; }

  bool operator!=( const Node &n) const
  { return _pm != n._pm || _vID != n._vID; }

  bool operator<( const Node &n) const {
    return _pm->pane()->id() < n._pm->pane()->id() ||
      _pm == n._pm && _vID < n._vID; 
  }

protected:
  const Pane_manifold_2 *_pm;      // Reference to Pane_manifold_2
  int                    _vID;     // Node id
  Access_Mode            _mode;    // Access mode
};

class Halfedge {
public:
  Halfedge() : _pm(NULL) {}
  Halfedge( const Pane_manifold_2 *pm, Edge_ID e, Access_Mode mode) 
    : _pm( pm), _eID( e), _mode(mode) {}

  Halfedge opposite() const 
  { return _pm->get_opposite_edge( _eID, _mode); }

  Halfedge prev() const { return _pm->get_prev_edge( _eID, _mode); }

  Halfedge next() const { return _pm->get_next_edge( _eID, _mode); }

  Node   origin() const
  { return _pm->get_primary( _pm->get_origin( _eID), _mode); }

  Node   destination() const
  { return _pm->get_primary( _pm->get_destination( _eID), _mode); }

  bool is_border() const 
  { return _pm->is_border_edge( _eID, _mode); }

  bool is_physical_border() const 
  { return _pm->is_physical_border_edge( _eID); }

  Vector_3<Real> tangent() const
  { return destination().point()-origin().point(); }

  Vector_3<Real> normal() const;

  Vector_3<Real> deformed_normal( const COM::Attribute *disp) const;

  void edge_center( Point_3<Real> &mid) const 
  { edge_average( _pm->pane()->attribute( COM::COM_NC), 
                  reinterpret_cast<Vector_3<Real> &>(mid)); }

  void face_center( Point_3<Real> &cnt) const 
  { face_average( _pm->pane()->attribute( COM::COM_NC), 
                  reinterpret_cast<Vector_3<Real> &>(cnt)); }

  template < class Value>
  void edge_average( const COM::Attribute *attr, Value &val) const {
    val = 0.5*(*reinterpret_cast<const Value *>(destination().addr(attr))+
               *reinterpret_cast<const Value *>(origin().addr(attr))); 
  }
  
  template < class Value>
  void face_average( const COM::Attribute *attr,  Value &val) const {
    // If this is a border halfedge, then return its edge center.
    if ( is_border()) { edge_average( attr, val); return; }
  
    // Otherwise, compute the center
    val = Value(0);
    
    int count=0;
    Halfedge h=*this, h0=h;

    do {
      val += *reinterpret_cast<const Value*>(h.origin().addr( attr));
      ++count;
    } while ( (h=h.next()) != h0);
    
    val/=count;
  }
  
  double dihedral_angle() const;

  const Pane_manifold_2 *pane_manifold() const { return _pm; }

  const COM::Pane *pane() const 
  { if ( _pm) return _pm->pane(); else return NULL; }

  Edge_ID id() const { return _eID; }

  const void *addr(  const COM::Attribute *attr) const;

  const COM::Attribute *attr( const COM::Attribute *a) const { 
    COM_assertion_msg( a, "Unexpected NULL pointer");
    return a->window()->pane( _pm->pane()->id()).attribute(a->id()); 
  }

  bool operator==( const Halfedge &h) const 
  { return _pm==h._pm && _eID==h._eID; }

  bool operator!=( const Halfedge &h) const 
  { return _pm!=h._pm || _eID!=h._eID; }

  bool operator<( const Halfedge &h) const {
    int pid1=_pm->pane()->id(), pid2=h._pm->pane()->id();
    return pid1<pid2 || pid1==pid2 && _eID<h._eID; 
  }

protected:
  const Pane_manifold_2 *_pm;      // Reference to Pane_manifold_2
  Edge_ID                _eID;     // Edge ID
  Access_Mode            _mode;    // Access mode
};

const Halfedge &Pane_manifold_2::get_counterpart( Edge_ID eID) const 
{
  COM_assertion( eID.is_border());
  return _bd_cnt[ eID.eid()-1];
}

const Vector_3<Real> &Pane_manifold_2::get_bd_normal( Edge_ID eID) const {
  COM_assertion_msg( eID.is_border() && !_bd_nrms.empty(),
                     "get_bd_normal must be called on border edges and be called after update_bd_normals() has been called.");
  return _bd_nrms[ eID.eid()-1];
}

int Pane_manifold_2::get_bd_flag( Edge_ID eID) const {
  COM_assertion_msg( eID.is_border() && !_bd_flgs.empty(),
                     "get_bd_flag must be called on border edges and be called after update_bd_flags() has been called.");
  return _bd_flgs[ eID.eid()-1];
}

bool Pane_manifold_2::get_bdedge_flag( Edge_ID eID) const {
  COM_assertion_msg( eID.is_border() && !_bd_bm.empty(),
                     "get_bdedge_flag must be called on border edges and be called after update_bdedge_bitmap() has been called.");
  return _bd_bm[ eID.eid()-1];
}

// Obtain the opposite edge in a given access mode
Halfedge Pane_manifold_2::
get_opposite_edge( Edge_ID eID, Access_Mode mode) const {
  Edge_ID opp;
  switch ( mode) {
  case REAL_PANE:  
    opp = get_opposite_real_edge( eID); break;
  case WHOLE_PANE: 
    opp = Simple_manifold_2::get_opposite_edge( eID); break;
  case ACROSS_PANE: {
    opp=get_opposite_real_edge( eID);
   
    if (  opp.is_border()) {
      Halfedge h = get_counterpart(opp);
      if ( h.pane_manifold()) return h;
    }

    break;
  }
  default:  COM_assertion_msg(false, "Should never reach here");
  }
  return Halfedge( this, opp, mode);
}

Node Pane_manifold_2::get_primary( int vID, Access_Mode mode) const 
{
  if ( mode != ACROSS_PANE) return Node(this, vID, mode);

  int bid = get_bnode_index( vID);

  COM_assertion( bid<int(_nd_prm.size()));
  if ( bid >=0 && _nd_prm[ bid].pane_manifold())
    return _nd_prm[ bid];
  else
    return Node( this, vID, ACROSS_PANE);
}

bool Pane_manifold_2::is_primary( int vID, Access_Mode mode) const 
{
  if ( mode != ACROSS_PANE) return true;

  int bid = get_bnode_index( vID);

  COM_assertion( bid<int(_nd_prm.size()));
  return bid <0 || !_nd_prm[ bid].pane_manifold();
}

// Is the edge on the boundary of a given mode?
bool Pane_manifold_2::is_border_edge( Edge_ID eID, Access_Mode mode) const 
{ 
  switch ( mode) {
  case REAL_PANE:  
    return is_real_border_edge( eID);
  case WHOLE_PANE: 
    // For whole pane, an edge is border only if it is pure real border edge.
    return is_pure_real_border_edge( eID);
  case ACROSS_PANE: 
    return eID.is_border() && get_counterpart( eID).pane_manifold() == NULL;
  default:  COM_assertion_msg(false, "Should never reach here");
    return false;
  }
}

bool Node::is_border() const { return halfedge().is_border(); }

Halfedge Node::halfedge() const 
{
  Edge_ID eID;

  switch ( _mode) {
  case REAL_PANE:
    eID = _pm->get_incident_real_edge( _vID); break;
  case WHOLE_PANE:
    eID = _pm->get_incident_edge( _vID); break;
  default: {
    eID =_pm->get_incident_real_edge( _vID);
    // the current node (vertex) is in general the origin of the edge or 
    // the edge center, unless the node is isolated. 

    if ( eID.is_border()) {
      Halfedge h=_pm->get_counterpart( eID);
      if ( h.pane_manifold()) return h;
    }
  }
  }

  return Halfedge(_pm, eID, _mode); 
}

Vector_3<Real> get_normal( const Halfedge &h, const Vector_2<Real> &nc);

Vector_3<Real> get_deformed_normal( const Halfedge &h, 
                                    const Vector_2<Real> &nc, 
                                    const COM::Attribute *disp);

inline Vector_3<Real> get_tangent( const Halfedge &h) 
{ return h.tangent(); }


SURF_END_NAMESPACE

#endif