include/shape/cell2d_algebraic_curve.hpp

Go to the documentation of this file.

/*****************************************************************************
 * M a t h e m a g i x
 *****************************************************************************
 * cell2d_algebraic_curve
 * 2008-03-20
 * Julien Wintz & Bernard Mourrain
 *****************************************************************************
 *               Copyright (C) 2008 INRIA Sophia-Antipolis
 *****************************************************************************
 * Comments :
 ****************************************************************************/
# ifndef shape_cell2d_algebraic_curve_hpp
# define shape_cell2d_algebraic_curve_hpp

# include <realroot/Interval.hpp>
# include <realroot/ring_bernstein_tensor.hpp>
# include <realroot/Seq.hpp>
# include <shape/point.hpp>
# include <shape/algebraic_curve.hpp>
# include <shape/solver_implicit.hpp>
# include <shape/cell2d.hpp>
# include <stack>

# define STACK std::stack<Point*>
# define TMPL  template<class C, class V>
# define TMPL1 template<class V>
# define SELF  cell2d_algebraic_curve<C,V>
# define REF   REF_OF(V)

# undef Cell_t
//====================================================================
namespace mmx { namespace shape {
//====================================================================
TMPL struct cell2d_subdivisor
{
    template<class CELL> static void
            subdivide(CELL & cl, CELL*&, CELL*&, int v, double s);
};

//--------------------------------------------------------------------
TMPL
template<class CELL> void 
cell2d_subdivisor<C,V>::subdivide(CELL& cl, CELL*& left, CELL*& right, int v, double s) {

    typedef typename topology<C,V>::Point Point;
    typedef solver_implicit<C,V> Solver;


    double eps =Approximate().eps;
    if(v==0) {//direction v=0
        left = new CELL(cl);   left ->set_xmax(s);
        right= new CELL(cl);   right->set_xmin(s);

        tensor::split(left->m_polynomial, right->m_polynomial, v);

        foreach(Point* p,cl.m_special) {

            if(p->x() < s+eps)
                left ->m_special << p ;
            //  else
            if(p->x() > s-eps)
                right->m_special << p ;
        }


    } else {//direction v==1

        left = new CELL(cl);  left->set_ymax(s);
        right= new CELL(cl);  right->set_ymin(s);

        tensor::split(left->m_polynomial, right->m_polynomial, v);

        foreach(Point* p,cl.m_special) {
            if(p->y() < s+eps)
                left->m_special << p ;
            //  else
            if(p->y() > s-eps)
                right->m_special << p ;
        }
    }

}
//====================================================================
template<class C, class V=default_env>
class cell2d_algebraic_curve : public cell2d<C,REF> {
public:
    typedef C                                          Scalar;
    typedef topology<C,REF>                            Topology;
    typedef typename Topology::Point                   Point;
    typedef Seq<Point*>                                Points;
    typedef typename use<point_def,C,V>::Point         Vector;
    typedef typename topology<C,REF>::Edge             Edge;
    typedef cell<C,REF>                                Cell;
    typedef cell<C,REF>                              CellBase;
    typedef algebraic_set<C,REF>                       AlgebraicSet;
    typedef algebraic_curve<C,REF>                     AlgebraicCurve;
    typedef polynomial< Interval<C>, with<Bernstein> > Polynomial;
    typedef bounding_box<C,REF>                        BoundingBox;

    typedef solver_implicit<C,V>                       Solver;

    cell2d_algebraic_curve(const Polynomial &, const BoundingBox&, bool inter=true);
    cell2d_algebraic_curve(const char*, const BoundingBox&);
    cell2d_algebraic_curve(const typename AlgebraicSet::Polynomial& , const BoundingBox&);
    cell2d_algebraic_curve(AlgebraicCurve*, const BoundingBox &x);
    cell2d_algebraic_curve(const SELF& cl)
        : cell2d<C,REF>(cl.boundingBox()), m_polynomial(cl.m_polynomial) {}

    Polynomial equation(void) { return  m_polynomial; }

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

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

    Vector  gradient(const Point& p);
    C       center_value() { return m_polynomial(0.5,0.5,0.); }

public:
    Polynomial      m_polynomial;
    Seq<Point *>    m_special;
    Point*          m_center;
};

//====================================================================
TMPL inline void 
print(SELF* cv) {
    typedef typename SELF::Point Point;
    std::cout <<"\nBox       :  ["<<cv->xmin()<<","<<cv->xmax()<<"] *"
             <<" ["<<cv->ymin()<<","<<cv->ymax()<<"]";
    std::cout <<"\npoint(s) s: ";
    foreach(Point* p, cv->s_intersections)
        std::cout <<" ["<<p->x()<<" "<<p->y()<<"]";
    std::cout << "\npoint(s) e: ";
    foreach(Point* p, cv->e_intersections)
        std::cout <<" ["<<p->x()<<" "<<p->y()<<"]";
    std::cout << "\npoint(s) n: ";
    foreach(Point* p, cv->n_intersections)
        std::cout <<" ["<<p->x()<<" "<<p->y()<<"]";
    std::cout << "\npoint(s) w: ";
    foreach(Point* p, cv->w_intersections)
        std::cout <<" ["<<p->x()<<" "<<p->y()<<"]";

    std::cout<<std::endl;

    std::cout << "Point(s) sing: ";
    foreach(Point* p, cv->m_special)
        std::cout <<" ["<<p->x()<<" "<<p->y()<<"]";
    std::cout<<std::endl;
}
//--------------------------------------------------------------------
TMPL
SELF::cell2d_algebraic_curve(const Polynomial& pol, const BoundingBox& bx, bool inter): cell2d<C,REF>(bx), m_polynomial(pol)
{
    //Solver::extremal(this->m_special, m_polynomial, bx);
}
//--------------------------------------------------------------------  
TMPL
SELF::cell2d_algebraic_curve(const typename AlgebraicSet::Polynomial& p, const BoundingBox & b): cell2d<C,REF>(b)
{
    //std::cout<<"Alegbraic cell: "<<p<<std::endl;
    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());
    
    tensor::bernstein<mmx::GMP::rational> polq(p.rep(),bx);

    let::assign(m_polynomial.rep(),polq);

    Solver::extremal(this->m_special, p, b);

    // foreach(Point* p,this-> m_special)
    //   std::cout<<"extremal point: "<<p->x()<<", "<<p->y() <<std::endl;

    //std::cout<<"Found "<<this->m_special.size()<<" extremal points" <<std::endl;
}

//--------------------------------------------------------------------
TMPL
SELF::cell2d_algebraic_curve(const char* s, const BoundingBox & b): cell2d<C,REF>(b)
{
    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());
    
    typename AlgebraicSet::Polynomial pol(s);
    tensor::bernstein<mmx::GMP::rational> polq(s,bx);

    let::assign(m_polynomial.rep(),polq);

    Solver::extremal(this->m_special, pol, b);
    foreach(Point* p,this-> m_special)
        std::cout<<"**extremal point: "<<p->x()<<" "<<p->y() <<std::endl;
}

//--------------------------------------------------------------------  
TMPL
SELF::cell2d_algebraic_curve(AlgebraicCurve* cv, const BoundingBox & b)
{
    new (this) SELF(cv->equation(),b );
}

//--------------------------------------------------------------------
TMPL bool 
SELF::is_regular() {
    // Assumption: cells with extremal points (dx=dy=0) that do not belong on the curve do not intersect the curve (assumes deep enough subdivision)

    if(this->m_special.size()>1) return false;

    if (this->m_special.size()==0) {
        if(!has_sign_variation(m_polynomial)) return true;

        Polynomial dx= diff(m_polynomial,0);
        if(!has_sign_variation(dx)) return true;

        Polynomial dy= diff(m_polynomial,1);
        if(!has_sign_variation(dy)) return true;

        return false;
//        int n= this->nb_intersect();
//        return   ( n==2 );
        //return   ( n%2==0);
        //return   (n == 4 || n == 2 || n == 0 );
    }

    //return true;
    //--------------------------------------------------------------------
    if(this->m_special.size()==1) {
        //    Polynomial
        //      dxf = diff(m_polynomial,0),
        //      dyf = diff(m_polynomial,1);
        //       std::cout<< "__________________"<<std::endl;
        //       std::cout<< m_polynomial;std::cout<< std::endl;
        //       std::cout<< "__________________"<<std::endl;
        //       std::cout<< dxf;std::cout<< std::endl;
        //       std::cout<< "__________________"<<std::endl;
        //       std::cout<< dyf;std::cout<< std::endl;
        //       std::cout<< "__________________"<<std::endl;

        //      int d;
        //      d = 1-2*topological_degree(dxf,dyf);
        //      d = 0;
        return false; //(d==ni);
    }
    
    return true;
}

//--------------------------------------------------------------------
TMPL bool 
SELF::is_active() {
    return (has_sign_variation(this->m_polynomial));
}
//--------------------------------------------------------------------
TMPL void
SELF::split_position(int& v, double& s) {
  double sx = (this->xmax()-this->xmin());
  double sy = (this->ymax()-this->ymin());
  if(sx<sy) {
    v=1;
    s=(this->ymax()+this->ymin())/2;
  } else {
    v=0;
    s=(this->xmax()+this->xmin())/2;
  }
}
//--------------------------------------------------------------------
TMPL void
SELF::subdivide(Cell*& Left, Cell*& Right, int v, double s) {
    cell2d_subdivisor<C,V>::subdivide( *this,(SELF*&)Left,(SELF*&)Right,v,s);
}
//--------------------------------------------------------------------
TMPL  typename SELF::Vector
SELF::gradient(const Point& p) {
    Polynomial
            dx= diff(m_polynomial,0),
            dy= diff(m_polynomial,1);
    double
            u0= (p[0]-this->xmin())/(this->xmax()-this->xmin()),
            u1= (p[1]-this->ymin())/(this->ymax()-this->ymin());
    Vector v;
    v[0] = dx(u0,u1);
    v[1] = dy(u0,u1);
    return v;
}
//====================================================================
} ; // namespace shape
} ; // namespace mmx
# undef Solver
# undef SELF
# undef REF
# undef TMPL
# undef TMPL1
# undef STACK
# endif // shape_cell2d_algebraic_curve_hpp

---