synaps/topology/sbd2d.h

/********************************************************************
 *   This file is part of the source code of the SYNAPS kernel.
 *   Author(s): B. Mourrain, GALAAD, INRIA
 *   $Id: topology2d.h,v 1.2 2005/09/28 06:40:47 mourrain Exp $
 ********************************************************************/
#ifndef synaps_topology_sbd2d_H
#define synaps_topology_sbd2d_H
//--------------------------------------------------------------------
#include <synaps/init.h>
#include <synaps/arithm/sign.h>
#include <synaps/upol/UPolDse.h>
#include <synaps/upol/bezier/BEZIER.m>
#include <synaps/linalg/VectDse.h>
#include <synaps/linalg/MatrDse.h>
#include <synaps/upol/bezier/rep1d.h>
#include <synaps/usolve/bezier/sbd1d.h>
#include <synaps/arithm/sign.h>
#include <synaps/topology/cell2d.h>
#include <synaps/topology/point_graph.h>
//--------------------------------------------------------------------
__BEGIN_NAMESPACE_SYNAPS
//--------------------------------------------------------------------
struct YX {
  template<class T>
  static bool less(const T& a, const T& b)
  {
    return (a[1]<b[1] || (a[1]==b[1] && a[0]< b[0]));
  }
  template<class T>
  bool operator()(const T& a, const T& b)
  {
    return less(a,b);
  }
};


struct XY {
  template<class T>
  static bool less(const T& a, const T& b)
  {
    return (a[0]<b[0] || (a[0]==b[0] && a[1]< b[1]));
  }
  template<class T>
  bool operator()(const T& a, const T& b)
  {
    return less(a,b);
  }
};
//====================================================================
namespace BEZIER {

//====================================================================
template<class MATR, class MPOL, class X>
void monomial2bezier2d(MATR& M, const MPOL& p, 
                       const X& x0, const X& x1, const X& y0, const X& y1)
{
  using let::assign;
  typedef typename MATR::coeff_t coeff_t;
  M.resize(degree(p,1)+1, degree(p,0)+1);
  for(typename  MPOL::const_iterator it=p.begin(); it != p.end(); it++)
    {
      assign(M((*it)[1],(*it)[0]),it->coeff());
    }
  for(unsigned i=0;i<M.nbrow();i++)
    {
      VectDse<coeff_t> V = row(M,i), BZ(M.nbcol(),AsSize());
      BEZIER::monomials2bezier(BZ,V,V.size(),x0,x1);
      //      row(M,i)=BZ;
      for(unsigned j=0;j<M.nbcol();j++) M(i,j)=BZ[j];
    }

  for(unsigned j=0;j<M.nbcol();j++)
    {
      VectDse<coeff_t> V= col(M,j), BZ(M.nbrow(),AsSize());
      BEZIER::monomials2bezier(BZ,V,V.size(),y0,y1);
      //      col(M,j)=BZ;
      for(unsigned i=0;i<M.nbrow();i++) M(i,j)=BZ[i];
    }
}
//--------------------------------------------------------------------
template<class matr_t>
void deCasteljau( matr_t& left, matr_t& right, ByRow M ) 
{
    typedef typename matr_t::value_type coeff_t;

    int col=left.nbcol(),row=left.nbrow();
    assert( col > 1 );
    coeff_t middle = coeff_t(0.5);
    for ( int i = 0 ; i < row ; ++i )
      {
        for ( int j = 1 ; j < col ; ++j )
          {
            left( i, j-1) = right( i, 0); 
            for ( int k = 0 ; k < col - j ; ++k ) 
              {
                right( i, k ) = ( right( i, k ) + right( i, k+1) ) * middle;
              }
          }
        left( i, col-1 ) = right( i, 0);
      }
}  
//--------------------------------------------------------------------
template<class matr_t>
void deCasteljau( matr_t & left, matr_t & right, ByCol M ) 
{

  typedef typename matr_t::value_type coeff_t;

  int col=right.nbcol(),row=right.nbrow();
  assert( row > 1 );
  coeff_t middle = coeff_t(0.5);
  for ( int i = 0 ; i < col ; ++i )
    {
      for ( int j = 1 ; j < row ; ++j )
        {
          left( j-1, i) = right( 0, i); 
          for ( int k = 0 ; k < row - j ; ++k ) 
            {
              right( k, i ) = ( right( k, i ) + right( k+1, i ) ) * middle;
            }
        }
      left( row-1, i) = right( 0, i);
    }
}  

  //--------------------------------------------------------------------
  template<class cell_t> 
  void split(std::vector<cell_t >& L, const cell_t& c)
  {
    typedef typename cell_t::matr_t matr_t;

    //std::cout<<"Split "<<" ["<<c.x0<<" "<<c.x1<<"; "<<c.y0<<" "<<c.y1<<"]"
    //         <<std::endl;
    //std::cout<<c.bz_pol<<std::endl;
    matr_t M1(c.bz_pol),M2(c.bz_pol);
    deCasteljau(M1,M2,ByRow());
    matr_t M11(M1),M21(M1);
    deCasteljau(M11,M21,ByCol());
    //  std::cout<<M11<<std::endl;
    //  std::cout<<M21<<std::endl;
    L.push_back(cell_t(c.x0,(c.x0+c.x1)/2,c.y0,(c.y0+c.y1)/2,M11));
    L.push_back(cell_t(c.x0,(c.x0+c.x1)/2,(c.y0+c.y1)/2,c.y1,M21));
    matr_t M12(M2),M22(M2);
    deCasteljau(M12,M22,ByCol());
    //  std::cout<<M12<<std::endl;
    //  std::cout<<M22<<std::endl;
    L.push_back(cell_t((c.x0+c.x1)/2,c.x1,c.y0,(c.y0+c.y1)/2,M12));
    L.push_back(cell_t((c.x0+c.x1)/2,c.x1,(c.y0+c.y1)/2,c.y1,M22));
  }
}  
//====================================================================
namespace bezier 
{
  template<class X>
  struct bzslv: public Isole<X> 
  {

    X eps;
    bzslv():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<coeff_t> bz(p.size());
      VECTOR::assign(bz,p);
      Seq<X> sol;
      this->Run(bz,eps*(x1-x0),false);
      this->GetMidPoints(sol);
      RConvert( sol, x0,x1);
      return sol;
    }  
  };

  template<class X, class C=X, class USLV= bezier::bzslv<C> >
  struct sbd2d  
  {
    typedef Seq<X>                  u_solution;
    typedef topology::point<X>        point_t;
    typedef topology::point_graph<X>  pointgraph_t;
    typedef C                       coeff_t;
    typedef cell2d<X,MatrDse<C> >   cell_t;

    USLV   slv;
    X eps, asr;
    sbd2d():eps(1e-3), asr(1) {}
    sbd2d(X e):eps(e), asr(1) {}
    sbd2d(X e, X a):eps(e), asr(a) {}
    
    cell_t pop(std::vector<cell_t>& L)
    {
      cell_t c=L[L.size()-1];
      L.resize(L.size()-1);
      return c;
    }

    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;
      monomial2bezier2d(fq,p,x0,x1,y0,y1);
      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;
    }

    void connect_xy(pointgraph_t & g, const cell_t & c, bool b);

    int c_sign(const cell_t&);
    int dx_sign(const cell_t&);
    int dy_sign(const cell_t&);
    int regularity(cell_t & c);
   
    template<typename mpol_t>
    void run(pointgraph_t& g, const mpol_t & p, 
             const X& x0, const X& x1, const X& y0, const X& y1);
  };
 
  //--------------------------------------------------------------------
  template<class X, class U, class S>
  void sbd2d<X,U,S>::connect_xy(pointgraph_t & g, const cell_t & c, bool b)
  {
    typedef typename pointgraph_t::point_t point_t;
    point_t pt(2);
    int n0, n1;
    Seq<point_t> l;
    
    Seq<X> sol0 =slv.solve((typename cell_t::row_type)row(c.bz_pol,0),c.x0,c.x1);
    pt[1]=c.y0;
    for(unsigned i=0;i<sol0.size();i++)
      {
        pt[0]=sol0[i]; 
        l.push_back(pt);
      }
    Seq<X> sol1 =slv.solve((typename cell_t::row_type)row(c.bz_pol,c.bz_pol.nbrow()-1),
                           c.x0,c.x1);
    pt[1]=c.y1;
    for(unsigned i=0;i<sol1.size();i++)
      {
        pt[0]=sol1[i]; l.push_back(pt);
    }
    Seq<X> sol2 =slv.solve((typename cell_t::col_type)col(c.bz_pol,0),c.y0,c.y1);
    pt[0]=c.x0;
    for(unsigned i=0;i<sol2.size();i++)
      {
        pt[1]=sol2[i]; l.push_back(pt);
      }
    Seq<X> sol3 =slv.solve((typename cell_t::col_type)col(c.bz_pol,c.bz_pol.nbcol()-1),
                           c.y0,c.y1);
    pt[0]=c.x1;
    for(unsigned i=0;i<sol3.size();i++)
      {
        pt[1]=sol3[i]; l.push_back(pt);
      }
    if(b)
      sort(l.begin(),l.end(), YX());
    else
      sort(l.begin(),l.end(), XY());
    if(l.size()%2 !=0) 
      {
        std::cout<<c.x0<<" "<<c.x1<<" "<<c.y0<<" "<<c.y1<<"\n  ";
        VECTOR::print(std::cout,l)<<" => ";
        Seq<point_t> ll;
        for(unsigned i=0;i<l.size()-1;i++)
          if(fabs(l[i][0]-l[i+1][0])> slv.eps || 
             fabs(l[i][1]-l[i+1][1])> slv.eps)
            ll.push_back(l[i]);
        ll.push_back(l[l.size()-1]);
        l=ll;
        VECTOR::print(std::cout,l)<<std::endl;
      }
    for(int i=0;i<l.size(); i+=2)
      {
        if(i<l.size()-1)
        {
          n0=g.insert(l[i]);
          n1=g.insert(l[i+1]);
          g.insert(n0,n1);
        }
      }
  }
  
  //--------------------------------------------------------------------
  template<class X, class U, class S> 
  template<class  mpol_t> 
  void sbd2d<X,U,S>::run(topology::point_graph<X>& r, const mpol_t & p, 
                         const X& x0, const X& x1, const X& y0, const X& y1)
  {
    typedef typename mpol_t::coeff_t       CF_T;
    typedef cell2d<coeff_t>                cell_t;
   
   
    cell_t C;
    assign(C,p,x0,x1,y0,y1);

    double sz = size(C)*asr;
    std::vector<cell_t> L, L_sing;
    L.push_back(C);
    while(L.size())
      {
        C=pop(L);
        int s=regularity(C);
        if(s==0)
          {
            if(size(C)<eps)
              {
                L_sing.push_back(C);
                addbox(r,C);
                addcenter(r,C);
              }
            else 
              BEZIER::split(L,C);
          }
        else
          {     
            if(s%4==0)
              {
                if(size(C)> sz)
                  split(L,C);
                else
                  {
                    connect_xy(r,C, (C.reg/4) %4);
                  }
              }
          }
      }
    //    std::cout <<"Reg. cell: "<<L_reg.size()<<std::endl;
    std::cout <<"Sing cell: "<<L_sing.size()<<std::endl;
  }
  
  //--------------------------------------------------------------------
  template<class X, class U, class S>
  int sbd2d<X,U,S>::c_sign(const cell_t & M)
  {
    typedef typename cell_t::coeff_t coeff_t;
    int s=0;
    bool run=true;
    for(unsigned i=0;i<M.nbrow() && run;i++)
      for(unsigned j=0;j<M.nbcol()&& run;j++)
        if(M(i,j)>(0))
          {
            s= 1;
            run=false;
          }
        else if(M(i,j)<(0))
          {
            s= -1;
            run=false;
          }
    // std::cout<<"s "<<s<<std::endl;
    
    for(unsigned i=0;i<M.nbrow();i++)
      for(unsigned j=0;j<M.nbcol();j++)
        if(sign(M(i,j))==-s)
          return 0;
    return s;
  }
  //--------------------------------------------------------------------
  template<class X, class U, class S>
  int sbd2d<X,U,S>::dx_sign(const cell_t & M)
  {
  typedef typename cell_t::coeff_t coeff_t;
  int s=0; bool run=true;
  for(int i=0;i<M.nbrow() && run;i++)
    for(int j=0;j<M.nbcol()-1 && run;j++)
      if(M(i,j+1)>M(i,j))
        {
          s= 1;
          run=false;
        }
      else if(M(i,j+1)<M(i,j))
        {
          s= -1;
          run=false;
        }
  for(int i=0;i<M.nbrow();i++)
    for(int j=0;j<M.nbcol()-1;j++)
      {
        if(sign(M(i,j+1)-M(i,j))==-s)
          return 0;
      }
  return s;
}
//--------------------------------------------------------------------
  template<class X, class U, class S>
  int sbd2d<X,U,S>::dy_sign(const cell_t & M)
  {
    typedef typename cell_t::coeff_t coeff_t;
    int s=0; bool run=true;
    for(unsigned j=0;j<M.nbcol()&& run;j++)
      for(int i=0;i<M.nbrow()-1&& run;i++)
        if(M(i+1,j)>M(i,j))
          {
            s= 1;
            run=false;
          }
        else if(M(i+1,j)<M(i,j))
          {
            s= -1;
            run=false;
          }
    
    for(unsigned j=0;j<M.nbcol();j++)
      for(int i=0;i<M.nbrow()-1;i++)
        {
        if(sign(M(i+1,j)-M(i,j))!= s)
          return 0;
      }
    return s;
  }
  //--------------------------------------------------------------------
  template<class X,class U,class S>
  int sbd2d<X,U,S>::regularity(cell_t& c)
  {
    int s0= c_sign(c);
    s0 += 4*dx_sign(c);
    s0 += 16*dy_sign(c);
    c.reg=s0;
    return s0;
  }    
  //--------------------------------------------------------------------  
}//namespace bezier
//====================================================================
template<class cell_t, class X>
void addbox(topology::point_graph<X>& g, const cell_t & c)
{
  typedef typename topology::point_graph<X> pointgraph_t;
  typedef typename pointgraph_t::point_t point_t;
  point_t p1(2), p2(2);
  unsigned n1, n2;

  X a0=c.x0,a1=c.x1,b0=c.y0,b1=c.y1;

  p1[0]=a0;p1[1]=b0;
  p2[0]=a0;p2[1]=b1;
  n1=g.insert(p1);
  n2=g.insert(p2);
  g.insert(n1,n2);

  p2[0]=a1;p2[1]=b0;
  n2=g.insert(p2);
  g.insert(n1,n2);

  p1[0]=a0;p1[1]=b1;
  p2[0]=a1;p2[1]=b1;
  n1=g.insert(p1);
  n2=g.insert(p2);
  g.insert(n1,n2);

  p1[0]=a1;p2[1]=b0;
  n1=g.insert(p2);
  g.insert(n1,n2);
}
//--------------------------------------------------------------------
template<class cell_t, class X>
void addplus(topology::point_graph<X>& g, const cell_t & c)
{
  topology::point<X> p1(2),p2(2);
  p1[0]=(c.x0+c.x1)/2;
  p1[1]=(c.y0+3*c.y1)/4;
  int n1 = g.insert(p1);
  p2[0]=(c.x0+c.x1)/2;
  p2[1]=(3*c.y0+c.y1)/4;
  int n2=  g.insert(p2);
  g.insert(n1,n2);

  p1[0]=(3*c.x0+c.x1)/4;
  p1[1]=(c.y0+c.y1)/2;
  n1 = g.insert(p1);
  p2[0]=(c.x0+3*c.x1)/4;
  p2[1]=(c.y0+c.y1)/2;
  n2=  g.insert(p2);
  g.insert(n1,n2);

}       

template<class cell_t, class X>
void addmoins(topology::point_graph<X>& g, const cell_t & c)
{
  topology::point<X> p1(2),p2(2);
  p1[0]=(3*c.x0+c.x1)/4;
  p1[1]=(c.y0+c.y1)/2;
  int n1 = g.insert(p1);
  p2[0]=(c.x0+3*c.x1)/4;
  p2[1]=(c.y0+c.y1)/2;
  int n2=  g.insert(p2);
  g.insert(n1,n2);

}       

template<class cell_t, class X>
void addcenter(topology::point_graph<X>& g, const cell_t & c)
{
  topology::point<X> pt(2);
  pt[0]=(c.x0+c.x1)/2;
  pt[1]=(c.y0+c.y1)/2;
  g.insert(pt);
}       
//--------------------------------------------------------------------
__END_NAMESPACE_SYNAPS
//====================================================================
#endif //synaps_topology_sbd2d_H