synaps/topology/TopSbdBdg2d.h

/********************************************************************
 *   This file is part of the source code of the SYNAPS kernel.
 *   Author(s): B. Mourrain, GALAAD, INRIA
 ********************************************************************/
#ifndef synaps_topology_TopSbdBdg2d_h_
#define synaps_topology_TopSbdBdg2d_h_

#include <synaps/init.h>
#include <synaps/arithm/gmp.h>
#include <synaps/mpol.h>
#include <synaps/topology/sbd2d.h>
#include <synaps/arithm/Interval.h>

__BEGIN_NAMESPACE_SYNAPS
//====================================================================
template<class O, class X, class I> struct TopSbdBdg2d;
//====================================================================
template<class X, int N>
void print_mx(const MatrDse<Interval<X,N> >& M)
{

  double dmx=0, imx=0,imn=1;
  for(unsigned i=0; i< M.nbrow();i++)
    for(unsigned j=0; j< M.nbcol();j++)
      {
        imx= std::max(imx,M(i,j).sup()-M(i,j).inf());
        imn= std::min(imn,M(i,j).sup()-M(i,j).inf());
        if( M(i,j).sup()<0)
          dmx = std::max(dmx,-M(i,j).sup());
        else
              dmx = std::max(dmx,M(i,j).sup());
          }
      std::cout<<"DMax "<<dmx<<std::endl;
      std::cout<<"IMax "<<imx<<std::endl;
      std::cout<<"IMin "<<imn<<std::endl;
}

template<class X, int N>
void print_sign(const MatrDse<Interval<X,N> >& M)
{
  for(unsigned i=0;i<M.nbrow();i++)
    {
      std::cout<<"[";
      for(unsigned j=0; j<M.nbcol();j++)
      {
        if(M(i,j).sup()<0)
          std::cout<<"-";
        else if(M(i,j).inf()>=0)
          std::cout<<"+";
        else if(M(i,j).inf()<0 && M(i,j).sup()>0)
          std::cout<<"?";
        else
          std::cout<<" ";
        //std::cout<<M(i,j).sup()-M(i,j).inf()<<std::endl;
      }
      std::cout<<"]"<<std::endl;
    }
}
template<class Int>
bool IsSmall (const Int& I)
{
  return (I.inf()<= 0 && I.sup() >= 0 && (I.inf()!=I.sup()));
}

template<class X, int N>
int sign(const Interval<X,N>& I)
{
  if(I.sup()< (0))
    return -1;
  else if(I.inf()>(0))
    return 1;
  else if(I.sup()==0 && I.inf()==0)
    return 0;
  else
    return 2;
}
//--------------------------------------------------------------------
namespace bezier 
{
  template<class X>
  struct bzislv: public Isole<X> 
  {
    X eps;
    bzislv():eps(1e-4){}

    template<class UPOL>
    Seq<X> solve(const UPOL& p, X x0, X x1)
    {
      typedef typename UPOL::value_type coeff_t;
      bezier::rep1d<X> ubp(p.size()), dbp(p.size());
      for(unsigned i=0;i<p.size();i++)
        {
          ubp[i]=p[i].sup();
          dbp[i]=p[i].inf();
          assert(dbp[i]<= ubp[i]);
        }
      this->Run(dbp,ubp,eps,false);
      Seq<X> sol; sol.resize(0);
      this->GetMidPoints(sol);
      for (unsigned i = 0; i < sol.size(); i ++ )
        {
          if ( this->m_res[i].t )
            std::cout<<" xxx"<<std::endl;
        }
      RConvert(sol,x0,x1);             
      return sol;
    }
  };

}//namespace bezier
//====================================================================
template<class O, class X=double, class I= Interval<X,3> >
struct TopSbdBdg2d : public bezier::sbd2d<X,I,bezier::bzislv<X> >
{
  typedef topology::point<X>            point_t;
  typedef topology::point_graph<X>      pointgraph_t;
  typedef I                             interval_t;
  typedef interval_t                    coeff_t;
  typedef cell2d<X,MatrDse<coeff_t> >   cell_t;
  bool draw_box;

  TopSbdBdg2d<O,X,I>(bool bx=true):bezier::sbd2d<X,I,bezier::bzislv<X> >(1.e-3,0.1), draw_box(bx) {}
  TopSbdBdg2d<O,X,I>(const X& e, const X& a, bool bx=true):bezier::sbd2d<X,I,bezier::bzislv<X> >(e,a), draw_box(bx) {}


  int c_sign(const cell_t&);
  int dx_sign(const cell_t&);
  int dy_sign(const cell_t&);
  int regularity(cell_t & c)
  {
      int s0=c_sign(c), s1=4*dx_sign(c);
      s1+=16*dy_sign(c);
      if(s1)
        {
          if(s0==1 || s0 == -1)
            s0+=s1;
          else 
            s0=s1;
        }
      c.reg=s0;
      //      std::cout<<s0<<std::endl;
      return s0;
  }    
  
  template<class mpol_t>
  void assign(cell_t& C, 
              const mpol_t& p,
              const X& x0, const X& x1, const X& y0, const X&y1)
  {
      typedef typename mpol_t::coeff_t Coeff_t;
      MatrDse<Coeff_t> fq;
      BEZIER::monomial2bezier2d(fq,p,Coeff_t(x0),Coeff_t(x1),Coeff_t(y0),Coeff_t(y1));
      Coeff_t mx=0;
      for(unsigned i=0; i< fq.nbrow();i++)
        for(unsigned j=0; j< fq.nbcol();j++)
          if( fq(i,j)<Coeff_t(0))
            mx = std::max(mx,Coeff_t(-fq(i,j)));
          else
            mx = std::max(mx,fq(i,j));
      std::cout<<"Max "<<mx<<std::endl;
      fq/=mx;
      mx=Coeff_t(0);
      C.bz_pol =MatrDse<coeff_t>(fq.nbrow(),fq.nbcol());
      MATRIX::assign(C.bz_pol,fq);
      C.x0=x0; C.x1=x1;C.y0=y0;C.y1=y1;
      //   std::cout<<C.bz_pol<<std::endl;
  }

  template<typename mpol_t>
  void run(topology::point_graph<double>& r, const mpol_t & p, 
           const X& x0, const X& x1, const X& y0, const X& y1);
};

//--------------------------------------------------------------------
template<class O, class X, class I> 
template<class  mpol_t> 
void TopSbdBdg2d<O,X,I>::run(topology::point_graph<double>& r, const mpol_t & p, 
                             const X& x0, const X& x1, const X& y0, const X& y1)
{
    cell_t C;
    assign(C,p,x0,x1,y0,y1);

    double min_sz=size(C)*this->eps,
      max_sz = size(C)*this->asr;
    std::cout<<"Min size "<<min_sz<<std::endl;
    std::cout<<"Max size "<<max_sz<<std::endl;
    if(draw_box) std::cout<<"Drawing singular boxes"<<std::endl;    
    std::vector<cell_t> L, L_sing;
    L.push_back(C);
    while(L.size())
      {
        C=this->pop(L);
        int s=regularity(C);
        //      std::cout<<"reg "<<s<<std::endl;
        if(s==0)
          {
            //print_sign(C.bz_pol);
            if(size(C)<min_sz)
              {
                L_sing.push_back(C);
                if(draw_box)
                  {
                    addbox(r,C);
                    addcenter(r,C);
                  }
              }
            else 
              BEZIER::split(L,C);
          }
        else
          {     
            if(s==2)
              {
                std::cout<<"Zoom [ "<<C.x0<<" "<<C.x1<<" "<<C.y0<<" "<<C.y1<<" ]"<<std::endl;
                assign(C,p,C.x0,C.x1,C.y0,C.y1);
                L.push_back(C);
                if(draw_box)
                  {
                    addbox(r,C);
                    // addcenter(r,C);
                  }
              }
            else if(s%4==0)
              {
                if(size(C)> max_sz)
                  BEZIER::split(L,C);
                else
                  {
                    connect_xy(r,C,(C.reg/4)%4);

                  }
              }
//          else
//            {
//              addbox(r,C);
//              // std::cout <<"Excluded box "<<s<<" "<<(s%4)<<std::endl;
//              if((s%4)==1 || (s%4)==-3)
//                addplus(r,C);
//              else if((s%4)==-1 || (s%4)==3)
//                addmoins(r,C);
//            }
          }
      }
    std::cout <<"Sing cell: "<<L_sing.size()<<std::endl;
    point_t pt(2);
  }
//--------------------------------------------------------------------
template<class O, class X, class I>
int TopSbdBdg2d<O,X,I>::c_sign(const cell_t& M)
{
    typedef typename cell_t::coeff_t coeff_t;
    int s=0;
    bool loop=true;
    unsigned i,j;

    for(i=0;i<M.nbrow() && IsSmall(M(i,0)); i++);
    if(i==M.nbrow()) return 2;
    for(i=0;i<M.nbrow() && IsSmall(M(i,M.nbcol()-1)); i++);
    if(i==M.nbrow()) return 2;

    for(j=0;j<M.nbcol() && IsSmall(M(0,j)); j++);
    if(j==M.nbcol()) return 2;
    for(j=0;j<M.nbcol() && IsSmall(M(M.nbrow()-1,j)); j++);
    if(j==M.nbcol()) return 2;
      
    for(i=0;i<M.nbrow() && loop;i++)
      for(j=0;j<M.nbcol()&& loop;j++)
        if(M(i,j).inf()>(0))
          {
            s= 1;
            loop=false;
          }
        else if(M(i,j).sup()<(0))
          {
            s= -1;
            loop=false;
          }
    for(i=0;i<M.nbrow();i++)
      for(j=0;j<M.nbcol();j++)
        if(s==1 && sign(M(i,j))==-1)
          return 0;
        else if(s==-1 && sign(M(i,j))==1)
          return 0;
    return s;
  }
 
  //--------------------------------------------------------------------
template<class O, class X,class I>
int TopSbdBdg2d<O,X,I>::dx_sign(const cell_t & M)
{
    typedef typename cell_t::coeff_t coeff_t;
    int s=0; bool loop=true;
    for(int i=0;i<M.nbrow() && loop;i++)
      for(int j=0;j<M.nbcol()-1 && loop;j++)
        if(M(i,j+1).inf()>M(i,j).sup())
          {
            s= 1;
            loop=false;
          }
        else if(M(i,j+1).sup()<M(i,j).inf())
          {
            s= -1;
            loop=false;
          }
        else if(IsSmall(interval_t(M(i,j+1)-M(i,j))))
          {
            s=0;
            loop=false;
          }
    for(int i=0;i<M.nbrow();i++)
      for(int j=0;j<M.nbcol()-1;j++)
        {
          interval_t D=M(i,j+1)-M(i,j);
          if(sign(D)!=s)
            return 0;
        }
    return s;
}
//--------------------------------------------------------------------
template<class O,class X, class I>
int TopSbdBdg2d<O,X,I>::dy_sign(const cell_t & M)
{
    typedef typename cell_t::coeff_t coeff_t;
    int s=0; bool loop=true;
    for(unsigned j=0;j<M.nbcol()&& loop;j++)
      for(int i=0;i<M.nbrow()-1&& loop;i++)
        if(M(i+1,j).inf()>M(i,j).sup())
          {
            s= 1;
            loop=false;
          }
        else if(M(i+1,j).sup()<M(i,j).sup())
          {
            s= -1;
            loop=false;
          }
        else if(IsSmall(interval_t(M(i+1,j)-M(i,j))))
          {
            s=0;
            loop=false;
          }
    for(unsigned j=0;j<M.nbcol();j++)
      for(int i=0;i<M.nbrow()-1;i++)
        {
          interval_t D=M(i+1,j)-M(i,j);
          if(sign(D)!= s)
            return 0;
        }
    return s;
}


namespace topology
{
  template<class C, class U>
  void assign(point_graph<C>& r,  const MPol<U>& s, 
              TopSbdBdg2d<C> mth,
              C x0, C x1, C y0, C y1,
              unsigned Z=1)
  {
    typedef QQ coeff_t;
    typedef MPol<coeff_t> mpol_t;
    mpol_t p; let::assign(p,s);
    double eps=mth.eps,a0=x0,b0=x1,a1=y0,b1=y1;
    C d0=(b0-a0)/Z, d1=(b1-a1)/Z;
    for(unsigned i=0;i<Z;i++)
      for(unsigned j=0;j<Z;j++)
        {
          mth.run(r,p,a0+d0*i,a0+(i+1)*d0,a1+d1*j,a1+d1*(j+1));
        }
  }
}
__END_NAMESPACE_SYNAPS 
//====================================================================
#endif //synaps_topology_TopSbdBdg2d_h_