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

Go to the documentation of this file.

/*****************************************************************************
 * M a t h e m a  g i x
 *****************************************************************************
 * solver_implicit
 * 2008-11-09
 * Bernard Mourrain
 *****************************************************************************
 *               Copyright (C) 2008 INRIA Sophia-Antipolis
 *****************************************************************************
 * Comments :
 ****************************************************************************/
#ifndef shape_implicit_solver_hpp
#define shape_implicit_solver_hpp
//===========================================================================
# include <realroot/Interval.hpp>
# include <realroot/polynomial_bernstein.hpp>
# include <realroot/GMP.hpp>
# include <realroot/solver_uv_sleeve.hpp>
# include <realroot/solver_mv_bernstein.hpp>
# include <shape/point.hpp>
# include <shape/bounding_box.hpp>
# include <shape/box_face.hpp>
# include <shape/topological_degree.hpp>

# define TMPL template<class C, class V>
# define SELF solver_implicit<C,V>

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

inline double lower(double x){return x;}
inline double upper(double x){return x;}

template<class S1, class S2>
point<double>* point_on_edge(const bounding_box<double>& bx, double u, const S1& s1, const S2& s2)  
{
  typedef point<double> Point;
  double p[3];
  p[S1::var] = S1::eval(bx);
  p[S2::var] = S2::eval(bx);
  int v= 3-S1::var-S2::var;
  p[v]= lower(bx,v)*(1-u) + upper(bx,v)*u;
  
  return new Point(p[0],p[1],p[2]);
  //    delete p;
}
//--------------------------------------------------------------------
TMPL struct solver_implicit {

  typedef REF_OF(V) Ref;
  typedef vertex<C,V> Point; 
  typedef bounding_box<C,Ref>     BoundingBox;
  typedef box_face<C,Ref>           BoxFace;
    
  static  BoxFace north_edge;
  static  BoxFace south_edge;
  static  BoxFace west_edge;
  static  BoxFace east_edge;
  static  BoxFace front_edge;
  static  BoxFace back_edge;

  static  BoxFace north_back_edge;
  static  BoxFace south_back_edge;
  static  BoxFace west_back_edge;
  static  BoxFace east_back_edge;
  
  static  BoxFace north_front_edge;
  static  BoxFace south_front_edge;
  static  BoxFace west_front_edge;
  static  BoxFace east_front_edge;
  
  static  BoxFace north_west_edge;
  static  BoxFace south_west_edge;
  static  BoxFace north_east_edge;
  static  BoxFace south_east_edge;

  static  BoxFace north_face;
  static  BoxFace south_face;
  static  BoxFace west_face;
  static  BoxFace east_face;
  static  BoxFace front_face;
  static  BoxFace back_face;
  
  static inline Point* new_point(double a, double b, double c) {
    return new Point(a,b,c);
  }
  
  //--------------------------------------------------------------------
  template<class MultivariateDenseBernstein> static
  void edge_point(Seq<Point *>& lp, 
                  const MultivariateDenseBernstein & p, 
                  const BoxFace & E,  
                  const BoundingBox & bx) {
      MultivariateDenseBernstein f;

      tensor::face(f,p,E.cvar(0),E.side(0));

      polynomial<double, with<Bernstein> > dw(1,f.size()-1), up(1,f.size()-1);
      for(unsigned i=0; i<f.size();i++) {
        up[i]= upper(f[i]); 
        dw[i]= lower(f[i]);
      }

      double 
        U = E.lower(bx),
        V = E.upper(bx);

      Seq<double> sol;
      solver< double, Sleeve<Approximate> >::solve(sol,dw,up,U,V);

      // std::cout<<"****U= "<<U<<std::endl;
      // std::cout<<"****V= "<<V<<std::endl;
      // std::cout<<"***dw= "<<dw<<std::endl;
      // std::cout<<"***up= "<<up<<std::endl;
      // std::cout<<"******Sols : "<<sol.size()<<std::endl;

      Approximate approx;
      std::vector<double> pt(3);
      for(unsigned i=0;i<sol.size(); i++) 
        if(E.is_valid(bx, sol[i], approx.eps)) {
          lp<< E.new_point(bx,sol[i]);
        }
    }
  //--------------------------------------------------------------------
  template<class MultivariateDenseBernstein> static
  int edge_sign_var(const MultivariateDenseBernstein & p, 
                    const BoxFace & E) {
      MultivariateDenseBernstein f;
      tensor::face(f,p,E.cvar(0),E.side(0));
      return sign_variation(f.begin(), f.end());

  }
  //--------------------------------------------------------------------
  template<class MultivariateDenseBernstein> static
  void face_point(Seq<Point *>& sol, 
                  const Seq<MultivariateDenseBernstein>& eqs, 
                  const BoxFace& F,
                  const BoundingBox& bx) {
      
      typedef polynomial< double, with<Bernstein> > Polynomial;

      Polynomial f, p;
      Seq<Polynomial> sys;
      for(unsigned i=0;i< eqs.size(); i++) {
        let::assign(p.rep(), eqs[i].rep());
        tensor::face(f.rep(), p.rep(), F.cvar(0), F.side(0));
        sys<<f;
      }

      solver<double, ProjRd<SBD_RDRDL> > slv;
      std::vector<double> res; 
      slv.solve_bernstein(res,sys);
      Approximate approx;

      for(int i=0;i<(int)res.size()-1;i+=2) {

//      if (res[i]<approx.eps) res[i]=0;
//      else if (res[i]>1-approx.eps) res[i]=1;

//      if (res[i+1]<approx.eps) res[i+1]=0;
//      else if (res[i+1]>1-approx.eps) res[i+1]=1;

        Point * p = F.new_point(bx, res[i], res[i+1]);
        if (F.is_valid(*p, bx))  sol<< p;
        //      else delete p;
      }
    }
  //--------------------------------------------------------------------
  template<class Polynomial>
  static void intersection(Seq<Point* >& sol, 
                           const Polynomial& f1, 
                           const Polynomial& f2, const BoundingBox & bx) {

      typedef polynomial< double, with<Bernstein> > MultivariateDenseBernstein;
      MultivariateDenseBernstein p1, p2;

      
      let::assign(p1.rep(), f1.rep() );
      let::assign(p2.rep(), f2.rep() );
      
      Seq<MultivariateDenseBernstein> sys;
      sys<<p1<<p2;

      solver<double, ProjRd<SBD_RDRDL> > slv;
      std::vector<double> res; 
      slv.solve_bernstein(res,sys);

      double x_sc= (bx.xmax()-bx.xmin()), y_sc= (bx.ymax()-bx.ymin());
      for(unsigned i=0;i<res.size();i+=2)
        {
          sol << new_point(bx.xmin()+res[i]*x_sc, bx.ymin()+res[i+1]*y_sc, 0.0);
        }
    }
  //--------------------------------------------------------------------
  template<class Polynomial>
  static void extremal(Seq<Point* >& sol, 
                       const Polynomial& pol, const BoundingBox& b) {
    
      typedef polynomial< double, with<Bernstein> >  MultivariateDenseBernstein;

      // typedef Polynomial MultivariateDenseBernstein;


      typedef polynomial< GMP::rational, with<Bernstein> > MultivariateDenseBernsteinQ;
      
      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> P(pol.rep(),bx);
      //typename Polynomial::Scalar mx = as<typename Polynomial::Scalar>(array::max_abs(P));
      //Polynomial Pol=pol/mx, dxpol = diff(Pol,0), dypol = diff(Pol,1);

      Polynomial dxpol = diff(pol,0), dypol = diff(pol,1);
      tensor::bernstein<mmx::GMP::rational>
         dxP(dxpol.rep(),bx),dyP(dypol.rep(),bx);
      
      MultivariateDenseBernstein dxp,dyp;
      let::assign(dxp.rep(), dxP);
      let::assign(dyp.rep(), dyP);

      Seq<MultivariateDenseBernstein> sys;
      //      Seq<MultivariateDenseBernsteinQ> sys;
      sys<<dxp<<dyp;

      solver<double, ProjRd<SBD_RDRDL> > slv;

      std::vector<double> res; 
      slv.solve_bernstein(res,sys);

      double x_sc= (b.xmax()-b.xmin()), y_sc= (b.ymax()-b.ymin());
      for(unsigned i=0;i<res.size();i+=2) {
        sol << new_point(b.xmin()+res[i]*x_sc, b.ymin()+res[i+1]*y_sc, 0.0);
      }
    }
  //--------------------------------------------------------------------
  template<class Polynomial>
  static void singular(Seq<Point* >& sol, 
                       const Polynomial& pol, const BoundingBox& b) {

      typedef polynomial< double, with<Bernstein> >        MultivariateDenseBernstein;
      typedef polynomial< GMP::rational, with<Bernstein> > MultivariateDenseBernsteinQ;
      
      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());

      Polynomial dxpol = diff(pol,0), dypol = diff(pol,1);
      tensor::bernstein<mmx::GMP::rational>
         P(pol.rep(),bx), dxP(dxpol.rep(),bx),dyP(dypol.rep(),bx);
      
      MultivariateDenseBernstein p,dxp,dyp;
      let::assign(p.rep()  , P);
      let::assign(dxp.rep(), dxP);
      let::assign(dyp.rep(), dyP);

      Seq<MultivariateDenseBernstein> sys;
      sys<<p<<dxp<<dyp;

      solver<double, ProjRd<SBD_RDRDL> > slv;
      //      solver<ring<double,Bernstein>, ProjRd<SBD_RDRDL> >::eps=1e-10;

      std::vector<double> res; 
      slv.solve_bernstein(res,sys);

      double x_sc= (b.xmax()-b.xmin()), y_sc= (b.ymax()-b.ymin());
      for(unsigned i=0;i<res.size();i+=2) {
        sol << new_point(b.xmin()+res[i]*x_sc, b.ymin()+res[i+1]*y_sc, 0.0);
      }
    }

  //--------------------------------------------------------------------
  template<class MultivariateDenseBernstein> static
  void common_edge_point(Seq<Point *>& lp, 
                         const MultivariateDenseBernstein & p, 
                         const BoxFace & E,  
                         const BoundingBox & bx1, 
                         const BoundingBox & bx2 ) {
      MultivariateDenseBernstein f;

      tensor::face(f,p,E.cvar(0),E.side(0));

      polynomial<double, with<Bernstein> > dw(1,f.size()-1), up(1,f.size()-1);
      for(unsigned i=0; i<f.size();i++) {
        up[i]= upper(f[i]); 
        dw[i]= lower(f[i]);
           //up[i]=dw[i]=f[i];
      }

      double 
        U = std::max( E.lower(bx1), E.lower(bx2) ) ,
        V = std::min( E.upper(bx1), E.upper(bx2) ) ;

      Seq<double> sol;

      solver< double, Sleeve<Approximate> >::solve(sol,dw,up,U,V);

      sol.sort();// sort from smaller to bigger

      Approximate approx;
      for(unsigned i=0;i<sol.size(); i++) 
          {
            lp<< E.new_point(bx1,sol[i]);
          }  
    }
};
//====================================================================
# define REF REF_OF(V)
TMPL box_face<C,REF> SELF::north_face(1,1);
TMPL box_face<C,REF> SELF::south_face(1,0);
TMPL box_face<C,REF> SELF::west_face (0,0);
TMPL box_face<C,REF> SELF::east_face (0,1);
TMPL box_face<C,REF> SELF::front_face(2,1);
TMPL box_face<C,REF> SELF::back_face(2,0);

// 2D box  edges;
TMPL box_face<C,REF> SELF::north_edge(1,1,2,0);
TMPL box_face<C,REF> SELF::south_edge(1,0,2,0);
TMPL box_face<C,REF> SELF::west_edge(0,0,2,0);
TMPL box_face<C,REF> SELF::east_edge(0,1,2,0);
  
TMPL box_face<C,REF> SELF::front_edge(2,1,0,0);
TMPL box_face<C,REF> SELF::back_edge (2,0,0,0);
  
  // 3D box  edges;
TMPL box_face<C,REF> SELF::north_back_edge (1,1,2,0);
TMPL box_face<C,REF> SELF::south_back_edge (1,0,2,0);
TMPL box_face<C,REF> SELF::west_back_edge  (0,0,2,0);
TMPL box_face<C,REF> SELF::east_back_edge  (0,1,2,0);
  
TMPL box_face<C,REF> SELF::north_front_edge(1,1,2,1);
TMPL box_face<C,REF> SELF::south_front_edge(1,0,2,1);
TMPL box_face<C,REF> SELF::west_front_edge (0,0,2,1);
TMPL box_face<C,REF> SELF::east_front_edge (0,1,2,1);
  
TMPL box_face<C,REF> SELF::north_west_edge (0,0,1,1);
TMPL box_face<C,REF> SELF::south_west_edge (0,0,1,0);
TMPL box_face<C,REF> SELF::north_east_edge (0,1,1,1);
TMPL box_face<C,REF> SELF::south_east_edge (0,1,1,0);
  
//===========================================================================
} //namespace shape
} //namespace mmx

# undef TMPL
# undef SELF
# undef REF
# undef BoxFace
# endif //shape_implicit_solver