synaps/topology/TopSbdSurf.h

/***************************************************************************
 *   Copyright (C) 2005 by Chen Liang                                      *
 *   cliang@cs.hku.hk                                                      *
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 *   This program is distributed in the hope that it will be useful,       *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
 *   GNU General Public License for more details.                          *
 *                                                                         *
 *   You should have received a copy of the GNU General Public License     *
 *   along with this program; if not, write to the                         *
 *   Free Software Foundation, Inc.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/

#ifndef synaps_topology_TopSbdSurf_h
#define synaps_topology_TopSbdSurf_h
#include <synaps/linalg/MatrDse.h>
#include <synaps/topology/Cell3d.h>
#include <synaps/topology/TopSbdCurve3d.h>
__BEGIN_NAMESPACE_SYNAPS 
//====================================================================
// cell3d_curve: the structure of the cell that is tailored for building topological graph of
//               the curve defined by two intersection implicit surfaces
template < class C,
           class CELL= Cell3d<C> ,
           class MSLV = SbdSolver<typename CELL::coeff_t,SBDSLV_RDL> >

struct TopSbdSurf3d : public TopSbdCurve3d<C,CELL,MSLV> 
{

  // type definitions //
  
  typedef TopSbdSurf3d<C, CELL,MSLV>          self_t;
  typedef TopSbdCurve3d<C,CELL,MSLV>   base_t;
  
  typedef typename CELL::coeff_t         coeff_t;
  typedef typename CELL::MPOL_T          MPOL_T;
  typedef topology::box3d<coeff_t>       box_t;

  typedef CELL                           cell_t;
  typedef OCTREE::Node<cell_t*>           node_t;
  typedef OCTREE::Octree<cell_t*>         octree_t;
  
  typedef topology::point<C>             point_t;
  typedef topology::point_graph<C>       pointgraph_t;
 
  typedef intersection3d<coeff_t>        intersection_t;
  typedef Seq<intersection3d<coeff_t>*>  intersectionlist_t;
  typedef shared_object<intersection_t*> intersection_ptr;
  typedef Seq<intersection_ptr>          intersection_ptr_list;
  typedef topology::point_graph<coeff_t, intersection_ptr_list> pointgraph_int_t;

  //  typedef topology::point<coeff_t>        point_t;
  //  typedef topology::point_graph<coeff_t>  pointgraph_t;

  typedef MPolDse< coeff_t, mpoldse::bernstein<coeff_t> > mpolb_t;

  static OCTREE::EDGE_TYPE EdgesOfFace[6][4];

  unsigned m_min_dp, m_max_dp;
  MatrDse< Seq<OCTREE::NODE_TYPE> > insert_matr;

  TopSbdSurf3d():base_t (), m_min_dp(0), m_max_dp(6) {} //init_insert_matr();}
  TopSbdSurf3d(unsigned mn, unsigned mx=6):base_t (), m_min_dp(mn), m_max_dp(mx) {} //init_insert_matr();}
  TopSbdSurf3d(coeff_t e):base_t(e), m_min_dp(0), m_max_dp(6) {} //init_insert_matr(); }
  TopSbdSurf3d(coeff_t e, int d): base_t(e,d), m_min_dp(0), m_max_dp(6) {} //init_insert_matr(); }
  //  ~TopSbdSurf3d();



  // member functions //

  bool     compute_point_inside(point_t& p, cell_t* c);

  void     computeintersections(intersectionlist_t& inters, cell_t& c);

  // compute intersection on the dedicated face for the given polynomial
  void     computeinter_face(intersectionlist_t& inters, cell_t& c, int nAxis, int nFace);

  // compute intersection on the edge corresponding to the two faces
  void     computeinter_face(intersectionlist_t& inters, cell_t& c, 
                             int n1Axis, int n1Face, int n2Axis, int n2Face);

  void     computeinter_edge(intersectionlist_t& inters, cell_t& c, 
                             OCTREE::EDGE_TYPE edge);

  void     computeinter_edges(intersectionlist_t& inters, cell_t& c);
  // compute the intersection of  the first polynomial with the edges
  //  void     computeinter_edges(Seq<point_t>& inters, cell_t& c);

  bool     regularity(CELL& c);

  template<class POINTGRAPH>
  void build_mesh(POINTGRAPH& g, node_t* n);
  
  
  void compute_points(node_t* n);

  template<class POINTGRAPH>
  void run(POINTGRAPH& g, CELL* rc);

protected:


  // this function subdivide a given node
  void subdivide(node_t*  n, std::vector<node_t*>& aNodeList);
  // this function subdivide a cell
  void subdivide(CELL* cell, std::vector<CELL*>& subCells);

};
//====================================================================
template< class C, class CELL, class MSLV>
OCTREE::EDGE_TYPE TopSbdSurf3d<C, CELL, MSLV>::EdgesOfFace[6][4]= 
  {{OCTREE::B_E,OCTREE::B_W,OCTREE::B_S,OCTREE::B_N},
   {OCTREE::F_E,OCTREE::F_W,OCTREE::F_S,OCTREE::F_N},
   {OCTREE::B_W,OCTREE::F_W,OCTREE::S_W,OCTREE::N_W},
   {OCTREE::B_E,OCTREE::F_E,OCTREE::S_E,OCTREE::N_E},
   {OCTREE::B_S,OCTREE::F_S,OCTREE::S_W,OCTREE::S_E},
   {OCTREE::B_N,OCTREE::F_N,OCTREE::N_W,OCTREE::N_E}};
//====================================================================
// template< class C, class CELL, class MSLV>
// void TopSbdSurf3d<C, CELL, MSLV>::init_insert_matr()
// {
//   using  namespace OCTREE;
//   insert_matr = MatrDse< Seq<OCTREE::NODE_TYPE> >(8,12);
//   //BS_W
//   insert_matr(B_SW,B_E).push_back(B_SE);
//   insert_matr(B_SW,B_N).push_back(B_SW);
//   insert_matr(B_SW,N_W).push_back(F_NW);
//   insert_matr(B_SW,N_E).push_back(B_SE).push_back(B_NW).push_back(B_NE);
//   insert_matr(B_SW,S_E).push_back(B_SE);
//   insert_matr(B_SW,F_N).push_back(B_NW).push_back(F_SW).push_back(F_NW);
//   insert_matr(B_SW,F_S).push_back(F_SW);
//   insert_matr(B_SW,F_E).push_back(F_SW).push_back(B_SE).push_back(F_SE);
//   insert_matr(B_SW,F_W).push_back(F_SW);

//   //B_SE
//   insert_matr(B_SE,B_N).push_back(B_NE);
//   insert_matr(B_SE,N_E).push_back(B_NE);
//   insert_matr(B_SE,F_N).push_back(B_NE).push_back(F_SE).push_back(F_NE);
//   insert_matr(B_SE,F_S).push_back(F_SE);
//   insert_matr(B_SE,F_E).push_back(F_SE);

//   //F_SW
//   insert_matr(F_SW,N_W).push_back(F_NW);
//   insert_matr(F_SW,N_E).push_back(F_SE).push_back(F_NW).push_back(F_NE);
//   insert_matr(F_SW,S_E).push_back(F_SE);
//   insert_matr(F_SW,F_N).push_back(F_NW);
//   insert_matr(F_SW,F_E).push_back(F_SE);
  
//   //F_SE
//   insert_matr(F_SE,N_E).push_back(F_NE);
//   insert_matr(F_SE,F_N).push_back(F_NE);

//   //B_NW
//   insert_matr(B_NW,B_E).push_back(B_NE);
//   insert_matr(B_NW,N_E).push_back(B_NE);
//   insert_matr(B_NW,F_N).push_back(F_NW);
//   insert_matr(B_NW,F_W).push_back(F_NW);
//   insert_matr(B_NW,F_E).push_back(F_NW).push_back(B_NE).push_back(F_NE);

//   //B_NE
//   insert_matr(B_NE,F_N).push_back(F_NE);
//   insert_matr(B_NE,F_E).push_back(F_NE);

//   //F_NW
//   insert_matr(F_NW,N_E).push_back(F_NE);
//   insert_matr(F_NW,F_E).push_back(F_NE);

// }

// member function definitions  //

//====================================================================

/***************************************/
/*               Z^                    */
/*                |_____               */
/*               /|    /|              */
/*     N  ^     /_|_B_/ |              */
/*        |    |  |___|_|___> Y        */
/*        |    | /F   | /              */
/*             |/_____|/               */
/*            |/_                      */
/*           X                         */
/*                                     */
/*    Fig.1  The coordinate system     */
/***************************************/

template< class C, class CELL, class MSLV>
void TopSbdSurf3d<C, CELL, MSLV>::subdivide(node_t* n, std::vector<node_t*>& aNodeList)
{
  using namespace OCTREE;

  std::vector<node_t*> aNodeNeigbrs;

    
  cell_t* cell = n->GetCell();
  std::vector<cell_t*> subCells;

  // subdivide a cell into 8 subcells
  subdivide(cell,subCells);
  

    // compute the intersections on the cell borders (Dynamic Programming applied here
    // to avoid repetitive computation, see arrangements of _cellidx/_faceidx/_shareidx)
//     for (int i = 0; i < 8*MAX_FACEIDX; i++)
//     {
//         intersectionlist_t aInters;
//         int nCellIdx=_cellidx[i], nFaceIdx = _faceidx[i], nShareIdx = _shareidx[i];

//         if(nFaceIdx >= 0)
//         {
//             //       computeinter_face(aInters, *subCells[nCellIdx], nFaceIdx>>1, nFaceIdx%2);
//             // the cell shares a face with this cell (i.e. subcells[_cellidx[i]]) should also
//             // share the some intersections (on a particular face) with this cell
//             // cell->mergelist(subCells[nCellIdx]->m_aInterList, aInters);
//             // cell->mergelist(subCells[nShareIdx]->m_aInterList, aInters);
//             // mark the this face (i.e. face nFaceIdx of cell nCellIdx) as 'computed'
//          //            subCells[nCellIdx]->m_FaceComputed[nFaceIdx] = true;
//             // for two subcells that shares a face, the shared face should of coz be parallel,
//             // so their face indices should have same rounded value when devided by 3 (since
//             // the face indices are (YZ_Plane_0, YZ_Plane_1, XZ_Plane0, XZ_Plane1, XY_Plane0,
//             // XY_Plane1)). And of coz, they have different mod against 2.
//             // by a simple deduction, f(x)=x+(x%2==1?-1:1) gives the correct face index for face x;
//             // now we mark the shared face also as 'computed'
//          // assert(subCells[nShareIdx]->m_FaceComputed[nFaceIdx+(nFaceIdx%2?-1:1)] == false);
//             // if (subCells[nShareIdx]->m_FaceComputed[nFaceIdx+(nFaceIdx%2?-1:1)])
//             //    std::cout << "assertion failed at depth " << cell->m_nDepth << " nCell " << nCellIdx << " nFace " << nFaceIdx << " nShare " << nShareIdx << " i"  << i << ", cell " << nShareIdx << " face " << (nFaceIdx+(nFaceIdx%2?-1:1)) << " has already been computed." << std::endl;
//          // subCells[nShareIdx]->m_FaceComputed[nFaceIdx+(nFaceIdx%2?-1:1)] = true;
//         }
//     }


  // encapsulate the cells with the node of the octree
  aNodeList.push_back(new node_t(B_SW, n, subCells[B_SW]));
  aNodeList.push_back(new node_t(F_SW, n, subCells[F_SW]));
  aNodeList.push_back(new node_t(B_SE, n, subCells[B_SE]));
  aNodeList.push_back(new node_t(F_SE, n, subCells[F_SE]));
  aNodeList.push_back(new node_t(B_NW, n, subCells[B_NW]));
  aNodeList.push_back(new node_t(F_NW, n, subCells[F_NW]));
  aNodeList.push_back(new node_t(B_NE, n, subCells[B_NE]));
  aNodeList.push_back(new node_t(F_NE, n, subCells[F_NE]));

//   for(unsigned i=0;i<8;i++)
//     {
//       intersectionlist_t aInters;
//       computeinter_edges(aInters, *subCells[i]);
//       cell->mergelist(subCells[i]->m_aInterList, aInters);
//     }


//   node_t* Nbg[6];
//   Nbg[B]=n->getBackNeighborNode();
//   Nbg[F]=n->getFrontNeighborNode();
//   Nbg[N]=n->getNorthNeighborNode();
//   Nbg[S]=n->getSouthNeighborNode();
//   Nbg[E]=n->getEastNeighborNode();
//   Nbg[W]=n->getWestNeighborNode();


//   MatrDse<node_t*> ngbrs(8,12);
//   ngbrs(B_SW,B_N)=Nbg[B];
//   ngbrs(B_SW,B_E)=Nbg[B];
//   ngbrs(B_SW,N_W)=Nbg[W];
//   ngbrs(B_SW,F_W)=Nbg[W];
//   ngbrs(B_SW,F_S)=Nbg[S];
//   ngbrs(B_SW,S_E)=Nbg[S];

//   ngbrs(B_SE,B_N)=Nbg[B];
//   ngbrs(B_SE,N_E)=Nbg[E];
//   ngbrs(B_SE,F_E)=Nbg[E];
//   ngbrs(B_SE,F_S)=Nbg[S];

//   ngbrs(F_SW,N_W)=Nbg[W];
//   ngbrs(F_SW,F_N)=Nbg[F];
//   ngbrs(F_SW,F_E)=Nbg[F];
//   ngbrs(F_SW,S_E)=Nbg[S];

//   ngbrs(F_SE,N_E)=Nbg[E];
//   ngbrs(F_SE,F_N)=Nbg[F];

//   ngbrs(B_NW,F_W)=Nbg[W];
//   ngbrs(B_NW,B_E)=Nbg[B];
//   ngbrs(B_NW,N_E)=Nbg[N];
//   ngbrs(B_NW,F_N)=Nbg[N];

//   ngbrs(F_NE,B_E)=Nbg[E];
//   ngbrs(F_NE,F_W)=Nbg[F];
//   ngbrs(F_NE,N_W)=Nbg[N];
//   ngbrs(F_NE,B_N)=Nbg[N];

//   cell_t* clnbr;

//   for(unsigned i=0;i<8;i++)
//     for(unsigned j=0;j<12;j++)
//       {
//      intersectionlist_t L;
//      computeinter_edge(L,*subCells[i],(EDGE_TYPE)j);
//      //      PRINT_DEBUG(" L sz "<<L.size());
//      if(ngbrs(i,j) != NULL && L.size())
//        {
//          clnbr = ngbrs(i,j)->GetCell();
//          clnbr->mergelist(clnbr->m_aInterList,L);
//        }
//       }
}
//--------------------------------------------------------------------
/***************************************/
/*          ^ Z                        */
/*          |                          */
/*          |                          */
/*          |                          */
/*          |4----6                    */
/*         /||___/|_______> Y          */
/*        5/-0--7 2                    */
/*        /|/   |/                     */
/*       / 1----3                      */
/*     \/_                             */
/*     X                               */
/*                                     */
/*    Fig.2  The layout of the         */
/*        subdivided cells             */
/***************************************/

/***************************************/
/*          ^ Z                        */
/*          |   5     0                */
/*          |  |     /                 */
/*          |  v   |/_                 */
/*          |------                    */
/*         /||___/|_______> Y          */
/*  2-->  -/ ---- | <--3               */
/*        /|/ 1 |/                     */
/*       / ------                      */
/*     \/_    ^                        */
/*     X      |                        */
/*             4                       */
/*    Fig.3  Faces index of a cell     */
/***************************************/
// this function spawns 8 evenly subdivided child cell of this call
template < class C, class CELL, class MSLV >
void TopSbdSurf3d<C, CELL, MSLV>::subdivide(CELL* cell,
                                            std::vector<CELL*>& subCells)
{
  typedef typename CELL::intersection_t intersection_t;
  typedef typename CELL::intersectionlist_t intersectionlist_t;
  typedef typename CELL::MPOL_T MPOL_T;
    // define the subdivision points
#ifdef SYNAPS_SBD3D_RNDFACTORWHENSUBDIVISION
    // a perturbation ranged from 0.005-0.1
    C t[] = {(C(rand())/C(RAND_MAX)-C(0.5))*C(0.095)+C(0.5),
             (C(rand())/C(RAND_MAX)-C(0.5))*C(0.095)+C(0.5),
             (C(rand())/C(RAND_MAX)-C(0.5))*C(0.095)+C(0.5)};
    t[0]+= (t[0]>0.5?0.005:-0.005);
    t[1]+= (t[1]>0.5?0.005:-0.005);
    t[2]+= (t[2]>0.5?0.005:-0.005);
    C mid[] = {cell->m_cBox.x0*(C(1)-t[0])+cell->m_cBox.x1*t[0], cell->m_cBox.y0*(C(1)-t[1])+cell->m_cBox.y1*t[1], cell->m_cBox.z0*(C(1)-t[2])+cell->m_cBox.z1*t[2]};
#else
    C mid[] = {(cell->m_cBox.x0+cell->m_cBox.x1)/C(2), (cell->m_cBox.y0+cell->m_cBox.y1)/C(2), (cell->m_cBox.z0+cell->m_cBox.z1)/C(2)};
#endif
    C box[] = {cell->m_cBox.x0, cell->m_cBox.x1, cell->m_cBox.y0, cell->m_cBox.y1, cell->m_cBox.z0, cell->m_cBox.z1};

    // create new cells
    subCells.push_back(new CELL(cell->nbpol(), box[0], mid[0], box[2], mid[1], box[4], mid[2]));
    subCells.push_back(new CELL(cell->nbpol(), mid[0], box[1], box[2], mid[1], box[4], mid[2]));
    subCells.push_back(new CELL(cell->nbpol(), box[0], mid[0], mid[1], box[3], box[4], mid[2]));
    subCells.push_back(new CELL(cell->nbpol(), mid[0], box[1], mid[1], box[3], box[4], mid[2]));
    subCells.push_back(new CELL(cell->nbpol(), box[0], mid[0], box[2], mid[1], mid[2], box[5]));
    subCells.push_back(new CELL(cell->nbpol(), mid[0], box[1], box[2], mid[1], mid[2], box[5]));
    subCells.push_back(new CELL(cell->nbpol(), box[0], mid[0], mid[1], box[3], mid[2], box[5]));
    subCells.push_back(new CELL(cell->nbpol(), mid[0], box[1], mid[1], box[3], mid[2], box[5]));

    // update the depth count for each subdivided cells
    for (int i = 0; i < 8; i++)
    {
        subCells[i]->m_nDepth = cell->m_nDepth+1;
    }

  //   // the subcells inherits the previous computed intersections on its parent cell
//     for (int i = 0; i < cell->m_aInterList.size(); i++)
//     {
//         intersection_t* inter = cell->m_aInterList[i];

//         // assign the intersections on the parent cell to subdivided cells, according to the plane
//         // they are on and the planar coordinates
//      int nidx = inter->nPlane; // testing the plane YZ (nidx=0) or XZ (nidx=1) or XZ (nidx=2)
//         int *cidx=subcellidx+(nidx<<3); // the right set of subcell indices to be assigned each intersection
//      int axisidx0 = plane_axisidx[nidx<<1], axisidx1 = plane_axisidx[(nidx<<1)+1]; // the first and the second axis indices, {1, 2} for YZ plane, {0, 2} for XZ plane and {0, 1} for XY plane

//         // assert(_is_verysmall(inter->pnt[nidx]-box[nidx<<1]) || _is_verysmall(box[(nidx<<1)+1]-inter->pnt[nidx]));
//      // assert(inter->pnt[axisidx0] >= box[axisidx0<<1] && inter->pnt[axisidx0] <= box[(axisidx0<<1)+1] &&
//         //        inter->pnt[axisidx1] >= box[axisidx1<<1] && inter->pnt[axisidx1] <= box[(axisidx1<<1)+1]);

//         if (inter->pnt[axisidx0] <= mid[axisidx0]) // if smaller than the first mid-point
//         {
//             if (inter->pnt[axisidx1] <= mid[axisidx1]) // if smaller than the second mid-point
//             {
//                 if (inter->pnt[nidx]-box[nidx<<1] <= box[(nidx<<1)+1]-inter->pnt[nidx]) // if the plane corresponds to x (or y, or z) = a, where x (or y, or z) is in [a, b]
//                     subCells[cidx[0]]->m_aInterList.push_back(inter);
//                 else // if the plane corresponds to x (or y, or z) = b, where x (or y, or z) is in [a, b]
//                     subCells[cidx[1]]->m_aInterList.push_back(inter);
//             }
//             if (inter->pnt[axisidx1] >= mid[axisidx1]) // if larger than the first mid-point
//             {
//                 if (inter->pnt[nidx]-box[nidx<<1] <= box[(nidx<<1)+1]-inter->pnt[nidx]) // if the plane corresponds to x (or y, or z) = a, where x (or y, or z) is in [a, b]
//                     subCells[cidx[2]]->m_aInterList.push_back(inter);
//                 else
//                     subCells[cidx[3]]->m_aInterList.push_back(inter);
//             }
//         }
//         if (inter->pnt[axisidx0] >= mid[axisidx0]) 
//        // if lager than the first mid-point
//         {
//             if (inter->pnt[axisidx1] <= mid[axisidx1]) 
//            // if smaller than the second  mid-point
//             {
//                 if (inter->pnt[nidx]-box[nidx<<1] <= box[(nidx<<1)+1]-inter->pnt[nidx]) // if the plane corresponds to x (or y, or z) = a, where x (or y, or z) is in [a, b]
//                     subCells[cidx[4]]->m_aInterList.push_back(inter);
//                 else // if the plane corresponds to x (or y, or z) = b, where x (or y, or z) is in [a, b]
//                     subCells[cidx[5]]->m_aInterList.push_back(inter);
//             }
//             if (inter->pnt[axisidx1] >= mid[axisidx1])  
//            // if larger than the second mid-point
//             {
//                 if (inter->pnt[nidx]-box[nidx<<1] <= box[(nidx<<1)+1]-inter->pnt[nidx])
//                     subCells[cidx[6]]->m_aInterList.push_back(inter);
//                 else // if the plane corresponds to x (or y, or z) = b, where x (or y, or z) is in [a, b]
//                     subCells[cidx[7]]->m_aInterList.push_back(inter);
//             }
//         }
//     }

  //   // update the 'computed faces' flag for all subcells
//     for (int i = 0; i < 8*MAX_FACEIDX; i++)
//     {
//         int nCellIdx=_cellidx[i], nInheritIdx = _inheritidx[i];
//         // assert that a face never inherit from more than one parent cell's face
//         assert((subCells[nCellIdx]->m_FaceComputed[nInheritIdx]) == false);
//         /*if ((subCells[nCellIdx]->m_FaceComputed[nInheritIdx]))
//             std::cout << "assertion failed at depth " << cell->m_nDepth << std::endl;*/
//      //        if (cell->m_FaceComputed[nInheritIdx]) // if the parent cell's face this cell inherited from has 'computed' flag
//      //            subCells[nCellIdx]->m_FaceComputed[nInheritIdx] = true;
//     }

    // split the representations for each surface
    for (int i = 0; i < cell->nbpol(); i++)
    {
        MPOL_T splited[8];
        int vr_idx[3] = {VARIDX(cell->m_F[i], 0), VARIDX(cell->m_F[i], 1), VARIDX(cell->m_F[i], 2)};

        subCells[7]->m_F[i] = cell->m_F[i];
#ifdef SYNAPS_SBD3D_RNDFACTORWHENSUBDIVISION
        // subdivide along z direction with a perturbation
        mpoldse::split(subCells[3]->m_F[i].rep(), subCells[7]->m_F[i].rep(), vr_idx[2], t[2]);

        // subdivide along y direction with a perturbation
        mpoldse::split(subCells[1]->m_F[i].rep(), subCells[3]->m_F[i].rep(), vr_idx[1], t[1]);
        mpoldse::split(subCells[5]->m_F[i].rep(), subCells[7]->m_F[i].rep(), vr_idx[1], t[1]);

        // subdivide along x direction with a perturbation
        mpoldse::split(subCells[0]->m_F[i].rep(), subCells[1]->m_F[i].rep(), vr_idx[0], t[0]);
        mpoldse::split(subCells[2]->m_F[i].rep(), subCells[3]->m_F[i].rep(), vr_idx[0], t[0]);
        mpoldse::split(subCells[4]->m_F[i].rep(), subCells[5]->m_F[i].rep(), vr_idx[0], t[0]);
        mpoldse::split(subCells[6]->m_F[i].rep(), subCells[7]->m_F[i].rep(), vr_idx[0], t[0]);
#else
        // subdivide along z direction
        mpoldse::split(subCells[3]->m_F[i].rep(), subCells[7]->m_F[i].rep(), vr_idx[2]);

        // subdivide along y direction
        mpoldse::split(subCells[1]->m_F[i].rep(), subCells[3]->m_F[i].rep(), vr_idx[1]);
        mpoldse::split(subCells[5]->m_F[i].rep(), subCells[7]->m_F[i].rep(), vr_idx[1]);

        // subdivide along x direction
        mpoldse::split(subCells[0]->m_F[i].rep(), subCells[1]->m_F[i].rep(), vr_idx[0]);
        mpoldse::split(subCells[2]->m_F[i].rep(), subCells[3]->m_F[i].rep(), vr_idx[0]);
        mpoldse::split(subCells[4]->m_F[i].rep(), subCells[5]->m_F[i].rep(), vr_idx[0]);
        mpoldse::split(subCells[6]->m_F[i].rep(), subCells[7]->m_F[i].rep(), vr_idx[0]);
#endif
    }
        
}
//====================================================================
template< class C, class CELL, class MSLV>
template< class POINTGRAPH>
void TopSbdSurf3d<C, CELL, MSLV>::run(POINTGRAPH& g, CELL* rc)
{
  typedef CELL cell_t;
    node_t* root = this->m_cOcTree.GetRoot();

    //    computeintersections(rc->m_aInterList, *rc);
    //    std::cout << "no. of roots at root: " << rc->m_aInterList.size() << std::endl;

    root->SetCell(rc);

    // step3:
    // subdivision
    //    PRINT_DEBUG("Subdivision");
    std::list<node_t*> aNodeQueue;
    std::list<node_t*> aNodeKept;
    std::vector<node_t*> aSingular;
    pointgraph_int_t g_tmp;
    int nCellProc = 0;

    aNodeQueue.push_front(root);

    while(aNodeQueue.size())
    {
        // process the first node in the queue
        node_t* n = *aNodeQueue.begin();
        aNodeQueue.pop_front();

        bool bSubdivide = false;
        bool bBuildGraph = false;

        // if the cell is regular, i.e. the topology of the surface within the cell is certified
        if (regularity(*n->GetCell()))
        {
            if (is_variety_in(*n->GetCell())) // it may contains the variety
            {
                if (n->GetCell()->m_nDepth < this->m_nMinSbd) // if we did less than the minimum number of subdivision
                {
                    // subdivide further
                    bSubdivide = true;
                }
                else // build the graph
                {
                    bBuildGraph = true;
                }
            }
            else // if regular yet does not contain the surface, of no interest for further subdivision
            {
                // its safe to say that not any facet of this cell should contain a intersection point
                memset(n->GetCell()->m_FaceComputed, true, sizeof(bool)*6);
            }
        }
        else // if the cell is not regular
        {
            if (n->GetCell()->m_nDepth> m_max_dp) // if the cell size is smaller than epsilon
            {
                // keep a count of singular cells
                if (is_variety_in(*n->GetCell()))
                {
                    bBuildGraph = true;
                    aSingular.push_back(n);
                }
                else
                  memset(n->GetCell()->m_FaceComputed, true, sizeof(bool)*6);
            }
            else // if the cell size is still larger than the set threshold
            {
                // further subdivision if it may contains the intersection
                if (is_variety_in(*n->GetCell()))
                    bSubdivide = true;
                else
                  memset(n->GetCell()->m_FaceComputed, true, sizeof(bool)*6);
            }
        }

        // in fact, bSubdivide and bBuildGraph cannot be both true (but can be both false)
        // assert(!(bSubdivide && bBuildGraph))
        if (bSubdivide)
        {
            // append the subdivided node at the back of the queue
            // in doing this, we can guarantee that larger nodes are processed first
            std::vector<node_t*> aNodeList;
            subdivide(n, aNodeList);
            for (int i = 0; i < (int)aNodeList.size(); i++)
                aNodeQueue.push_back(aNodeList[i]);
        }
        else if (bBuildGraph)
        {
          aNodeKept.push_front(n);
        }

        // a counter
        nCellProc++;
    }
    PRINT_DEBUG("Computing points");
    typedef typename std::list<node_t*>::const_iterator node_iterator_t;
    for(node_iterator_t it=aNodeKept.begin(); it != aNodeKept.end(); it++)
      {
        compute_points(*it);
      }
    //Once all cells are computed, we build the mesh.
    PRINT_DEBUG("Building mesh");
    for(node_iterator_t it=aNodeKept.begin(); it != aNodeKept.end(); it++)
      {
        build_mesh(g,(*it));
      }

    std::cout << "Number of Cells     : " << nCellProc << std::endl;
    std::cout << "Number of Sing Cells: " << aSingular.size() << std::endl;
}

//====================================================================
template <class C, class CELL, class MSLV >
bool TopSbdSurf3d<C, CELL, MSLV>::compute_point_inside(point_t& p, cell_t* cell)
{
  using namespace OCTREE;
  std::vector<cell_t*> subCells;
  subdivide(cell,subCells);
  
  intersectionlist_t L;
  computeinter_edge(L,*subCells[B_SW],N_E);
  if(L.size())
    {
      for(unsigned j=0;j<3;j++) p[j]=L[0]->pnt[j]; return true; 
    }

  computeinter_edge(L,*subCells[B_SW],F_N);
  if(L.size())
    {
      for(unsigned j=0;j<3;j++) p[j]=L[0]->pnt[j]; return true; 
    }

  computeinter_edge(L,*subCells[B_SW],F_E);
  if(L.size())
    {
      for(unsigned j=0;j<3;j++) p[j]=L[0]->pnt[j]; return true; 
    }

  computeinter_edge(L,*subCells[F_NE],B_W);
  if(L.size())
    {
      for(unsigned j=0;j<3;j++) p[j]=L[0]->pnt[j]; return true; 
    }

  computeinter_edge(L,*subCells[F_NE],B_S);
  if(L.size())
    {
      for(unsigned j=0;j<3;j++) p[j]=L[0]->pnt[j]; return true; 
    }

  computeinter_edge(L,*subCells[F_NE],S_W);
  if(L.size())
    {
      for(unsigned j=0;j<3;j++) p[j]=L[0]->pnt[j]; return true; 
    }

  return false;
}
//--------------------------------------------------------------------
template <class C, class CELL, class MSLV >
void TopSbdSurf3d<C, CELL, MSLV>::computeintersections(intersectionlist_t& inters, CELL& c)
{

  // for each facet of the cell, compute its intersections with the curve
  for (int v = 0; v < 3; v++) // for each axis
    {
      for (int f = 0; f < 2; f++) // for each facet perpendicular to that axis
        {         // compute the intersections on this facet if it has not yet been computed
          if (!c.m_FaceComputed[v*2+f])
            {
              computeinter_face(inters, c, v, f);
              //c.m_FaceComputed[v*2+f] = true;
            }
        }
    }
}
//--------------------------------------------------------------------
template < class C, class CELL, class MSLV >
void TopSbdSurf3d<C, CELL, MSLV>::computeinter_face(intersectionlist_t& inters, 
                                                    CELL& c, 
                                                    int nAxis, int nFace)
{
    PRINT_DEBUG("Points on face "<< nAxis<<" "<<nFace);
    unsigned N=c.m_nT;
    if (N==1)      return;

    // build a system of equations
    std::vector<MPOL_T> facepol(N+3); // = {m_F[0], m_F[1]};
    std::list<MPOL_T> eqnList;
    unsigned sz=0;
    for(unsigned i=0;i<N;i++)
      {
        // obtain the polynomial representing the face
        face(facepol[i].rep(), c.m_F[i].rep(), nAxis, nFace);
        eqnList.push_back(facepol[i]);
        sz++;
      }
    MPOL_T dp[3];
    for (unsigned i = 0; i < 3; i++)
      {
        diff(dp[i],c.m_F[0],VARIDX(c.m_F[0], i));
        face(facepol[N+i].rep(), dp[i].rep(), nAxis, nFace);
        sz++;
      }
    PRINT_DEBUG("Solving "<<sz<<" equations");
    // solve the equation using a 2-variate solver
    MSLV solver(&facepol[0], &facepol[0]+sz, 1e-4);
    PRINT_DEBUG("Solving "<<sz<<" equations");
    std::vector<C> results;
    int nSol, nVar;
    solver.solve(results, nSol, nVar, eqnList.begin(), eqnList.end());
    PRINT_DEBUG("Nb of sol: "<<nSol<<" Nb of var: "<<nVar);
    //    assert(nVar==2);
    // rescale the solution to the coordinates of the cell
    // also update the multiplicity accordingly
    for (int s = 0; s < nSol; s++)
    {
        C x[3] = {(C)nFace,(C)nFace,(C)nFace}; // the 3D coord of the point

        switch(nAxis)
        {
        case PLANEYZ: // yz-planes
            x[1] = (C)(results[nVar*s]);
            x[2] = (C)(results[nVar*s+1]);
            break;
        case PLANEXZ: // xz-planes
            x[0] = (C)(results[nVar*s]);
            x[2] = (C)(results[nVar*s+1]);
            break;
        case PLANEXY: // xy-planes
            x[0] = (C)(results[nVar*s]);
            x[1] = (C)(results[nVar*s+1]);
            break;
        }
        // determine the multiplicity
        //int nMaxDeg = c.m_T[nAxis].rep().env.sz(VARIDX(c.m_T[nAxis], nAxis)); // the degree of nAxis_th variable of nAxis_th component of t = {t1, t2, t3} = gradF cross gradG
        int nMulti =  2;

        // currently, only account for multiplicity 1 or 2...
        //for (int v = 0; v< nMaxDeg; v++)
        if(false)
          {
            C val[3] = {-1,-1,-1};
            eval(val[0], c.m_F[c.m_nT], x); // evaluate the t_nAxis(x0, y0, z0) at the intersection (x0, y0, z0)
            eval(val[1], c.m_F[c.m_nT+1], x); // evaluate the t_nAxis(x0, y0, z0) at the intersection (x0, y0, z0)
            eval(val[2], c.m_F[c.m_nT+2], x); // evaluate the t_nAxis(x0, y0, z0) at the intersection (x0, y0, z0)

            /*if (abs(val[nAxis]) < 0.1)
            {
                std::cout << "sin val: " << val[nAxis];
                val[nAxis] /= sqrt(val[0]*val[0]+val[1]*val[1]+val[2]*val[2]); // normalize the vector
                std::cout << " nor: " << val[nAxis] << std::endl;
            }
            else
            {
                val[nAxis] /= sqrt(val[0]*val[0]+val[1]*val[1]+val[2]*val[2]); // normalize the vector
            }*/

            val[nAxis] /= sqrt(val[0]*val[0]+val[1]*val[1]+val[2]*val[2]); // normalize the vector

            if (_is_verysmall(std::abs(val[nAxis])))
            {
                nMulti = 2;
                //std::cout << "multiple point: val =" << val << ",  f(" << nAxis << ")=" << c.m_T[nAxis] << std::endl;
            }
          }


        // scale to the cell's coordinate
        point_t l(3u);
        l[0] = c.m_cBox.x0*((C)1-x[0]) + c.m_cBox.x1*x[0];
        l[1] = c.m_cBox.y0*((C)1-x[1]) + c.m_cBox.y1*x[1];
        l[2] = c.m_cBox.z0*((C)1-x[2]) + c.m_cBox.z1*x[2];


        intersection_t *ni = new intersection_t(l, nAxis, nMulti);
        inters.push_back(ni);
        ni->nDepth = c.m_nDepth;
        // PRINT_DEBUG("new point");
    }
    PRINT_DEBUG("End point on face");
}
//--------------------------------------------------------------------
template < class C, class CELL, class MSLV >
void TopSbdSurf3d<C, CELL, MSLV>::computeinter_face(intersectionlist_t& inters, 
                                                CELL& c, 
                                                int n1Axis, int n1Face,
                                                int n2Axis, int n2Face)
{
  assert(n1Axis<n2Axis);
  MPOL_T pf,p[1];
  //PRINT_DEBUG(" ");
  //PRINT_DEBUG("FaceFace "<<n1Axis<<" "<<n1Face<<" "<<n2Axis<<" "<<n2Face);

  face(pf.rep(), c.m_F[0].rep(),n1Axis,n1Face); 
  face(p[0].rep(), pf.rep(),n2Axis-1,n2Face); 

  //  PRINT_DEBUG(p[0]);

  MSLV slv(p, p+1, 1e-4);

  std::vector<coeff_t> sol;
  std::list<MPOL_T> eqnList; eqnList.push_back(p[0]);
  int nSol, nVar;
  slv.solve(sol, nSol, nVar, eqnList.begin(), eqnList.end());
  assert(nVar==1);
  //  if(sol.size()) { VECTOR::print(std::cout,sol)<<std::endl;}
  point_t pt(3u);
  for(unsigned j=0;j<sol.size();j++)
    {
      for(unsigned i=0;i<3;i++)
        {
          if(i==n1Axis)
            switch(i)
              {
              case 0:
                pt[0]=(n1Face?c.m_cBox.x1:c.m_cBox.x0);
              case 1:
                pt[1]=(n1Face?c.m_cBox.y1:c.m_cBox.y0);
              case 2:
                pt[2]=(n1Face?c.m_cBox.z1:c.m_cBox.z0);
              }
          else if(i==n2Axis)
            switch(i)
              {
              case 0:
                pt[0]=(n2Face?c.m_cBox.x1:c.m_cBox.x0);
              case 1:
                pt[1]=(n2Face?c.m_cBox.y1:c.m_cBox.y0);
              case 2:
                pt[2]=(n2Face?c.m_cBox.z1:c.m_cBox.z0);
              }
          else
            switch(i)
              {
              case 0:
                pt[0]=c.m_cBox.x0*((coeff_t)1-sol[j]) + c.m_cBox.x1*sol[j];
              case 1:
                pt[1]=c.m_cBox.y0*((coeff_t)1-sol[j]) + c.m_cBox.y1*sol[j];
              case 2:
                pt[2]=c.m_cBox.z0*((coeff_t)1-sol[j]) + c.m_cBox.z1*sol[j];
              }
        }
      int nMulti = 1;
      C val;
      unsigned sz=p[0].rep().size()-1;
      if(sz>1)
        {
          C dp[sz];
          for(unsigned u=0;u<sz;u++) dp[u]=p[0][u+1]-p[0][u];
          mpoldse::decasteljau(dp,sz,sol[j]);
          if(_is_singular(dp[0]))
            {
              //              PRINT_DEBUG(p[0]);
              PRINT_DEBUG("Multiple point: "<<pt<<" "<<dp[0]<<
                          " ["<<c.m_cBox.x0<<" "<<c.m_cBox.x1<<" "<<c.m_cBox.y0<<" "<<c.m_cBox.y1<<" "<<c.m_cBox.z0<<" "<<c.m_cBox.z1<<"]");

              nMulti=2;
            }
        }
      intersection_t *ni = new intersection_t(pt, n1Axis , nMulti);
      inters.push_back(ni);
      ni->nDepth = c.m_nDepth;
      //      nbpt++;
    }


}
//====================================================================
template < class C, class CELL, class MSLV >
void TopSbdSurf3d<C, CELL, MSLV>::computeinter_edge(intersectionlist_t& inters, 
                                                     CELL& c,
                                                     OCTREE::EDGE_TYPE edge)
{
  using namespace OCTREE;
  switch(edge)
    {
    case B_W:
      computeinter_face(inters,c,0,0,1,0);
      break;
    case B_E:
      computeinter_face(inters,c,0,0,1,1);
      break;
    case B_S:
      computeinter_face(inters,c,0,0,2,0);
      break;
    case B_N:
      computeinter_face(inters,c,0,0,2,1);
      break;
    case F_W:
      computeinter_face(inters,c,0,1,1,0);
      break;
    case F_E:
      computeinter_face(inters,c,0,1,1,1);
      break;
    case F_S:
      computeinter_face(inters,c,0,1,2,0);
      break;
    case F_N:
      computeinter_face(inters,c,0,1,2,1);
      break;
    case S_W:
      computeinter_face(inters,c,1,0,2,0);
      break;
    case N_W:
      computeinter_face(inters,c,1,0,2,1);
      break;
    case S_E:
      computeinter_face(inters,c,1,1,2,0);
      break;
    case N_E:
      computeinter_face(inters,c,1,1,2,1);
      break;
    }
}
//====================================================================
template < class C, class CELL, class MSLV >
void TopSbdSurf3d<C, CELL, MSLV>::computeinter_edges(intersectionlist_t& inters, 
                                                 CELL& c)
{
  //  PRINT_DEBUG("-------------------------------");
  //  PRINT_DEBUG(c.m_cBox.x0<<" "<<c.m_cBox.x1<<" "<<c.m_cBox.y0<<" "<<c.m_cBox.y1<<" "<<c.m_cBox.z0<<" "<<c.m_cBox.z1);
//   for(unsigned i0=0;i0<3;i0++)
//     for(unsigned i1=i0+1;i1<3;i1++)
//       for(unsigned j0=0;j0<2;j0++)
//      for(unsigned j1=0;j1<2;j1++)
//        computeinter_face(inters,c,i0,j0,i1,j1);
  using namespace OCTREE;
  for(unsigned i=0;i<12;i++)
    computeinter_edge(inters,c,(EDGE_TYPE)i);
}
//--------------------------------------------------------------------
template < class POINT, class CELL>
unsigned Computeinter_edges(VectDse<Seq<POINT> >& inters, 
                            const CELL& c)
{

  typedef typename CELL::coeff_t         coeff_t;
  typedef typename CELL::MPOL_T          MPOL_T;
  typedef POINT                          point_t;

  unsigned nbpt=0;
  MPOL_T lp[3],px0,px1,pz0,pz1;
  int nSol, nVar;
  point_t pt(3u);

  PRINT_DEBUG("================================");
  //  PRINT_DEBUG(c.m_cBox.x0<<" "<<c.m_cBox.x1<<" "<<c.m_cBox.y0<<" "<<c.m_cBox.y1<<" "<<c.m_cBox.z0<<" "<<c.m_cBox.z1);
  face(px0.rep(), c.m_F[0].rep(),0,0); 
  face(px1.rep(), c.m_F[0].rep(),0,1); 


  face(lp[2].rep(), px0.rep(),0,0); 
  //  PRINT_DEBUG(" ");
  //  PRINT_DEBUG("X0Y0 "<<lp[2]);

  SbdSolver<coeff_t,SBDSLV_RDL> slv2(lp+2, lp+3, 1e-4);
  std::vector<coeff_t> sol;
  std::list<MPOL_T> eqnList; eqnList.push_back(lp[2]);
  slv2.solve(sol, nSol, nVar, eqnList.begin(), eqnList.end());

  assert(nVar==1);
  //  if(sol.size()) { VECTOR::print(std::cout,sol)<<std::endl;}

  pt[0]=c.m_cBox.x0;
  pt[1]=c.m_cBox.y0;
  for(unsigned j=0;j<sol.size();j++)
    {
       pt[2]= c.m_cBox.z0*((coeff_t)1-sol[j]) + c.m_cBox.z1*sol[j];
       inters[0].push_back(pt);
       //  PRINT_DEBUG(pt<<" |"<<sol[j]<<"|  "<<c.m_cBox.x0<<" "<<c.m_cBox.x1<<" "<<c.m_cBox.y0<<" "<<c.m_cBox.y1<<" "<<c.m_cBox.z0<<" "<<c.m_cBox.z1);
       //  PRINT_DEBUG(pt);
       nbpt++;
    }

  face(lp[2].rep(), px0.rep(),0,1); 
  //  PRINT_DEBUG(" ");
  //  PRINT_DEBUG("X0Y1 "<<lp[2]);
  (*eqnList.begin())=lp[2];
  slv2.solve(sol, nSol, nVar, eqnList.begin(), eqnList.end());
  //  if(sol.size()) { VECTOR::print(std::cout,sol)<<std::endl;}
  pt[1]=c.m_cBox.y1;
  for(unsigned j=0;j<sol.size();j++)
    {
      pt[2]= c.m_cBox.z0*((coeff_t)1-sol[j]) + c.m_cBox.z1*sol[j];
      inters[1].push_back(pt);
      //      PRINT_DEBUG(pt<<" |"<<sol[j]<<"|  "<<c.m_cBox.x0<<" "<<c.m_cBox.x1<<" "<<c.m_cBox.y0<<" "<<c.m_cBox.y1<<" "<<c.m_cBox.z0<<" "<<c.m_cBox.z1);
      //      PRINT_DEBUG(pt);
      nbpt++;
    }

  face(lp[2].rep(), px1.rep(),0,0); 
  //  PRINT_DEBUG(" ");
  //  PRINT_DEBUG("X1Y0 "<<lp[2]);
  (*eqnList.begin())=lp[2];
  slv2.solve(sol, nSol, nVar, eqnList.begin(), eqnList.end());
  //  if(sol.size()) { VECTOR::print(std::cout,sol)<<std::endl;}

  pt[0]=c.m_cBox.x1;
  pt[1]=c.m_cBox.y0;
  for(unsigned j=0;j<sol.size();j++)
    {
      pt[2]= c.m_cBox.z0*((coeff_t)1-sol[j]) + c.m_cBox.z1*sol[j];
      inters[2].push_back(pt);
      //      PRINT_DEBUG(pt);
      nbpt++;
    }

  face(lp[2].rep(), px1.rep(),0,1); 
  //  PRINT_DEBUG(" ");
  //  PRINT_DEBUG("X1Y1 "<<lp[2]);
  (*eqnList.begin())=lp[2];
  slv2.solve(sol, nSol, nVar, eqnList.begin(), eqnList.end());
  //  if(sol.size()) { VECTOR::print(std::cout,sol)<<std::endl;}
  pt[1]=c.m_cBox.y1;
  for(unsigned j=0;j<sol.size();j++)
    {
      pt[2]= c.m_cBox.z0*((coeff_t)1-sol[j]) + c.m_cBox.z1*sol[j];
      inters[3].push_back(pt);
      //      PRINT_DEBUG(pt<<" |"<<sol[j]<<"|  "<<c.m_cBox.x0<<" "<<c.m_cBox.x1<<" "<<c.m_cBox.y0<<" "<<c.m_cBox.y1<<" "<<c.m_cBox.z0<<" "<<c.m_cBox.z1);
      //      PRINT_DEBUG(pt);
      nbpt++;
    }

  face(lp[1].rep(), px0.rep(),1,0); 
  SbdSolver<coeff_t,SBDSLV_RDL> slv1(lp+1, lp+2, 1e-4);

  //  PRINT_DEBUG(" ");
  //  PRINT_DEBUG("X0Z0 "<<lp[1]);
  (*eqnList.begin())=lp[1];
  slv1.solve(sol, nSol, nVar, eqnList.begin(), eqnList.end());
  //  if(sol.size()) { VECTOR::print(std::cout,sol)<<std::endl;}
  pt[0]=c.m_cBox.x0;
  pt[2]=c.m_cBox.z0;
  for(unsigned j=0;j<sol.size();j++)
    {
      pt[1]= c.m_cBox.y0*((coeff_t)1-sol[j]) + c.m_cBox.y1*sol[j];
      inters[4].push_back(pt);
      //      PRINT_DEBUG(pt<<" |"<<sol[j]<<"|  "<<c.m_cBox.x0<<" "<<c.m_cBox.x1<<" "<<c.m_cBox.y0<<" "<<c.m_cBox.y1<<" "<<c.m_cBox.z0<<" "<<c.m_cBox.z1);
      //      PRINT_DEBUG(pt);
      nbpt++;
    }
  
  face(lp[1].rep(), px0.rep(),1,1); 
  //  PRINT_DEBUG(" ");
  //  PRINT_DEBUG("X0Z1 "<<lp[1]);
  (*eqnList.begin())=lp[1];
  slv1.solve(sol, nSol, nVar, eqnList.begin(), eqnList.end());
  //  if(sol.size()) { VECTOR::print(std::cout,sol)<<std::endl;}
  pt[0]=c.m_cBox.x0;
  pt[2]=c.m_cBox.z1;
  for(unsigned j=0;j<sol.size();j++)
    {
      pt[1]= c.m_cBox.y0*((coeff_t)1-sol[j]) + c.m_cBox.y1*sol[j];
      inters[5].push_back(pt);
      //      PRINT_DEBUG(pt<<" |"<<sol[j]<<"|  "<<c.m_cBox.x0<<" "<<c.m_cBox.x1<<" "<<c.m_cBox.y0<<" "<<c.m_cBox.y1<<" "<<c.m_cBox.z0<<" "<<c.m_cBox.z1);
      //      PRINT_DEBUG(pt);
      nbpt++;
    }
  
  face(lp[1].rep(), px1.rep(),1,0); 
  //  PRINT_DEBUG(" ");
  //  PRINT_DEBUG("X1Z0 "<<lp[1]);
  (*eqnList.begin())=lp[1];
  slv1.solve(sol, nSol, nVar, eqnList.begin(), eqnList.end());
  //  if(sol.size()) { VECTOR::print(std::cout,sol)<<std::endl;}
  pt[0]=c.m_cBox.x1;
  pt[2]=c.m_cBox.z0;
  for(unsigned j=0;j<sol.size();j++)
    {
      pt[1]= c.m_cBox.y0*((coeff_t)1-sol[j]) + c.m_cBox.y1*sol[j];
      inters[6].push_back(pt);
      //      PRINT_DEBUG(pt<<" |"<<sol[j]<<"|  "<<c.m_cBox.x0<<" "<<c.m_cBox.x1<<" "<<c.m_cBox.y0<<" "<<c.m_cBox.y1<<" "<<c.m_cBox.z0<<" "<<c.m_cBox.z1);
      //      PRINT_DEBUG(pt);
      nbpt++;
    }
  
  face(lp[1].rep(), px1.rep(),1,1); 
  //  PRINT_DEBUG(" ");
  //  PRINT_DEBUG("X1Z1 "<<lp[1]);
  (*eqnList.begin())=lp[1];
  slv1.solve(sol, nSol, nVar, eqnList.begin(), eqnList.end());
  //  if(sol.size()) { VECTOR::print(std::cout,sol)<<std::endl;}
  pt[0]=c.m_cBox.x1;
  pt[2]=c.m_cBox.z1;
  for(unsigned j=0;j<sol.size();j++)
    {
      pt[1]= c.m_cBox.y0*((coeff_t)1-sol[j]) + c.m_cBox.y1*sol[j];
      inters[7].push_back(pt);
      //      PRINT_DEBUG(pt<<" |"<<sol[j]<<"|  "<<c.m_cBox.x0<<" "<<c.m_cBox.x1<<" "<<c.m_cBox.y0<<" "<<c.m_cBox.y1<<" "<<c.m_cBox.z0<<" "<<c.m_cBox.z1);
      //      PRINT_DEBUG(pt);
      nbpt++;
    }

  face(pz0.rep(), c.m_F[0].rep(),2,0); 
  face(pz1.rep(), c.m_F[0].rep(),2,1); 

  face(lp[0].rep(), pz0.rep(),1,0); 
  SbdSolver<coeff_t,SBDSLV_RDL> slv0(lp, lp+1, 1e-4);

  //  PRINT_DEBUG(" ");
  //  PRINT_DEBUG("Y0Z0 "<<lp[0]);
  (*eqnList.begin())=lp[0];
  slv0.solve(sol, nSol, nVar, eqnList.begin(), eqnList.end());
  //  if(sol.size()) { VECTOR::print(std::cout,sol)<<std::endl;}
  pt[1]=c.m_cBox.y0;
  pt[2]=c.m_cBox.z0;
  for(unsigned j=0;j<sol.size();j++)
    {
      pt[0]= c.m_cBox.x0*((coeff_t)1-sol[j]) + c.m_cBox.x1*sol[j];
      inters[8].push_back(pt);
      //      PRINT_DEBUG(pt<<" |"<<sol[j]<<"|  "<<c.m_cBox.x0<<" "<<c.m_cBox.x1<<" "<<c.m_cBox.y0<<" "<<c.m_cBox.y1<<" "<<c.m_cBox.z0<<" "<<c.m_cBox.z1);
      //      PRINT_DEBUG(pt);
      nbpt++;
    }
  
  face(lp[0].rep(), pz0.rep(),1,1); 
  //  PRINT_DEBUG(" ");
  //  PRINT_DEBUG("Y1Z0 "<<lp[0]);
   (*eqnList.begin())=lp[0];
  slv0.solve(sol, nSol, nVar, eqnList.begin(), eqnList.end());
  //  if(sol.size()) { VECTOR::print(std::cout,sol)<<std::endl;}
  pt[1]=c.m_cBox.y1;
  pt[2]=c.m_cBox.z0;
  for(unsigned j=0;j<sol.size();j++)
    {
      pt[0]= c.m_cBox.x0*((coeff_t)1-sol[j]) + c.m_cBox.x1*sol[j];
      inters[9].push_back(pt);
      //      PRINT_DEBUG(pt<<" |"<<sol[j]<<"|  "<<c.m_cBox.x0<<" "<<c.m_cBox.x1<<" "<<c.m_cBox.y0<<" "<<c.m_cBox.y1<<" "<<c.m_cBox.z0<<" "<<c.m_cBox.z1);
      //      PRINT_DEBUG(pt);
      nbpt++;
    }

  face(lp[0].rep(), pz1.rep(),1,0); 
  //  PRINT_DEBUG(" ");
  //  PRINT_DEBUG("Y0Z1 "<<lp[0]);
  (*eqnList.begin())=lp[0];
  slv0.solve(sol, nSol, nVar, eqnList.begin(), eqnList.end());
  //  if(sol.size()) { VECTOR::print(std::cout,sol)<<std::endl;}
  pt[1]=c.m_cBox.y0;
  pt[2]=c.m_cBox.z1;
  for(unsigned j=0;j<sol.size();j++)
    {
      pt[0]= c.m_cBox.x0*((coeff_t)1-sol[j]) + c.m_cBox.x1*sol[j];
      inters[10].push_back(pt);
      //      PRINT_DEBUG(pt<<" |"<<sol[j]<<"|  "<<c.m_cBox.x0<<" "<<c.m_cBox.x1<<" "<<c.m_cBox.y0<<" "<<c.m_cBox.y1<<" "<<c.m_cBox.z0<<" "<<c.m_cBox.z1);
      //      PRINT_DEBUG(pt);
      nbpt++;
    }

  face(lp[0].rep(), pz1.rep(),1,1); 
  //  PRINT_DEBUG(" ");
  //  PRINT_DEBUG("Y1Z1 "<<lp[0]);
  (*eqnList.begin())=lp[0];
  slv0.solve(sol, nSol, nVar, eqnList.begin(), eqnList.end());
  //  if(sol.size()) {VECTOR::print(std::cout,sol)<<std::endl;}
  pt[1]=c.m_cBox.y1;
  pt[2]=c.m_cBox.z1;
  for(unsigned j=0;j<sol.size();j++)
    {
      pt[0]= c.m_cBox.x0*((coeff_t)1-sol[j]) + c.m_cBox.x1*sol[j];
      inters[11].push_back(pt);
      //      PRINT_DEBUG(pt<<" |"<<sol[j]<<"|  "<<c.m_cBox.x0<<" "<<c.m_cBox.x1<<" "<<c.m_cBox.y0<<" "<<c.m_cBox.y1<<" "<<c.m_cBox.z0<<" "<<c.m_cBox.z1);
      //      PRINT_DEBUG(pt);
      nbpt++;
    }

  return nbpt;

}
//--------------------------------------------------------------------
template <class C, class CELL, class MSLV >
void TopSbdSurf3d<C, CELL, MSLV>::compute_points(node_t* n)
{
 


  std::vector<node_t*> Nbg[6];
  n->getAllBackNeighbors(Nbg[OCTREE::B]);
  n->getAllFrontNeighbors(Nbg[OCTREE::F]);
  n->getAllSouthNeighbors(Nbg[OCTREE::S]);
  n->getAllNorthNeighbors(Nbg[OCTREE::N]);
  n->getAllEastNeighbors(Nbg[OCTREE::E]);
  n->getAllWestNeighbors(Nbg[OCTREE::W]);
  cell_t* c=n->GetCell();
  //  PRINT_DEBUG(">>>>>["<<c->m_cBox.x0<<" "<<c->m_cBox.x1<<" "<<c->m_cBox.y0<<" "<<c->m_cBox.y1<<" "<<c->m_cBox.z0<<" "<<c->m_cBox.z1<<"]");
  //  PRINT_DEBUG(">>>>>|"<<c->m_cBox.x1-c->m_cBox.x0<<" "<<c->m_cBox.y1-c->m_cBox.y0<<" "<<c->m_cBox.z1-c->m_cBox.z0<<"|");
  unsigned f;
  for(unsigned i=0;i<6;i++)
    {
      f=i+(i%2?-1:1);
      //      PRINT_DEBUG(i<<" "<<f<<" sz "<<Nbg[i].size());
      c->m_FaceComputed[i]=false;
      for(unsigned j=0;j<Nbg[i].size();j++)
        {
          if(Nbg[i][j]->GetCell()->m_FaceComputed[f])
            {
              //            PRINT_DEBUG(i<<" "<<j<<" merge");
              c->mergelist(c->m_FaceLstPtIdx[i], Nbg[i][j]->GetCell()->m_FaceLstPtIdx[f]);
              c->m_FaceComputed[i]=true;
            }
          else
            {
              // There can be only one neightbor (of the same size), which faces are not computed
              // since the size of the cells is increasing in the list aNodeKept.
              assert(Nbg[i].size()==1);
              //                  cell_t* cf=Nbg[i][j]->GetCell();
              //                  PRINT_DEBUG("Face "<< f<<" of box "<<j<<" not computed");
              //                  PRINT_DEBUG("Contains variety? "<<is_variety_in(*Nbg[i][j]->GetCell()));
              //                  PRINT_DEBUG("["<<cf->m_cBox.x0<<" "<<cf->m_cBox.x1<<" "<<cf->m_cBox.y0<<" "<<cf->m_cBox.y1<<" "<<cf->m_cBox.z0<<" "<<cf->m_cBox.z1<<"]");
            }
          for(unsigned e=0;e<4;e++)
            if(is_edge_in_face(Nbg[i][j]->GetCell(),EdgesOfFace[f][e],c,i))
              c->mergelist(c->m_FaceLstPtIdx[i], Nbg[i][j]->GetCell()->m_EdgeLstPtIdx[EdgesOfFace[f][e]]);
        }
      if(!c->m_FaceComputed[i])
        {
          //  PRINT_DEBUG(i<<" compute points");
          if (c->m_nRegularIdx ==-1) 
            {
              computeinter_face(c->m_FaceLstPtIdx[i],*c,i/2,i%2);
              PRINT_DEBUG("End of face comp.");
            }
          c->m_FaceComputed[i]=true;
        }
    }
  for(unsigned j=0;j<12;j++)
    computeinter_edge(c->m_EdgeLstPtIdx[j],*c,OCTREE::EDGE_TYPE(j));

}
//--------------------------------------------------------------------
template <class C, class CELL, class MSLV >
template < class POINTGRAPH>
void TopSbdSurf3d<C, CELL, MSLV>::build_mesh(POINTGRAPH& g, node_t* n)
{
  typedef typename CELL::coeff_t         coeff_t;
  typedef typename CELL::MPOL_T          MPOL_T;
  typedef intersection3d<coeff_t>        intersection_t;
  typedef topology::point<coeff_t>       point_t;
  typedef shared_object<intersection_t*> intersection_ptr;

  cell_t* c= n->GetCell();

  // determine the regularity if it is not known
  if (c->m_nRegularIdx == -2)
    regularity(*c);
  
  if (c->m_nRegularIdx == -1) // a singular cell
    {
      build_box(g,c);
    }
  else // a regular cell
    {
      unsigned nbp=0;
      point_t m(3u);
      m[0]=0;m[1]=0;m[2]=0;
      for(unsigned i=0;i<6;i++)
        {
          nbp+=c->m_FaceLstPtIdx[i].size();
          for(unsigned j=0;j<c->m_FaceLstPtIdx[i].size();j++)
            {
              m[0]+=c->m_FaceLstPtIdx[i][j]->pnt[0];
              m[1]+=c->m_FaceLstPtIdx[i][j]->pnt[1];
              m[2]+=c->m_FaceLstPtIdx[i][j]->pnt[2];
            }
        }
      for(unsigned i=0;i<12;i++)
        {
          nbp+=c->m_EdgeLstPtIdx[i].size();
          for(unsigned j=0;j<c->m_EdgeLstPtIdx[i].size();j++)
            {
              m[0]+=c->m_EdgeLstPtIdx[i][j]->pnt[0];
              m[1]+=c->m_EdgeLstPtIdx[i][j]->pnt[1];
              m[2]+=c->m_EdgeLstPtIdx[i][j]->pnt[2];
            }
        }
      if(nbp)
        {
          m[0]/= nbp;
          m[1]/= nbp;
          m[2]/= nbp;
          //      Seq<intersection_t*> inters = c->m_aInterList;
          // intersections of the polar curve and the cell facets
          
          // Computeinter_edges(boundary,c);
          //      unsigned nbp = inters.size();
          // intersection of the surface with the edges of the box.
          //    PRINT_DEBUG(nbp);

          //        PRINT_DEBUG("Treating cell of reg "<<c->m_nRegularIdx<<" sz = "<<inters.size());
          //      for(unsigned i=0;i<inters.size();i++) PRINT_DEBUG(inters[i]->pnt);
          // build_box(g,c);
          
          //        PRINT_DEBUG("Ordering points on the border");
          std::vector<point_t> lf[6];
          unsigned nbp=0;
          for(unsigned i=0;i<6;i++)
            {
              for(unsigned j=0;j< c->m_FaceLstPtIdx[i].size();j++)
                lf[i].push_back(c->m_FaceLstPtIdx[i][j]->pnt);
              for(unsigned e=0;e<4;e++)
                for(unsigned j=0;j< c->m_EdgeLstPtIdx[EdgesOfFace[i][e]].size();j++)
                  {
                    lf[i].push_back(c->m_EdgeLstPtIdx[EdgesOfFace[i][e]][j]->pnt);
                    if(c->m_EdgeLstPtIdx[EdgesOfFace[i][e]][j]->nMulti==2)
                      lf[i].push_back(c->m_EdgeLstPtIdx[EdgesOfFace[i][e]][j]->pnt);
                  }
            }
//        PRINT_DEBUG(">>>>>["<<c->m_cBox.x0<<" "<<c->m_cBox.x1<<" "<<c->m_cBox.y0<<" "<<c->m_cBox.y1<<" "<<c->m_cBox.z0<<" "<<c->m_cBox.z1<<"]");
//        PRINT_DEBUG(">>>>>|"<<c->m_cBox.x1-c->m_cBox.x0<<" "<<c->m_cBox.y1-c->m_cBox.y0<<" "<<c->m_cBox.z1-c->m_cBox.z0<<"|");
//        PRINT_DEBUG("=======================");
//        for(unsigned i=0;i<6;i++) 
//          for(unsigned j=0;j< c->m_FaceLstPtIdx[i].size();j++)
//            PRINT_DEBUG("Face "<<i<<" "<<c->m_FaceLstPtIdx[i][j]->pnt);
//        PRINT_DEBUG(".......................");
//        for(unsigned i=0;i<12;i++) 
//          for(unsigned j=0;j< c->m_EdgeLstPtIdx[i].size();j++)
//            PRINT_DEBUG("Edge "<<i<<" "<<c->m_EdgeLstPtIdx[i][j]->pnt);
          
          //      PRINT_DEBUG("----------------------");
          //      PRINT_DEBUG("Computing point inside");
          compute_point_inside(m,c);

            //  VECTOR::print(std::cout,m)<<std::endl;
            //  int im=g.push_back(intersection_ptr(new intersection_t(m,0,1))),
            // sort according to the regular direction


          //        PRINT_DEBUG("Connecting the points by triangles");
          int i0, i1,
            im=g.push_back(m);
          //        PRINT_DEBUG(im<<" "<<m);
          // g.insert(im,iv);
          for(unsigned i=0;i<6;i++) 
            {
              if(lf[i].size()==2)
                {
                  i0=g.push_back(lf[i][0]);
                  i1=g.push_back(lf[i][1]);
                  g.insert_face(im,i0,i1);
                  g.insert(im,i0);
                  g.insert(im,i1);
//                PRINT_DEBUG("f "<<i<<" sz "<<lf[i].size());
//                for(unsigned l=0;l<lf[i].size();l++) PRINT_DEBUG(lf[i][l]);
                }
              else if(lf[i].size()>0)
                {
                  int reg= c->m_nRegularIdx %3;
                  int rf =reg;
                  //              if(rf>3) rf/=3;
                  if((i/2)==reg)
                    if(i%2)
                      rf+=c->m_nRegularIdx/9;
                    else
                      rf+=(c->m_nRegularIdx/3)%3;
                  else
                    rf+=(i/2);
                  //PRINT_DEBUG("face "<<i<<" sz "<<lf[i].size()<<" RIDX "<<c->m_nRegularIdx<<" reg "<<reg<<" rf "<<rf);
                  //for(unsigned l=0;l<lf[i].size();l++) PRINT_DEBUG(lf[i][l]);
                  switch(rf)
                    {
                    case 3:
                      //                      PRINT_DEBUG("Sort X");
                      sort(lf[i].begin(), lf[i].end(), XPT());
                      break;
                    case 2:
                      //                      PRINT_DEBUG("Sort Y");
                      sort(lf[i].begin(), lf[i].end(), YPT());
                      break;
                    case 1:
                      //                      PRINT_DEBUG("Sort Z");
                      sort(lf[i].begin(), lf[i].end(), ZPT());
                      break;
                    }
//                PRINT_DEBUG("f "<<i<<" sz "<<lf[i].size());
//                for(unsigned l=0;l<lf[i].size();l++) PRINT_DEBUG(lf[i][l]);
                  
                  int j=0, f=i/2;

                  i0=g.push_back(lf[i][0]);
                  g.insert(im,i0);
                  while(j<lf[i].size()-1)
                    {
                      if(same_point(lf[i][j],lf[i][j+1]))
                        j++;
                      else
                        {
                          i0=g.push_back(lf[i][j+1]);
                          g.insert(im,i0);
                          g.insert_face(im,i0-1,i0);
                          if(is_on_border(lf[i][j+1],c,f)) 
                            {
                              j+=2;
                              if(j<lf[i].size()) i0=g.push_back(lf[i][j]);
                              g.insert(im,i0);
                            }
                          else
                            j++;
                        }
                    }
                }
            }
        }
      
      //         // construct the graph
      //         bool bConnectWithLast = false;
      //         int nMulti = 1;
      //         int nLastIdx = -1;
      //        point_t ipt(3u);
      //         for (int v = 0; v < (int)inters.size(); v++)
      //         {
      //          // int nBIdx = g.push_back(intersection_ptr(inters[v]));
      //          ipt[0]=inters[v]->pnt[0];
      //          ipt[1]=inters[v]->pnt[1];
      //          ipt[2]=inters[v]->pnt[2];
      //          int nBIdx = g.push_back(ipt);
      
      //             if (bConnectWithLast)
      //                 g.insert(nLastIdx, nBIdx);
      
      //             if (!bConnectWithLast && nMulti < inters[v]->nMulti)
      //                 bConnectWithLast = false;
      //             else
      //                 bConnectWithLast = !bConnectWithLast;
      
      //             if (nMulti < inters[v]->nMulti)
      //             {
      //                 nMulti++;
      //                 v--;
      //             }
      //             else
      //             {
      //                 nMulti = 1;
      //             }
      
      //             nLastIdx = nBIdx;
      //         }

    }
}
//====================================================================
template<class C, class MP>
bool is_edge_in_face(Cell3d<C,MP>* c0, unsigned e, Cell3d<C,MP>* c, unsigned f)
{
  typedef topology::point<C>             point_t;
  using namespace OCTREE;
  point_t m(3u);
  switch(OCTREE::EDGE_TYPE(e))
    {
    case B_E:
      m[0]= c0->m_cBox.x0; m[1]= c0->m_cBox.y1; m[2]= (c0->m_cBox.z0 + c0->m_cBox.z1)/2; break;
    case B_W:
      m[0]= c0->m_cBox.x0; m[1]= c0->m_cBox.y0; m[2]= (c0->m_cBox.z0 + c0->m_cBox.z1)/2; break;
    case B_S:
      m[0]= c0->m_cBox.x0; m[1]= (c0->m_cBox.y0 + c0->m_cBox.y1)/2; m[2]= c0->m_cBox.z0; break;
    case B_N:
      m[0]= c0->m_cBox.x0; m[1]= (c0->m_cBox.y0 + c0->m_cBox.y1)/2; m[2]= c0->m_cBox.z1; break;

    case F_E:
      m[0]= c0->m_cBox.x1; m[1]= c0->m_cBox.y1; m[2]= (c0->m_cBox.z0 + c0->m_cBox.z1)/2; break;
    case F_W:
      m[0]= c0->m_cBox.x1; m[1]= c0->m_cBox.y0; m[2]= (c0->m_cBox.z0 + c0->m_cBox.z1)/2; break;
    case F_S:
      m[0]= c0->m_cBox.x1; m[1]= (c0->m_cBox.y0 + c0->m_cBox.y1)/2; m[2]= c0->m_cBox.z0; break;
    case F_N:
      m[0]= c0->m_cBox.x1; m[1]= (c0->m_cBox.y0 + c0->m_cBox.y1)/2; m[2]= c0->m_cBox.z1; break;

    case S_E:
      m[0]= (c0->m_cBox.x0 + c0->m_cBox.x1)/2; m[1]= c0->m_cBox.y1; m[2]= c0->m_cBox.z0; break;
    case S_W:
      m[0]= (c0->m_cBox.x0 + c0->m_cBox.x1)/2; m[1]= c0->m_cBox.y0; m[2]= c0->m_cBox.z0; break;
    case N_E:
      m[0]= (c0->m_cBox.x0 + c0->m_cBox.x1)/2; m[1]= c0->m_cBox.y1; m[2]= c0->m_cBox.z1; break;
    case N_W:
      m[0]= (c0->m_cBox.x0 + c0->m_cBox.x1)/2; m[1]= c0->m_cBox.y0; m[2]= c0->m_cBox.z1; break;
    }
  switch(f)
    {
    case 0:
      return (m[0]== c->m_cBox.x0 && m[1]>c->m_cBox.y0 && m[1]<c->m_cBox.y1 && m[2]>c->m_cBox.z0 && m[2]<c->m_cBox.z1);
    case 1:
      return (m[0]== c->m_cBox.x1 && m[1]>c->m_cBox.y0 && m[1]<c->m_cBox.y1 && m[2]>c->m_cBox.z0 && m[2]<c->m_cBox.z1);
    case 2:
      return (m[1]== c->m_cBox.y0 && m[0]>c->m_cBox.x0 && m[0]<c->m_cBox.x1 && m[2]>c->m_cBox.z0 && m[2]<c->m_cBox.z1);
    case 3:
      return (m[1]== c->m_cBox.y1 && m[0]>c->m_cBox.x0 && m[0]<c->m_cBox.x1 && m[2]>c->m_cBox.z0 && m[2]<c->m_cBox.z1);
    case 4:
      return (m[2]== c->m_cBox.z0 && m[0]>c->m_cBox.x0 && m[0]<c->m_cBox.x1 && m[1]>c->m_cBox.y0 && m[1]<c->m_cBox.y1);
    case 5:
      return (m[2]== c->m_cBox.z1 && m[0]>c->m_cBox.x0 && m[0]<c->m_cBox.x1 && m[1]>c->m_cBox.y0 && m[1]<c->m_cBox.y1);
    }
  return false;
}

//====================================================================
template < class C >
bool is_variety_in(Cell3d<C>& c)
{
  return c.containsVariety(0);
}
//--------------------------------------------------------------------
template <class C, class CELL, class MSLV >
bool TopSbdSurf3d<C, CELL, MSLV>::regularity(CELL& c)
{
  typedef typename CELL::MPOL_T mpol_t;
  if(c.m_nDepth<m_min_dp)
     return false;
  else if(c.m_nDepth>m_max_dp)
    {
      c.m_nRegularIdx=-1;
      return false;
    }
  mpol_t dp[3];
  mpol_t f0,f1,df0,df1;
  for (int i = 0; i < 3; i++)
    {
      diff(dp[i],c.m_F[0],VARIDX(c.m_F[0], i));
      //      PRINT_DEBUG(i);
      if(!has_sign_change(dp[i]))
        {
          //      PRINT_DEBUG("df/dx"<<i<<" no sign change");
          bool b0= true, b1=true;
          c.m_nRegularIdx = i;
          for(int j=0;j <3 && (b0 ||b1) ;j++)
            if(j!=i)
              {
                int v=j;
                if(j>i) v--;
                if(b0)
                  {
                    face(f0.rep(),c.m_F[0].rep(),i,0);
                    //              PRINT_DEBUG("f0 "<<v<<" "<<df0);
                    if(!has_sign_change(f0))
                      {
                        c.m_nRegularIdx += 3*j;
                        //PRINT_DEBUG("f0 "<<j<<" "<<c.m_nRegularIdx);
                        b0=false;
                        //                  return true;
                      }
                    else 
                      {
                        diff(df0, f0,  v);
                        if(!has_sign_change(df0))
                          {
                            c.m_nRegularIdx += 3*j;
                            //PRINT_DEBUG("f0 "<<j<<" "<<c.m_nRegularIdx);
                            b0=false;
                            //              return true;
                          }
                      }
                  }
                if(b1)
                  {
                    face(f1.rep(),c.m_F[0].rep(),i,1);
                    //              PRINT_DEBUG("f1 "<<v<<" "<<df1);
                    if (!has_sign_change(f1))
                      {
                        c.m_nRegularIdx += 9*j;
                        b1=false;
                        //PRINT_DEBUG("f1 "<<j<<" "<<c.m_nRegularIdx);
                      }
                    else
                      {
                        diff(df1, f1, v);
                        if (!has_sign_change(df1))
                          {
                            c.m_nRegularIdx += 9*j;
                            b1=false;
                            //PRINT_DEBUG("f1 "<<j<<" "<<c.m_nRegularIdx);
                          }
                      }
                  }
              }
          if(!(b0 || b1))
            {
              //PRINT_DEBUG("reg "<<  c.m_nRegularIdx);
              //PRINT_DEBUG("reg for x"<<i);
              return true;
            }
//        else
//          PRINT_DEBUG("reg not for x"<<i);
        }
    }
  c.m_nRegularIdx =-1;
  //  PRINT_DEBUG("reg "<<  c.m_nRegularIdx);
  return false;

  // check if the curve tangent, i.e. t = (t0, t1, t2) = grad F cross grad G
  // has no sign change for at least one of its component (or say, along one axis)
  for (unsigned i = c.m_nT; i < c.nbpol(); i++)
    {
      C* coeff = c.m_F[i].begin();
      C lastC = coeff[0];
      bool bSignChange = false;
      
      for (int j = 1; j < c.m_F[i].rep().env.sz() && !bSignChange; j++)
        {
          C cirC = coeff[j];
          if (_is_singular(cirC)) // if cirC is a zero value
            {
              continue;
            }
          else
            {
              if (lastC*cirC < 0) // a sign change is detected
                {
                  bSignChange = true;
                  lastC = cirC;
                  break;
                }
              
              lastC = cirC; // always remember the last non-singular coefficient
            }
        }
      
      if (_is_singular(lastC))  //all coefficients are singular (or 0)
        continue;
      else if (!bSignChange) // no sign changes
        {
            c.m_nRegularIdx = i;
            PRINT_DEBUG("REG "<<  c.m_nRegularIdx);
            return TRUE;
        }
    }
  c.m_nRegularIdx = -1;
  PRINT_DEBUG("reg "<<  c.m_nRegularIdx);
  return FALSE;
}
//====================================================================
namespace topology 
{
  template<class C, class POINTGRAPH, class MPOL, class MTH>
  void assign(POINTGRAPH& g, 
              const MPOL & p0,
              MTH mth,
              const C& x0, const C& x1, 
              const C& y0, const C& y1, 
              const C& z0, const C& z1)
  {
    typedef typename MTH::cell_t  cell_t;
    typedef typename MTH::mpolb_t mpolb_t;

    std::cout<<p0<<std::endl;
    
    // step1:
    // computing grad F1 cross grad F2
    MPOL p1 = diff(p0,0u);
    MPOL p2 = diff(p0,1u);

    MPOL d0 = diff(p1,1u)*diff(p2,2u)-diff(p1,2u)*diff(p2,1u);
    MPOL d1 = diff(p1,0u)*diff(p2,2u)-diff(p1,2u)*diff(p2,0u);
    MPOL d2 = diff(p1,0u)*diff(p2,1u)-diff(p1,1u)*diff(p2,0u);
    std::cout<<"Tangent field:"<<std::endl;

    std::cout<<d0<<"\n"<<d1<<"\n"<<d2<<std::endl;
    mpolb_t p0dse;

    toBernstein(p0, p0dse, x0,x1,y0,y1,z0,z1);
//     mpolb_t  p1dse, p2dse;
//     toBernstein(p1, p1dse, x0,x1,y0,y1,z0,z1);
//     toBernstein(p2, p2dse, x0,x1,y0,y1,z0,z1);

//     // debug only
//     std::cout << "F("
//               << degree(p1,0) << ", "
//               << degree(p1,1) << ", "
//               << degree(p1,2) << "): "
//               << p1 << std::endl;

//     std::cout << "G("
//               << degree(p2,0) << ", "
//               << degree(p2,1) << ", "
//               << degree(p2,2) << "): "
//               << p2 << std::endl;

//     std::cout << "Fb("
//               << p1dse.rep().env.sz(VARIDX(pd1dse, 0)) << ", "
//               << p1dse.rep().env.sz(VARIDX(pd1dse, 1)) << ", "
//               << p1dse.rep().env.sz(VARIDX(pd1dse, 2)) << "): "
//               << p1dse << std::endl;

//     std::cout << "Gb("
//               << p2dse.rep().env.sz(VARIDX(pd2dse, 0)) << ", "
//               << p2dse.rep().env.sz(VARIDX(pd2dse, 1)) << ", "
//               << p2dse.rep().env.sz(VARIDX(pd2dse, 2)) << "): "
//               << p2dse << std::endl;


//     diff(p1_dy, p1dse, VARIDX(pd1dse, 1));
//     diff(p1_dz, p1dse, VARIDX(pd1dse, 2));
//     diff(p2_dx, p2dse, VARIDX(pd2dse, 0));
//     diff(p2_dy, p2dse, VARIDX(pd2dse, 1));
//     diff(p2_dz, p2dse, VARIDX(pd2dse, 2));

//     std::cout << "grad Fb = {" << p1_dx << ", " <<  std::endl
//               << "           " << p1_dy << ", " <<  std::endl
//               << "           " << p1_dz << "} " <<  std::endl;
//     std::cout << "grad Gb = {" << p2_dx << ", " <<  std::endl
//               << "           " << p2_dy << ", " <<  std::endl
//               << "           " << p2_dz << "} " <<  std::endl;


/*   
     mpolb_t c0,c1,c2;
     toBernstein(d0, c0, x0,x1,y0,y1,z0,z1);
     toBernstein(d1, c1, x0,x1,y0,y1,z0,z1);
     toBernstein(d2, c2, x0,x1,y0,y1,z0,z1);
     std::cout<<c0<<"\n"<<c1<<"\n"<<c2<<std::endl;
*/
//     mpolb_t c1, c2, c3; // grad F1 cross grad F2
//     using arithm::operator*;
//     c1 = p1_dy * p2_dz - p1_dz * p2_dy;
//     c2 = p2_dx * p1_dz - p2_dz * p1_dx;
//     c3 = p1_dx * p2_dy - p1_dy * p2_dx;

//     // debug only
//     std::cout << "grad(F) X grad(G) = {" << std::endl
//               << "                     " << c1 << "," << std::endl
//               << "                     " << c2 << "," << std::endl
//               << "                     " << c3 << "}" << std::endl;

    // step2:
    // assign to the octree
    cell_t* rc = new cell_t(p0dse,x0,x1,y0,y1,z0,z1);

    //    rc->m_F =std::vector<mpolb_t>(7);
    //    rc->m_F[0] = p0dse;
    //    rc->m_F[1] = p0_dx;
    //    rc->m_F[2] = p0_dy;
    //    rc->m_F[3] = p0_dz;
    //    rc->m_F[4] = c1;
    //    rc->m_F[5] = c2;
    //    rc->m_F[6] = c3;
    //    rc->m_F[1] = c1;
    //    rc->m_F[2] = c2;
    //    rc->m_F[3] = c3;
    mth.run(g,rc);
  }
}// namespace topology
__END_NAMESPACE_SYNAPS
/********************************************************************/
#endif //CELL3DCURVE_H