/Users/mourrain/Devel/mmx/shape/include/shape/cell3d_algebraic_curve.hpp

Go to the documentation of this file.

/*****************************************************************************
 * M a t h e m a  g i x
 *****************************************************************************
 * Cell3d for Algebraic Curve
 * 2008-03-20
 * Bernard Mourrain & Julien Wintz
 *****************************************************************************
 *               Copyright (C) 2008 INRIA Sophia-Antipolis
 *****************************************************************************
 * Comments :
 ****************************************************************************/

# ifndef shape_cell3d_algebraic_curve_hpp
# define shape_cell3d_algebraic_curve_hpp

# include <realroot/Interval.hpp>
# include <realroot/polynomial_bernstein.hpp>
# include <realroot/Seq.hpp>
# include <shape/cell3d.hpp>
# include <shape/box_face.hpp>
# include <shape/topology.hpp>
# include <shape/algebraic_curve.hpp>
# include <shape/solver_implicit.hpp>

# define TMPL template<class C, class V>
# define AlgebraicCurve algebraic_curve<C,V>
# define SELF cell3d_algebraic_curve<C,V>

//====================================================================
namespace mmx { namespace shape {
//====================================================================

template<class C, class V= default_env>
struct cell3d_algebraic_curve : public cell3d<C,V> {
public:

    typedef topology<C,V>  Topology;
    typedef tpl3d<C,V>     Topology3d;
    typedef vertex<C,V>    Point;  

    typedef polynomial< double, with<Bernstein> > Polynomial;

    typedef cell3d<C,V>    Cell;  
    typedef cell3d<C,V>    Cell3d;
    typedef typename  Cell3d::BoundingBox  BoundingBox;

    typedef solver_implicit<C,V> Solver;

    cell3d_algebraic_curve(const SELF&);
    cell3d_algebraic_curve(const Seq<Polynomial>&, const BoundingBox& );
    cell3d_algebraic_curve(AlgebraicCurve *, const BoundingBox& bx);

    bool is_regular(void);
    bool is_active (void) const;

    virtual bool is_intersected (void) {return true;}

    template<class TOPOLOGY>
    bool inserted_regular_in(TOPOLOGY*) ;

    template<class TOPOLOGY> 
    bool inserted_singular_in(TOPOLOGY*) ;

    virtual bool insert_regular (Topology*) ;
    virtual bool insert_singular(Topology*) ;

    virtual void polygonise (Topology3d*) {};

    virtual int  subdivide(SELF*& left, SELF*& right) ;
    virtual void subdivide(Cell*& left, Cell*& right, int v, double s);

    int nbeq() const { return m_polynomials.size(); }
    Polynomial equation(int i=0) const { return  m_polynomials[i]; }

    struct along {
      int m_dir;
      along(int i):m_dir(i) {}
      bool operator()(Point *p1, Point* p2) {
        return (*p1)[m_dir]<(*p2)[m_dir];
      }
    };

  bool topology_regular(Topology*) { return true;}
    //    void check_cluster(Seq<Point*>& , double eps=0.00001);
    void check_insert (Point*, Seq<Point*>& , double eps=0.00001);

  public:
    Seq<Polynomial>   m_polynomials;
    AlgebraicCurve*   m_shape;
    int               m_reg;
};
  
//====================================================================
TMPL
SELF::cell3d_algebraic_curve(const SELF& c): 
  m_polynomials(c.m_polynomials), m_shape(c.m_shape), m_reg(-1)  {};
    
TMPL
SELF::cell3d_algebraic_curve(const Seq<Polynomial>& s, const BoundingBox& b): 
  Cell3d(b), m_polynomials(s), m_reg(-1)
{

}
  
TMPL
SELF::cell3d_algebraic_curve(AlgebraicCurve* cv, const BoundingBox& b): Cell3d(b), m_reg(-1) {
    Seq<mmx::GMP::rational> bx;
    bx<<as<mmx::GMP::rational>(b.xmin());
    bx<<as<mmx::GMP::rational>(b.xmax());
    bx<<as<mmx::GMP::rational>(b.ymin());
    bx<<as<mmx::GMP::rational>(b.ymax());
    bx<<as<mmx::GMP::rational>(b.zmin());
    bx<<as<mmx::GMP::rational>(b.zmax());
    
    Polynomial tmp;
    for(int i=0;i<cv->nbequation();i++) {
      tensor::bernstein<mmx::GMP::rational> polq(cv->equation(i).rep(),bx);
      let::assign(tmp.rep(),polq);
      m_polynomials<<tmp;
    }

    Solver::face_point(this->m_boundary, m_polynomials, Solver::north_face, b);
    Solver::face_point(this->m_boundary, m_polynomials, Solver::south_face, b);
    Solver::face_point(this->m_boundary, m_polynomials, Solver::west_face , b);
    Solver::face_point(this->m_boundary, m_polynomials, Solver::east_face , b);
    Solver::face_point(this->m_boundary, m_polynomials, Solver::front_face, b);
    Solver::face_point(this->m_boundary, m_polynomials, Solver::back_face , b);

    foreach(Point* p, this->m_boundary) std::cout<<"n "<<p->x()<<" "<<p->y()<<" "<<p->z() <<std::endl;

    //    check_cluster(this->m_boundary);
    //     Solver::singular(m_singular, m_polynomial, b); 
  }

  TMPL bool 
  SELF::is_active() const {
    //    std::cout<<"Is active "<<m_polynomials[0].rep()<<std::endl;

    if(!has_sign_variation(m_polynomials[0])) 
      return false;

    //    std::cout<<"Is active "<<m_polynomials[1]<<std::endl;

    if(!has_sign_variation(m_polynomials[1]))
      //m_polynomials[1].begin(),m_polynomials[1].end())) 
      return false;

       //    std::cout<<"Is active end"<<std::endl;
    
    return true; 
  }

  TMPL bool 
  SELF::is_regular() {
    if(this->m_singular.size()>1) return false;

    Polynomial 
      tx=diff(equation(0),1)*diff(equation(1),2)-diff(equation(1),1)*diff(equation(0),2);
    if(!has_sign_variation(tx)) 
      { m_reg=0; return true; }
    
    Polynomial 
        ty=diff(equation(0),0)*diff(equation(1),2)-diff(equation(1),0)*diff(equation(0),2);
    if(!has_sign_variation(ty)) 
      { m_reg=1; return true; }
    
    Polynomial 
      tz=diff(equation(0),0)*diff(equation(1),1)-diff(equation(1),0)*diff(equation(0),1);
    if(!has_sign_variation(tz)) 
      { m_reg=2; return true; }

    //    int n=this->m_boundary.size();
    //    if(n > 2) 
    return false;
    //    else {
    //      return true;
    //    }
  }
  //--------------------------------------------------------------------
  TMPL int
  SELF::subdivide(SELF*& Left, SELF*& Right) {
    int v; double s;
    this->split_position(v,s);
    this->subdivide((Cell*&)Left,(Cell*&)Right,v,s);
    return v;
  }
  //--------------------------------------------------------------------
  TMPL void
  SELF::subdivide(Cell*& Left, Cell*& Right, int v, double c) {

    SELF * left, * right;
   //  if (this->m_boundary.size()>2){
//       std::cout<<"             split cell with "<<this->m_boundary.size()<<" points"<<std::endl;
//       foreach(Point* p, this->m_boundary)
//      std::cout<<" "<<p->x()<<" "<<p->y()<<" "<<p->z()<<std::endl;
//     }
    if(v==0){
      c=(this->xmin()+this->xmax())/2;
      left = new SELF(m_polynomials, BoundingBox(this->xmin(),c, this->ymin(),this->ymax(), this->zmin(),this->zmax()));
      right= new SELF(m_polynomials, BoundingBox(c,this->xmax(), this->ymin(),this->ymax(), this->zmin(),this->zmax()));
    } else if(v==1) {
      c=(this->ymin()+this->ymax())/2;
      left = new SELF(m_polynomials, BoundingBox(this->xmin(),this->xmax(), this->ymin(),c, this->zmin(),this->zmax()));
      right= new SELF(m_polynomials, BoundingBox(this->xmin(),this->xmax(), c,this->ymax(), this->zmin(),this->zmax()));
    } else {
      c=(this->zmin()+this->zmax())/2;
      left = new SELF(m_polynomials, BoundingBox(this->xmin(),this->xmax(), this->ymin(),this->ymax(), this->zmin(),c));
      right= new SELF(m_polynomials, BoundingBox(this->xmin(),this->xmax(), this->ymin(),this->ymax(), c,this->zmax()));
    }

    for(int i=0;i<nbeq();i++) {
      tensor::split(left->m_polynomials[i], right->m_polynomials[i], v);
    }

    if(v==0) {
      foreach(Point* p, this->m_boundary) {
        if (Solver::east_face.is_valid(*p,left->boundingBox())) 
          left->m_boundary << p ;
        if (Solver::west_face.is_valid(*p,right->boundingBox())) 
          right->m_boundary << p ;
      }
    } else if (v==1) {
      foreach(Point* p, this->m_boundary) {
        if (Solver::north_face.is_valid(*p,left->boundingBox())) 
          left->m_boundary << p ;
        if (Solver::south_face.is_valid(*p,right->boundingBox())) 
          right->m_boundary << p ;
      }
    } else {
      foreach(Point* p, this->m_boundary) {
        if (Solver::front_face.is_valid(*p,left->boundingBox())) 
          left->m_boundary << p ;
        if (Solver::back_face.is_valid(*p,right->boundingBox())) 
          right->m_boundary << p ;
      }
    }

 

    Seq<Point*> tmp;
    // Points on vertical faces of back subcells
    if(v==0){
      Solver::face_point(tmp, left->m_polynomials, Solver::east_face, left->boundingBox());   
      foreach(Point *p, tmp){
        check_insert(p, left->m_boundary);
        check_insert(p, right->m_boundary);
      }
    } else if(v==1) {
      Solver::face_point(tmp, left->m_polynomials, Solver::north_face, left->boundingBox());   
      foreach(Point *p, tmp){
        check_insert(p, left->m_boundary);
        check_insert(p, right->m_boundary);
      }
    } else {
      Solver::face_point(tmp, left->m_polynomials, Solver::front_face, left->boundingBox());   
      foreach(Point *p, tmp){
        check_insert(p, left->m_boundary);
        check_insert(p, right->m_boundary);
      }
    }
//  if (tmp.size()>0){
//       std::cout<<"      insert  "<<tmp.size()<<" points"<<std::endl;
//       foreach(Point* p, tmp)
//      std::cout<<"            "<<p->x()<<" "<<p->y()<<" "<<p->z()<<std::endl;
//     }
//     if (left->m_boundary.size()>2){
//       std::cout<<"["<<left->xmin()<<" "<<left->xmax()<<" "<<left->ymin()<<" "<<left->ymax()<<" "<<left->zmin()<<" "<<left->zmax()<<"]"<<std::endl;
//       std::cout<<"      left cell with "<<left->m_boundary.size()<<" points"<<std::endl;
//       foreach(Point* p, left->m_boundary)
//      std::cout<<"            "<<p->x()<<" "<<p->y()<<" "<<p->z()<<std::endl;
//     }
//     if (right->m_boundary.size()>2){
//       std::cout<<"      right cell with "<<right->m_boundary.size()<<" points"<<std::endl;
//       std::cout<<"["<<right->xmin()<<" "<<right->xmax()<<" "<<right->ymin()<<" "<<right->ymax()<<" "<<right->zmin()<<" "<<right->zmax()<<"]"<<std::endl;
//       foreach(Point* p, right->m_boundary)
//      std::cout<<"            "<<p->x()<<" "<<p->y()<<" "<<p->z()<<std::endl;
//     }
    //    left ->m_boundary<<tmp;
    //    right->m_boundary<<tmp;

   
    Left=left; Right=right;
  }

  //--------------------------------------------------------------------
  TMPL template<class TOPOLOGY> bool SELF::inserted_regular_in(TOPOLOGY* t) {

    typedef typename TOPOLOGY::Edge Edge; 
    if(this->m_boundary.size()==0) return true;

//     if(this->m_boundary.size()==1) { 
//       std::cout<<"cell with 1 boundary point with reg "<<m_reg<<std::endl;      //      return false;
//     }
    if(this->m_boundary.size()==2) {
      t->insert(this->m_boundary[0]); 
      t->insert(this->m_boundary[1]); 
      t->insert(new Edge(this->m_boundary[0],this->m_boundary[1])); 
    } else {
      //      std::cout<<"cell with "<<this->m_boundary.size()<<" boundary points with reg "<<m_reg<<std::endl;
      BoundingBox bx=*this;
      //      std::cout<<"["<<bx.xmin()<<" "<<bx.xmax()<<" "<<bx.ymin()<<" "<<bx.ymax()<<" "<<bx.zmin()<<" "<<bx.zmax()<<"]"<<std::endl;

      std::sort(this->m_boundary.begin(), this->m_boundary.end(), along(m_reg) );

      foreach(Point* p, this->m_boundary) t->insert(p);
      
      for(unsigned i=0;i<this->m_boundary.size();i++) {
        unsigned j=i+1, i0=i;
        //      std::cout<<"i0 "<<i<<" "<<j<<std::endl;
        for(j=i+1; j<this->m_boundary.size() 
              && (this->m_boundary[i0] == this->m_boundary[j]);i++, j++) ;
        if(j<this->m_boundary.size())
          t->insert(new Edge( this->m_boundary[i0], this->m_boundary[j])); 
        //      std::cout<<"i0 "<<i<<" "<<j<<std::endl;
      }
    }

    return true; 
  }
//--------------------------------------------------------------------
TMPL bool SELF::insert_regular(Topology* t) {
  return this->inserted_regular_in<Topology>(t);
}
//--------------------------------------------------------------------
//this insert_singular is used for detecting the singularities in the computed topology from tpl3d
TMPL template<class TOPOLOGY> bool SELF::inserted_singular_in(TOPOLOGY* t) { 
  if(this->m_boundary.size()>0) { 
    t->insert((BoundingBox *)this);
    foreach(Point* p, this->m_boundary)
      t->insert(p);
  }       
  return true; 
} 
//--------------------------------------------------------------------
TMPL bool SELF::insert_singular(Topology* t) { 
  return this->inserted_singular_in<Topology>(t);
}
//--------------------------------------------------------------------
//  TMPL void
//   SELF::check_cluster(Seq<Point*>& l, double eps) {
//     //    Seq<Point*> res;
//     for(unsigned i=0; i<l.size(); i++) {
//       unsigned j=0;
//       for(j=0; j<i && distance(*l[j],*l[i])>eps;j++) ;
//       if(j<i) l[j]=l[i];
//     }
//   } 

  TMPL void
  SELF::check_insert(Point* p, Seq<Point*>& l, double eps) {
    //    l<<p; return ;
    unsigned i=0;
    for(i=0; i<l.size() && distance(*l[i],*p)>eps; i++) ;
    if(i==l.size()) l<<p;
  }


} ; // namespace shape
} ; // namespace mmx

# undef TMPL
# undef Point
# undef Cell
# undef AlgebraicCurve 
# undef SELF
# undef Solver 
# endif // SHAPE_IMPLICITCELL_H