include/shape/cell3d_algebraic_surface.hpp

Go to the documentation of this file.

/*****************************************************************************
 * M a t h e m a  g i x
 *****************************************************************************
 * bcell of algebraic surface
 * 2008-03-20
 * Bernard Mourrain & Julien Wintz
 *****************************************************************************
 *               Copyright (C) 2008 INRIA Sophia-Antipolis
 *****************************************************************************
 * Comments :
 ****************************************************************************/

# ifndef shape_cell3d_algebraic_surface_hpp
# define shape_cell3d_algebraic_surface_hpp

# include <realroot/Interval.hpp>
# include <realroot/ring_bernstein_tensor.hpp>
# include <realroot/sign_variation.hpp>
# include <realroot/Seq.hpp>
# include <shape/bounding_box.hpp>
# include <shape/algebraic_surface.hpp>
# include <shape/cell3d.hpp>
# include <shape/solver_implicit.hpp>
# include <shape/marching_cube.hpp>

# define TMPL template<class C, class V>
# define TMPL1 template<class C>
# define SELF cell3d_algebraic_surface<C,V>
# undef Topology
//====================================================================
namespace mmx { namespace shape {
//====================================================================

//TMPL struct topology;
//TMPL struct tpl3d;
//TMPL struct cell3d_algebraic_surface;

//--------------------------------------------------------------------
struct cell3d_algebraic_surface_def {};

template<> struct use<cell3d_algebraic_surface_def>
        //   :public use<topology_def>
        //   ,public use<algebraic_surface_def>
{
    typedef Interval< double >                    Scalar;
    typedef polynomial< double, with<Bernstein> > Polynomial;
    //  typedef use<cell3d_def>::Cell3d           Cell3d;
    //  typedef use<cell3d_def>::BoundingBox      BoundingBox;
    //  typedef cell3d_algebraic_surface<double,double> Cell3dAlgebraicSurface;
}; 

//--------------------------------------------------------------------
template<class C, class V=default_env>
struct cell3d_algebraic_surface : public cell3d<C,V> {
public:

    typedef REF_OF(V)  Ref;
    typedef tpl3d<C,V>                       Topology3d;
    typedef topology<C,V>                    Topology;

    typedef typename use<mesh3d_def,C,V>::Point         Point;
    typedef typename Topology::Edge          Edge;
    typedef typename Topology::Face          Face;

    typedef cell3d<C,V>                      Cell;
    typedef cell3d<C,V>                      CellBase;
    typedef cell3d<C,V>                      Cell3d;

    typedef bounding_box<C,V>              BoundingBox;

    typedef algebraic_surface<C,V>         AlgebraicSurface;

    typedef typename use<cell3d_algebraic_surface_def,C,V>::Polynomial  Polynomial;

    //  typedef polynomial< Interval<C>, with<Bernstein> > Polynomial;
    //  typedef polynomial< Scalar, with<Bernstein> > Polynomial;

    typedef solver_implicit<C,V>             Solver;

    cell3d_algebraic_surface(const SELF&);
    cell3d_algebraic_surface(const Polynomial&, const BoundingBox& );
    cell3d_algebraic_surface(char*, const BoundingBox& );

    cell3d_algebraic_surface(const AlgebraicSurface&, const BoundingBox&);
    cell3d_algebraic_surface(AlgebraicSurface *, const BoundingBox&);

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

    template<class CELL>
    void split    (CELL*& left, CELL*& right, int v, double s);
    virtual void subdivide(Cell*& left, Cell*& right, int v, double s);

    int          mc_index(void)      const ;
    bool         edge_point(Point**, int ) const ;

    C            value(const C& x, const C& y, const C& z) const;
    C            center_value() const { return m_polynomial(0.5,0.5,0.5); }
    C            lower_bound() const  {return min_value<C>(m_polynomial.begin(), m_polynomial.end());}
    C            upper_bound() const  {return max_value<C>(m_polynomial.begin(), m_polynomial.end());}

    double       vertex_eval(unsigned sx, unsigned sy, unsigned sz) const;
    const Polynomial&   equation() const {return m_polynomial;}
    Polynomial   get_polynomial() const {return m_polynomial;}

    void         point(Point& r, const Point& A, const Point& B) const;
public:
    Polynomial        m_polynomial;
    Point*            m_center;
    C                 m_center_value;
    int               m_idx;


};

//====================================================================
TMPL
SELF::cell3d_algebraic_surface(const SELF& s)
    : Cell3d((BoundingBox)s), m_polynomial(s.equation()), m_center(0)
{

}
//--------------------------------------------------------------------
TMPL
SELF::cell3d_algebraic_surface(const Polynomial& pol, const BoundingBox& bx)
    : Cell3d(bx), m_polynomial(pol), m_center(0)
{

}
//--------------------------------------------------------------------
TMPL
SELF::cell3d_algebraic_surface(char* pol, const BoundingBox& bx)
    : Cell3d(bx), m_polynomial(pol), m_center(0)
{
}
//--------------------------------------------------------------------
TMPL
SELF::cell3d_algebraic_surface(const AlgebraicSurface& sf, const BoundingBox& b): Cell3d(b), m_center(0)
{
    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;
    tensor::bernstein<mmx::GMP::rational> polq(sf.equation().rep(),bx);
    let::assign(tmp.rep(),polq);
    m_polynomial=tmp;

}
//--------------------------------------------------------------------
TMPL bool 
SELF::is_active()  const {
    return has_sign_variation(m_polynomial.begin(),m_polynomial.end());
}
//--------------------------------------------------------------------
TMPL bool 
SELF::is_regular() { 
    Polynomial dxf, dyf, dzf;
    Polynomial fp0, fp1;
    Polynomial bk, ft;


    tensor::face(bk, m_polynomial, 2, 0);
    if(Solver::edge_sign_var(bk, Solver::north_back_edge) >1) return false;
    if(Solver::edge_sign_var(bk, Solver::south_back_edge) >1) return false;
    if(Solver::edge_sign_var(bk, Solver::east_back_edge ) >1) return false;
    if(Solver::edge_sign_var(bk, Solver::west_back_edge ) >1) return false;


    tensor::face(ft, m_polynomial, 2, 1);
    if(Solver::edge_sign_var(ft, Solver::north_front_edge) >1) return false;
    if(Solver::edge_sign_var(ft, Solver::south_front_edge) >1) return false;
    if(Solver::edge_sign_var(ft, Solver::east_front_edge ) >1) return false;
    if(Solver::edge_sign_var(ft, Solver::west_front_edge ) >1) return false;

    tensor::face(ft, equation(), 1, 1);
    if(Solver::edge_sign_var(ft, Solver::north_west_edge) >1) return false;
    if(Solver::edge_sign_var(ft, Solver::north_east_edge) >1) return false;

    tensor::face(bk, equation(), 1, 0);
    if(Solver::edge_sign_var(bk, Solver::south_west_edge) >1) return false;
    if(Solver::edge_sign_var(bk, Solver::south_east_edge) >1) return false;

    tensor::diff(dxf.rep(),m_polynomial.rep(),0);
    tensor::diff(dyf.rep(),m_polynomial.rep(),1);
    tensor::diff(dzf.rep(),m_polynomial.rep(),2);

    if(!has_sign_variation(dxf.begin(),dxf.end())) {

        this->m_idx=0;
        return true;

        tensor::face(fp0, dyf, 0, 0);
        tensor::face(fp1, dyf, 0, 1);
        if(!has_sign_variation(fp0.begin(),fp0.end()) &&
                !has_sign_variation(fp1.begin(),fp1.end()) ) return true;

        tensor::face(fp0, dzf, 0, 0);
        tensor::face(fp1, dzf, 0, 1);
        if(!has_sign_variation(fp0.begin(),fp0.end()) &&
                !has_sign_variation(fp1.begin(),fp1.end())) return true;

        return false;
    }

    //tensor::diff(dyf.rep(),m_polynomial.rep(),1);
    if(!has_sign_variation(dyf.begin(),dyf.end())) {
        this->m_idx=1;
        return true;
    }

    //tensor::diff(dzf.rep(),m_polynomial.rep(),2);
    if(!has_sign_variation(dzf.begin(),dzf.end())) {
        this->m_idx=2;
        return true;
    }

    return false;
}
//--------------------------------------------------------------------
TMPL template<class CELL> void 
SELF::split(CELL*& left, CELL*& right, int v, double c) {

    //  cout<<"Split begin"<<endl;
    if(v==0) {
        left = new CELL(m_polynomial, BoundingBox(this->xmin(),c, this->ymin(), this->ymax(), this->zmin(), this->zmax()));
        right= new CELL(m_polynomial, BoundingBox(c,this->xmax(), this->ymin(), this->ymax(), this->zmin(), this->zmax()));
    } else if (v==1) {
        left = new CELL(m_polynomial, BoundingBox(this->xmin(),this->xmax(), this->ymin(), c, this->zmin(), this->zmax()));
        right= new CELL(m_polynomial, BoundingBox(this->xmin(),this->xmax(), c, this->ymax(), this->zmin(), this->zmax()));
    } else {//v==2
        left = new CELL(m_polynomial, BoundingBox(this->xmin(),this->xmax(),this->ymin(),this->ymax(), this->zmin(), c));
        right= new CELL(m_polynomial, BoundingBox(this->xmin(),this->xmax(),this->ymin(),this->ymax(), c, this->zmax()));
    }
    //cell3d_split(left,right,this,v,c);
    //cout<<"Subdivide end "<<*left<<" "<<*right<<endl;
    tensor::split(left->m_polynomial, right->m_polynomial, v);
}

//--------------------------------------------------------------------
TMPL void 
SELF::subdivide(Cell*& Left, Cell*& Right, int v, double c) {
    SELF * left=0, * right=0;
    this->split(left,right,v,c);
    Left=left; Right=right;
}

//--------------------------------------------------------------------
TMPL  C
SELF::value(const C& x, const C& y, const  C& z) const {
    return m_polynomial((x-this->xmin())/(this->xmax()-this->xmin()),
                        (y-this->ymin())/(this->ymax()-this->ymin()),
                        (z-this->zmin())/(this->zmax()-this->zmin()));
}
//--------------------------------------------------------------------
TMPL double
SELF::vertex_eval(unsigned sx, unsigned sy, unsigned sz) const {

    int s=0;
    s+= sx*(m_polynomial.rep().env.sz(0)-1)*m_polynomial.rep().env.st(0);
    s+= sy*(m_polynomial.rep().env.sz(1)-1)*m_polynomial.rep().env.st(1);
    s+= sz*(m_polynomial.rep().env.sz(2)-1)*m_polynomial.rep().env.st(2);
    return m_polynomial[s];
}

//--------------------------------------------------------------------
TMPL int
SELF::mc_index() const {

    double isolevel = 0;
    int CubeIndex=0;

    if (vertex_eval(0,1,0) <= isolevel) CubeIndex |= 1;
    if (vertex_eval(1,1,0) <= isolevel) CubeIndex |= 2;
    if (vertex_eval(1,0,0) <= isolevel) CubeIndex |= 4;
    if (vertex_eval(0,0,0) <= isolevel) CubeIndex |= 8;
    if (vertex_eval(0,1,1) <= isolevel) CubeIndex |= 16;
    if (vertex_eval(1,1,1) <= isolevel) CubeIndex |= 32;
    if (vertex_eval(1,0,1) <= isolevel) CubeIndex |= 64;
    if (vertex_eval(0,0,1) <= isolevel) CubeIndex |= 128;

    return CubeIndex;
}

//--------------------------------------------------------------------
TMPL void
SELF::point (Point& R, const Point& A, const Point& B) const {

    Point  U0((A.x()-this->xmin())/(this->xmax()-this->xmin()),
              (A.y()-this->ymin())/(this->ymax()-this->ymin()),
              (A.z()-this->zmin())/(this->zmax()-this->zmin())),
            U1((B.x()-this->xmin())/(this->xmax()-this->xmin()),
               (B.y()-this->ymin())/(this->ymax()-this->ymin()),
               (B.z()-this->zmin())/(this->zmax()-this->zmin()));

    C f0 = this->m_polynomial(U0.x(), U0.y(), U0.z()),
            f1 = this->m_polynomial(U1.x(), U1.y(), U1.z());

    Point M;
    if (abs(f0) < 0.00001) {
        M=U0;
    } else if (abs(f1) < 0.00001) {
        M=U1;
    } else if (abs(f1-f0) < 0.00001) {
        M=U0;
    } else {

        C mu = - f0 / (f1 - f0);

        M.x() = U0.x() + mu * (U1.x() - U0.x());
        M.y() = U0.y() + mu * (U1.y() - U0.y());
        M.z() = U0.z() + mu * (U1.z() - U0.z());

        C  fm = this->m_polynomial(M.x(),M.y(),M.z());

        for(int i=0;i<10;i++) {
            if(f0*fm<0) {
                U1 = M; f1 = fm;
            } else {
                U0 = M; f0 = fm;
            }
            mu = - f0 / (f1 - f0);
            M.x() = U0.x() + mu * (U1.x() - U0.x());
            M.y() = U0.y() + mu * (U1.y() - U0.y());
            M.z() = U0.z() + mu * (U1.z() - U0.z());

            fm = this->m_polynomial(M.x(),M.y(),M.z());
        }
    }

    R.x()=this->xmin()+ M.x()*(this->xmax()-this->xmin());
    R.y()=this->ymin()+ M.y()*(this->ymax()-this->ymin());
    R.z()=this->zmin()+ M.z()*(this->zmax()-this->zmin());

}
//--------------------------------------------------------------------
TMPL bool
SELF::edge_point(Point** CELL, int cube_index) const {

    BoundingBox* bx = (BoundingBox*)this;
    Seq<Point*>  F;
    Polynomial bk, ft;

    tensor::face(bk, equation(), 2, 0);
    if (cube_index & 1) {
        Solver::edge_point(F, bk, Solver::north_back_edge, *bx);
        if(F.size()==0) {
            std::cout<<"Problem in MC0"<<std::endl;
            return false;
        } else
            CELL[0] = F[0];
        F.resize(0);
    }
    if (cube_index & 2) {
        Solver::edge_point(F, bk, Solver::east_back_edge , *bx);
        if(F.size()==0) {
            std::cout<<"Problem in MC1"<<std::endl;
            return false;
        } else
            CELL[1] = F[0];
        F.resize(0);
    }
    if (cube_index & 4) {
        Solver::edge_point(F, bk, Solver::south_back_edge, *bx);
        if(F.size()==0) {
            std::cout<<"Problem in MC2"<<std::endl;
            return false;
        } else
            CELL[2] = F[0];
        F.resize(0);
    }
    if (cube_index & 8) {
        Solver::edge_point(F, bk, Solver::west_back_edge , *bx);
        if(F.size()==0) {
            std::cout<<"Problem in MC3"<<std::endl;
        } else
            CELL[3] = F[0];
        F.resize(0);
    }

    tensor::face(ft,equation(), 2, 1);

    if (cube_index & 16) {
        Solver::edge_point(F, ft, Solver::north_front_edge, *bx);
        if(F.size()==0) {
            std::cout<<"Problem in MC4 "<<std::endl;
            return false;
        } else
            CELL[4] = F[0];
        F.resize(0);
    }

    if (cube_index & 32) {
        Solver::edge_point(F, ft, Solver::east_front_edge , *bx);
        if(F.size()==0) {
            std::cout<<"Problem in MC5"<<std::endl;
            return false;
        } else
            CELL[5] = F[0];
        F.resize(0);
    }

    if (cube_index & 64) {
        Solver::edge_point(F, ft, Solver::south_front_edge, *bx);
        if(F.size()==0) {
            std::cout<<"Problem in MC6"<<std::endl;
            return false;
        } else
            CELL[6] = F[0];
        F.resize(0);
    }

    if (cube_index & 128) {
        Solver::edge_point(F, ft, Solver::west_front_edge , *bx);
        if(F.size()==0)  {
            std::cout<<"Problem in MC7"<<std::endl;
            return false;
        } else
            CELL[7] = F[0];
        F.resize(0);
    }

    tensor::face(ft, equation(), 1, 1);
    if (cube_index & 256) {
        Solver::edge_point(F, ft, Solver::north_west_edge, *bx);
        if(F.size()==0)  {
            std::cout<<"Problem in MC8"<<std::endl;
            return false;
        } else
            CELL[8] = F[0];
        F.resize(0);
    }

    if (cube_index & 512) {
        Solver::edge_point(F, ft, Solver::north_east_edge, *bx);
        if(F.size()==0)  {
            std::cout<<"Problem in MC9"<<std::endl;
            return false;
        } else
            CELL[9] = F[0];
        F.resize(0);
    }

    tensor::face(bk, equation(), 1, 0);

    if (cube_index & 1024) {
        Solver::edge_point(F, bk, Solver::south_east_edge, *bx);
        if(F.size()==0)  {
            std::cout<<"Problem in MC10"<<std::endl;
            return false;
        } else
            CELL[10] = F[0];
        F.resize(0);
    }

    if (cube_index & 2048) {
        Solver::edge_point(F, bk, Solver::south_west_edge, *bx);
        if(F.size()==0)  {
            std::cout<<"Problem in MC11"<<std::endl;
            return false;
        } else
            CELL[11] = F[0];
        F.resize(0);
    }

    return true;
}

//====================================================================
} ; // namespace shape
              } ; // namespace mmx
# undef TMPL
# undef TMPL1
# undef Point
# undef Edge
# undef AlgebraicSurface
# undef SELF 
# endif // SHAPE_IMPLICITCELL_H

---