synaps/topology/inter_points.h

/*******************************************************************
 *   This file is part of the source code of the SYNAPS kernel.
 *   Author(s): B. Mourrain, GALAAD, INRIA
 ********************************************************************/
#ifndef synaps_topology_inter_points_h
#define synaps_topology_inter_points_h
//===================================================================
#include <synaps/init.h>
#include <synaps/topology/point_graph.h>
//===================================================================

__BEGIN_NAMESPACE_SYNAPS 
// for sorting
struct XMAJOR
{
  template<class T>
  static bool less(const T& a, const T& b)
  {
    return (a->pnt[0]<b->pnt[0]);
  }
  template<class T>
  bool operator()(const T& a, const T& b)
  {
    return less(a,b);
  }
};
struct YMAJOR
{
  template<class T>
  static bool less(const T& a, const T& b)
  {
    return (a->pnt[1]<b->pnt[1]);
  }
  template<class T>
  bool operator()(const T& a, const T& b)
  {
    return less(a,b);
  }
};
struct ZMAJOR
{
  template<class T>
  static bool less(const T& a, const T& b)
  {
    return (a->pnt[2]<b->pnt[2]);
  }
  template<class T>
  bool operator()(const T& a, const T& b)
  {
    return less(a,b);
  }
};

//--------------------------------------------------------------------
struct XPT
{
  template<class T>
  static bool less(const T& a, const T& b)
  {
    return (a[0]<b[0]);
  }
  template<class T>
  bool operator()(const T& a, const T& b)
  {
    return less(a,b);
  }
};
struct YPT
{
  template<class T>
  static bool less(const T& a, const T& b)
  {
    return (a[1]<b[1]);
  }
  template<class T>
  bool operator()(const T& a, const T& b)
  {
    return less(a,b);
  }
};
struct ZPT
{
  template<class T>
  static bool less(const T& a, const T& b)
  {
    return (a[2]<b[2]);
  }
  template<class T>
  bool operator()(const T& a, const T& b)
  {
    return less(a,b);
  }
};
//--------------------------------------------------------------------
template<class point_t>
struct XYPLANE
{
  point_t center;

  XYPLANE(const point_t& pt): center(pt){};

  template<class T>
  bool less(const T& a, const T& b)
  {
    return ((a[0]-center[0]) * (b[1]-center[1]) <= (b[0]-center[0]) * (a[1]-center[1]));
  }
  template<class T>
  bool operator()(const T& a, const T& b)
  {
    return less(a,b);
  }
};


template<class point_t>
struct XZPLANE
{
  point_t center;

  XZPLANE(const point_t& pt): center(pt){};

  template<class T>
  bool less(const T& a, const T& b)
  {
    return ((a[0]-center[0]) * (b[2]-center[2]) <= (b[0]-center[0]) * (a[2]-center[2]));
  }
  template<class T>
  bool operator()(const T& a, const T& b)
  {
    return less(a,b);
  }
};

template<class point_t>
struct YZPLANE
{
  point_t center;

  YZPLANE(const point_t& pt): center(pt){};

  template<class T>
  bool less(const T& a, const T& b)
  {
    return ((a[1]-center[1]) * (b[2]-center[2]) <= (b[1]-center[1]) * (a[2]-center[2]));
  }
  template<class T>
  bool operator()(const T& a, const T& b)
  {
    return less(a,b);
  }
};
//--------------------------------------------------------------------
#define _is_singular(x) ((x) < 1e-6 && (x) > -1e-6)
#define _is_verysmall(x) (std::abs(x)<1e-3)
template<class POINT>
bool same_point(const POINT & a, const POINT & b)
{
  return (_is_singular(a[0]-b[0]) && _is_singular(a[1]-b[1]) && _is_singular(a[2]-b[2]));
}
//--------------------------------------------------------------------
// maximum number of faces that will be computed for each cell
#define MAX_FACEIDX 3
// index of a cell to be computed per computation for a face
static int _cellidx[] =  { 0, 0, 0,  1, 1, 1,  2, 2, 2,   3, 3, 3,
                           4, 4, 4,  5, 5, 5,  6, 6, 6,   7, 7, 7};
// index of a face that needed to compute the intersections; a cell has at most
// 2 faces needed to be computed; if a cell gets less than 3 faces to be computed,
// a -1 index will hold the place
static int _faceidx[] =  { 1, 3, 5,  3, 5,-1,  1, 5,-1,   5,-1,-1,
                           1, 3,-1,  3,-1,-1,  1,-1,-1,  -1,-1,-1};
// index of a subcell that sharing a faces, for each face indicated by _faceidx
static int _shareidx[] = { 1, 2, 4,  3, 5,-1,  3, 6,-1,  7,-1,-1,
                           5, 6,-1,  7,-1,-1,  7,-1,-1, -1,-1,-1};
// index of faces that the subcells inherits from their parent cells
static int _inheritidx[]={ 0, 2, 4,  1, 2, 4,  0, 3, 4,  1, 3, 4,
                           0, 2, 5,  1, 2, 5,  0, 3, 5,  1, 3, 5};

// index of axis that represent the plane (of a cell face), see also inter->nPlane
static int plane_axisidx[] = {1, 2,
                              0, 2,
                              0, 1}; // first and second indices for x-, y-, or z-major intersection
// index of subcell that should inherit the intersection from its parent cell (after various test)
static int subcellidx[]    = {0, 1, 4, 5, 2, 3, 6, 7,  // x-major (intersection lies on YZ plane)
                              0, 2, 4, 6, 1, 3, 5, 7,  // y-major (intersection lies on XZ plane)
                              0, 4, 2, 6, 1, 5, 3, 7}; // z-major (intersection lies on XY plane)


// enum PLANE indices
enum PLANES
{
    PLANEYZ = 0,
    PLANEXZ = 1,
    PLANEXY = 2
};

//====================================================================
// custom structures

// intersection3d: the class for representing the intersections between the curve and cell facets
//                 and a variety of information that is related to those intersections
template < class C = double >
struct intersection3d
{
    // type definitions:

    typedef topology::point_graph<C>   pointgraph_t;
    typedef topology::point<C>         point_t;


    // constructors:

    intersection3d() : nMulti(1), nPlane(-1), nCount(0) {}
    intersection3d(const point_t& p, int plane = -1, int nMul = 1) : nMulti(nMul), nPlane(plane), nCount(0) {pnt = p;}
    intersection3d(const intersection3d& i) : nCount(0) {assign(i);}


    // member data:

    point_t pnt;    // the intersection point
    point_t der;    // kth order derivative if derivative k-1 of curve at pnt vanishes;
    int nMulti;     // its multiplicity
    int nPlane;     // the id of the plane (yz=0, xz=1, xy=2)

    // ephemeral data, won't be duped by assign
    int nFlag;      // temporary flags, for buiding graphs, house keeping etc.
    int nCount;     // an incremental count for traversal
    int nDepth;

    // member functions:

    void assign(const intersection3d& i)
    {
        pnt = i.pnt;
        der = i.der;
        nMulti = i.nMulti;
        nPlane = i.nPlane;
    }

    intersection3d&  operator = (const intersection3d& i)
    {
        assign(i);
        return *this;
    }
  C operator[](unsigned i) {return pnt[i];}


};
//====================================================================
// for comparison between two 'intersections'
int ptcmp(const shared_object<intersection3d<double>*>& p1, 
          const shared_object<intersection3d<double>*>& p2)
{
    return ptcmp(rep(p1)->pnt,rep(p2)->pnt);
}

//===================================================================
__END_NAMESPACE_SYNAPS
#endif //synaps_topology_inter_points_h