shape: include/shape/bcell3d.hpp Source File

00001 /*****************************************************************************
00002  * M a t h e m a  g i x
00003  *****************************************************************************
00004  * Cell
00005  * 2008-03-20
00006  * Julien Wintz & Bernard Mourrain & Angelos Mantzaflaris
00007  *****************************************************************************
00008  *               Copyright (C) 2008 INRIA Sophia-Antipolis
00009  *****************************************************************************
00010  * Comments :
00011  ****************************************************************************/
00012 # ifndef shape_bcell3d_hpp
00013 # define shape_bcell3d_hpp
00014 
00015 # include <shape/vertex.hpp>
00016 # include <shape/edge.hpp>
00017 # include <shape/face.hpp>
00018 # include <shape/bounding_box.hpp>
00019 # include <shape/cell3d.hpp>
00020 # include <shape/graph.hpp>
00021 # include <shape/topology.hpp>
00022 # include <shape/solver_implicit.hpp>
00023 //# include <shape/tpl3d.hpp>
00024 # define TMPL template<class C, class V>
00025 # define TMPL1 template<class K>
00026 # define SELF bcell3d<C,V>
00027 //====================================================================
00028 namespace mmx { namespace shape {
00029 
00030     // TMPL struct bcell3d;
00031     // TMPL struct tpl3d;
00032 
00033 TMPL class bcell3d;
00034 TMPL class tpl3d;
00035 
00036 //--------------------------------------------------------------------
00037 struct bcell3d_def {};
00038 template<> struct use<bcell3d_def>
00039   :public use<bcell_def>
00040 {
00041   typedef bcell3d<double,double> Cell3d;
00042 };
00043 
00044 //--------------------------------------------------------------------
00045 template<class C, class V=default_env>
00046 class bcell3d : public  cell3d<C,V>
00047 {
00048 public:
00049   typedef topology<C,REF_OF(V)>         Topology;  
00050   typedef tpl3d<C,V>                    Topology3d;  
00051   typedef typename Topology::Point      Point;  
00052   typedef typename Topology::Edge       Edge;
00053 
00054   typedef cell3d<C,V>                   CellBase;
00055   typedef bcell<C,V>                    Cell;
00056   typedef typename Cell::BoundingBox    BoundingBox;
00057 
00058   typedef solver_implicit<C,V> Solver;
00059 
00060   bcell3d(void) {};
00061   bcell3d(const BoundingBox& bx): CellBase(bx) {} ;
00062   virtual ~bcell3d(void) {};
00063 
00064   virtual bool        is_regular(void) = 0 ;
00065   virtual bool        is_active (void) const = 0 ;
00066   virtual bool        is_intersected(void)   = 0 ;
00067 //  virtual int         subdivide(SELF*& left, SELF*& right);
00068 //  virtual void        subdivide(CellBase*& left, CellBase*& right, int v, double s);
00069 
00070   virtual void        split_position(int &v, double& s) ;
00071   virtual bool        insert_regular (Topology*) = 0 ;
00072   virtual bool        insert_singular(Topology*) = 0 ;
00073   
00074   virtual Point       center(void) const { 
00075     // XXX ???
00076     return Point(this->xmax()-this->xmin(),
00077                  this->ymax()-this->ymin(),
00078                  this->zmax()-this->zmin()); 
00079   }
00080 
00081   BoundingBox boundingBox() const { return (BoundingBox)*this; }
00082 
00083 
00084   virtual bool topology_regular (Topology*) = 0 ;
00085   virtual void polygonise (Topology3d*) =0;
00086 
00087   inline void disconnect();
00088   inline void connect0(SELF * a, SELF * b);
00089   inline void connect1(SELF * a, SELF * b) ;
00090   inline void connect2(SELF * a, SELF * b) ;
00091   inline void join0( SELF * b) ;
00092   inline void join1( SELF * b) ;
00093   inline void join2( SELF * b) ;
00094 
00095   Seq<Point *>  m_boundary;
00096   Seq<Point *>  m_singular ;
00097   int           m_type;
00098 
00099   /*Pointer to graph instance of this bcell*/
00100   gNode<bcell3d*>* m_gnode;
00101 
00102   /*Neighbor bcells*/
00103   Seq<bcell3d *>  neighbors(){ 
00104     Seq<bcell3d *> t;
00105     t<< s_neighbors ;
00106     t<< e_neighbors ;
00107     t<< n_neighbors ;
00108     t<< b_neighbors ;
00109     t<< w_neighbors ;
00110     t<< f_neighbors ;
00111     return t;     } ;
00112 
00113   Seq<bcell3d *>  s_neighbors ;
00114   Seq<bcell3d *>  e_neighbors ;
00115   Seq<bcell3d *>  n_neighbors ;
00116   Seq<bcell3d *>  w_neighbors ;
00117   Seq<bcell3d *>  f_neighbors ;
00118   Seq<bcell3d *>  b_neighbors ;
00119 
00120   inline Seq<Point *> intersections(int i) const;  
00121 
00122 };
00123 
00124 //====================================================================
00125 //TMPL int
00126 //SELF::subdivide(SELF*& Left, SELF*& Right) {
00127 //  int v; double s;
00128 //  this->split_position(v,s);
00129 //  this->subdivide(Left,Right,v,s);
00130 //  return v;
00131 //}
00132 //====================================================================
00133 TMPL void 
00134 SELF::split_position(int& v, double& s) {
00135   double sx = (this->xmax()-this->xmin());
00136   double sy = (this->ymax()-this->ymin());
00137   double sz = (this->zmax()-this->zmin());
00138   
00139   if(sx<sy)
00140       if(sy<sz) {
00141         v=2;
00142         s=(this->zmax()+this->zmin())/2;
00143       } else {
00144         v=1;
00145         s=(this->ymax()+this->ymin())/2;
00146       }
00147     else
00148       if(sx<sz) {
00149         v=2;
00150         s=(this->zmax()+this->zmin())/2;
00151       } else {
00152         v=0;
00153         s=(this->xmax()+this->xmin())/2;
00154       }
00155 }
00156 
00157 TMPL inline void 
00158 SELF::join0(SELF * b) //left,right
00159 {
00160     this->e_neighbors << b; 
00161     b->w_neighbors << this; 
00162 }
00163     
00164 TMPL inline void 
00165 SELF::join1(SELF * b) //up,down
00166 {
00167     b->s_neighbors << this;
00168     this->n_neighbors << b;
00169 }
00170 
00171 TMPL inline void 
00172 SELF::join2(SELF * b) //front,back
00173 {
00174     b->f_neighbors << this;
00175     this->b_neighbors << b;
00176 }
00177 
00178 TMPL inline void 
00179 SELF::connect0(SELF * a, SELF * b) 
00180 {
00181     int i;
00182     bool flag;
00183 
00184     //copy horizontally
00185     b->e_neighbors= this->e_neighbors ;
00186     foreach(SELF* cl,b->e_neighbors) {
00187       i= cl->w_neighbors.search(this);
00188       cl->w_neighbors[i]= b;
00189     }
00190 
00191     a->w_neighbors= this->w_neighbors ;
00192     foreach(SELF* cl,a->w_neighbors) {
00193       i= cl->e_neighbors.search(this);
00194       cl->e_neighbors[i]= a;
00195     }
00196 
00197     //update vertically
00198     foreach(SELF* cl,this->s_neighbors) {
00199       flag=false;
00200       if ( check_overlap(cl,a,0) //)
00201            && ( check_overlap(cl,a,1) || check_overlap(cl,a,2)) )
00202         {
00203           //assert( cl->ymax()== a->ymin() );
00204           a->s_neighbors<< cl;
00205           i= cl->n_neighbors.search(this);
00206           cl->n_neighbors[i]= a;
00207           flag=true;
00208         }
00209       if ( check_overlap(cl,b,0) //)
00210            && ( check_overlap(cl,a,1) || check_overlap(cl,a,2)) )
00211         {
00212           //assert( cl->ymax()== b->ymin() );
00213           b->s_neighbors<< cl;
00214           if (!flag)
00215             {
00216               i= cl->n_neighbors.search(this);
00217               cl->n_neighbors[i]= b;
00218             }
00219           else
00220             cl->n_neighbors << b;
00221         }
00222     }
00223     foreach(SELF* cl,this->n_neighbors) {
00224         flag=false;
00225         if ( check_overlap(cl,a,0) //)
00226              && ( check_overlap(cl,a,1) || check_overlap(cl,a,2)) ) 
00227         {
00228             a->n_neighbors<< cl;
00229             i= cl->s_neighbors.search(this);
00230             cl->s_neighbors[i]= a;
00231             flag=true;
00232         }
00233         if ( check_overlap(cl,b,0) //)
00234              && ( check_overlap(cl,a,1) || check_overlap(cl,a,2)) )
00235         {
00236             b->n_neighbors<< cl;
00237             if (!flag)
00238             {
00239               i= cl->s_neighbors.search(this);
00240               cl->s_neighbors[i]= b;
00241             }
00242             else
00243               cl->s_neighbors << b;
00244         }
00245     } 
00246 
00247     //update depth
00248     foreach(SELF* cl,this->f_neighbors) {
00249       flag=false;
00250       if ( check_overlap(cl,a,0) //)
00251            && ( check_overlap(cl,a,1) || check_overlap(cl,a,2)) )
00252         {
00253           //assert( cl->ymax()== a->ymin() );
00254           a->f_neighbors<< cl;
00255           i= cl->b_neighbors.search(this);
00256           cl->b_neighbors[i]= a;
00257           flag=true;
00258         }
00259       if ( check_overlap(cl,b,0) //)
00260            && ( check_overlap(cl,a,1) || check_overlap(cl,a,2)) )
00261         {
00262           //assert( cl->ymax()== b->ymin() );
00263           b->f_neighbors<< cl;
00264           if (!flag)
00265             {
00266               i= cl->b_neighbors.search(this);
00267               cl->b_neighbors[i]= b;
00268             }
00269           else
00270             cl->b_neighbors << b;
00271         }
00272     }
00273     foreach(SELF* cl,this->b_neighbors) {
00274         flag=false;
00275         if ( check_overlap(cl,a,0) //)
00276            && ( check_overlap(cl,a,1) || check_overlap(cl,a,2)) ) 
00277         {
00278             a->b_neighbors<< cl;
00279             i= cl->f_neighbors.search(this);
00280             cl->f_neighbors[i]= a;
00281             flag=true;
00282         }
00283         if ( check_overlap(cl,b,0) //)
00284            && ( check_overlap(cl,a,1) || check_overlap(cl,a,2)) )
00285         {
00286             b->b_neighbors<< cl;
00287             if (!flag)
00288             {
00289               i= cl->f_neighbors.search(this);
00290               cl->f_neighbors[i]= b;
00291             }
00292             else
00293               cl->f_neighbors << b;
00294         }
00295     }
00296 }
00297 
00298 TMPL inline void 
00299 SELF::connect1(SELF * a, SELF * b) 
00300 {
00301     int i;
00302     bool flag;
00303 
00304     //copy vertically
00305     a->s_neighbors= this->s_neighbors ;
00306     foreach(SELF* cl,a->s_neighbors) {
00307       i= cl->n_neighbors.search(this);
00308       cl->n_neighbors[i]= a;
00309     }
00310     b->n_neighbors= this->n_neighbors ;
00311     foreach(SELF* cl,b->n_neighbors) {
00312       i= cl->s_neighbors.search(this);
00313       cl->s_neighbors[i]= b;
00314     }
00315     
00316     //update horizontally
00317     foreach(SELF* cl,this->w_neighbors) {
00318       flag=false;
00319       if ( check_overlap(cl,a,1) //) 
00320         && ( check_overlap(cl,a,0) || check_overlap(cl,a,2)) ) 
00321         {
00322           //assert( cl->xmax()== a->xmin() );
00323           a->w_neighbors<< cl;
00324           i= cl->e_neighbors.search(this);
00325           cl->e_neighbors[i]= a;
00326           flag=true;
00327         }
00328       if ( check_overlap(cl,b,1) //)
00329         && ( check_overlap(cl,a,0) || check_overlap(cl,a,2)) ) 
00330         {
00331           //assert( cl->xmax()== b->xmin() );
00332           b->w_neighbors<< cl;
00333           if (!flag)
00334             {
00335               i= cl->e_neighbors.search(this);
00336               cl->e_neighbors[i]= b;
00337             }
00338           else
00339             cl->e_neighbors << b;
00340         }
00341     }
00342     foreach(SELF* cl,this->e_neighbors) {
00343       flag=false;
00344       if ( check_overlap(cl,a,1) //)
00345         && ( check_overlap(cl,a,0) || check_overlap(cl,a,2)) )  
00346         {
00347           a->e_neighbors<< cl;
00348           i= cl->w_neighbors.search(this);
00349           cl->w_neighbors[i]= a;
00350           flag=true;
00351         }
00352       if ( check_overlap(cl,b,1) //)
00353         && ( check_overlap(cl,a,0) || check_overlap(cl,a,2)) ) 
00354         {
00355           b->e_neighbors<< cl;
00356           if (!flag)
00357             {
00358               i= cl->w_neighbors.search(this);
00359               cl->w_neighbors[i]= b;
00360             }
00361           else
00362             cl->w_neighbors << b;
00363         }
00364     }
00365 
00366 
00367     //update depth
00368     foreach(SELF* cl,this->f_neighbors) {
00369       flag=false;
00370       if ( check_overlap(cl,a,1) //)
00371         && ( check_overlap(cl,a,0) || check_overlap(cl,a,2)) )  
00372         {
00373           //assert( cl->xmax()== a->xmin() );
00374           a->f_neighbors<< cl;
00375           i= cl->b_neighbors.search(this);
00376           cl->b_neighbors[i]= a;
00377           flag=true;
00378         }
00379       if ( check_overlap(cl,b,1) //)
00380         && ( check_overlap(cl,a,0) || check_overlap(cl,a,2)) ) 
00381         {
00382           //assert( cl->xmax()== b->xmin() );
00383           b->f_neighbors<< cl;
00384           if (!flag)
00385             {
00386               i= cl->b_neighbors.search(this);
00387               cl->b_neighbors[i]= b;
00388             }
00389           else
00390             cl->b_neighbors << b;
00391         }
00392     }
00393     foreach(SELF* cl,this->b_neighbors) {
00394       flag=false;
00395       if ( check_overlap(cl,a,1) //)
00396         && ( check_overlap(cl,a,0) || check_overlap(cl,a,2)) )  
00397         {
00398           a->b_neighbors<< cl;
00399           i= cl->f_neighbors.search(this);
00400           cl->f_neighbors[i]= a;
00401           flag=true;
00402         }
00403       if ( check_overlap(cl,b,1) //)
00404         && ( check_overlap(cl,a,0) || check_overlap(cl,a,2)) ) 
00405         {
00406           b->b_neighbors<< cl;
00407           if (!flag)
00408             {
00409               i= cl->f_neighbors.search(this);
00410               cl->f_neighbors[i]= b;
00411             }
00412           else
00413             cl->f_neighbors << b;
00414         }
00415     }
00416 }
00417 
00418 //--------------------------------------------------------------------
00419 TMPL inline void 
00420 SELF::connect2(SELF * a, SELF * b) 
00421 {
00422     int i;
00423     bool flag;
00424 
00425     //copy vertically
00426     a->f_neighbors= this->f_neighbors ;
00427     foreach(SELF* cl,a->f_neighbors) {
00428       i= cl->b_neighbors.search(this);
00429       cl->b_neighbors[i]= a;
00430     }
00431     b->b_neighbors= this->b_neighbors ;
00432     foreach(SELF* cl,b->b_neighbors) {
00433       i= cl->f_neighbors.search(this);
00434       cl->f_neighbors[i]= b;
00435     }
00436     
00437     //update horizontally
00438     foreach(SELF* cl,this->w_neighbors) {
00439       flag=false;
00440       if ( check_overlap(cl,a,2) //)
00441         && ( check_overlap(cl,a,0) || check_overlap(cl,a,1)) )  
00442         {
00443           //assert( cl->xmax()== a->xmin() );
00444           a->w_neighbors<< cl;
00445           i= cl->e_neighbors.search(this);
00446           cl->e_neighbors[i]= a;
00447           flag=true;
00448         }
00449       if ( check_overlap(cl,b,2) //)
00450         && ( check_overlap(cl,a,0) || check_overlap(cl,a,1)) )  
00451         {
00452           //assert( cl->xmax()== b->xmin() );
00453           b->w_neighbors<< cl;
00454           if (!flag)
00455             {
00456               i= cl->e_neighbors.search(this);
00457               cl->e_neighbors[i]= b;
00458             }
00459           else
00460             cl->e_neighbors << b;
00461         }
00462     }
00463     foreach(SELF* cl,this->e_neighbors) {
00464       flag=false;
00465       if ( check_overlap(cl,a,2) //)
00466         && ( check_overlap(cl,a,0) || check_overlap(cl,a,1)) )   
00467         {
00468           a->e_neighbors<< cl;
00469           i= cl->w_neighbors.search(this);
00470           cl->w_neighbors[i]= a;
00471           flag=true;
00472         }
00473       if ( check_overlap(cl,b,2) //)
00474         && ( check_overlap(cl,a,0) || check_overlap(cl,a,1)) )  
00475         {
00476           b->e_neighbors<< cl;
00477           if (!flag)
00478             {
00479               i= cl->w_neighbors.search(this);
00480               cl->w_neighbors[i]= b;
00481             }
00482           else
00483             cl->w_neighbors << b;
00484         }
00485     }
00486     //update vertically
00487     foreach(SELF* cl,this->s_neighbors) {
00488       flag=false;
00489       if ( check_overlap(cl,a,0) //)
00490         && ( check_overlap(cl,a,0) || check_overlap(cl,a,1)) )  
00491         {
00492           //assert( cl->ymax()== a->ymin() );
00493           a->s_neighbors<< cl;
00494           i= cl->n_neighbors.search(this);
00495           cl->n_neighbors[i]= a;
00496           flag=true;
00497         }
00498       if ( check_overlap(cl,b,0) //)
00499         && ( check_overlap(cl,a,0) || check_overlap(cl,a,1)) )  
00500         {
00501           //assert( cl->ymax()== b->ymin() );
00502           b->s_neighbors<< cl;
00503           if (!flag)
00504             {
00505               i= cl->n_neighbors.search(this);
00506               cl->n_neighbors[i]= b;
00507             }
00508           else
00509             cl->n_neighbors << b;
00510         }
00511     }
00512     foreach(SELF* cl,this->n_neighbors) {
00513         flag=false;
00514         if ( check_overlap(cl,a,2) //)
00515         && ( check_overlap(cl,a,0) || check_overlap(cl,a,1)) )  
00516         {
00517           a->n_neighbors<< cl;
00518           i= cl->s_neighbors.search(this);
00519           cl->s_neighbors[i]= a;
00520           flag=true;
00521         }
00522         if ( check_overlap(cl,b,2) //)
00523         && ( check_overlap(cl,a,0) || check_overlap(cl,a,1)) )  
00524           {
00525             b->n_neighbors<< cl;
00526             if (!flag) {
00527               i= cl->s_neighbors.search(this);
00528               cl->s_neighbors[i]= b;
00529             }
00530             else
00531               cl->s_neighbors << b;
00532         }
00533     }
00534 }
00535 //--------------------------------------------------------------------
00536 TMPL inline void 
00537 SELF::disconnect()
00538 {
00539     this->e_neighbors.clear();  
00540     this->w_neighbors.clear();  
00541     this->n_neighbors.clear();  
00542     this->s_neighbors.clear();  
00543     this->b_neighbors.clear();  
00544     this->f_neighbors.clear();  
00545 }
00546 //--------------------------------------------------------------------
00547 TMPL
00548 inline bool check_overlap(SELF * a, SELF * b, int v) //left,right or up,down
00549 {// check if a,b overlap wrt v-th coordinate
00550     if (v==0)  //direction v=0
00551       return !( a->xmax()<= b->xmin() ||
00552                 b->xmax()<= a->xmin() );
00553     if (v==1)  //direction v=1
00554       return !( a->ymax()<= b->ymin() ||
00555                 b->ymax()<= a->ymin() );
00556     if (v==2)  //direction v=2
00557       return !( a->zmax()<= b->zmin() ||
00558                 b->zmax()<= a->zmin() );
00559     return false;
00560 }
00561 //--------------------------------------------------------------------
00562 TMPL 
00563 inline Seq<typename SELF::Point *>
00564 SELF::intersections(int i) const
00565       {
00566         Seq<Point *> l;
00567         switch(i) {
00568         case 0://south
00569           foreach(Point* p, this->m_boundary) 
00570             if (Solver::south_face.is_valid(*p,this->boundingBox())) 
00571               l << p ;
00572             return l;
00573           case 1://east
00574           foreach(Point* p, this->m_boundary) 
00575             if (Solver::east_face.is_valid(*p,this->boundingBox())) 
00576               l << p ;
00577             return l;
00578           case 2://north
00579           foreach(Point* p, this->m_boundary) 
00580             if (Solver::north_face.is_valid(*p,this->boundingBox())) 
00581               l << p ;
00582             return l;
00583           case 3://west
00584           foreach(Point* p, this->m_boundary) 
00585             if (Solver::west_face.is_valid(*p,this->boundingBox())) 
00586               l << p ;
00587             return l;
00588           case 4://front
00589             foreach(Point* p, this->m_boundary) 
00590               if (Solver::front_face.is_valid(*p,this->boundingBox())) 
00591                 l << p ;
00592               return l;
00593             case 5://back
00594           foreach(Point* p, this->m_boundary) 
00595             if (Solver::back_face.is_valid(*p,this->boundingBox())) 
00596               l << p ;
00597             return l;
00598           default:
00599             return (Seq<Point *>());
00600                     }
00601           }
00602 
00603   //--------------------------------------------------------------------
00604 template<class CELL> inline void 
00605 bcell3d_split(CELL * left, CELL * right, CELL* cl, int v, double s) {
00606  if(v==0) {
00607    left = new CELL(*cl);
00608    right = new CELL(*cl);
00609    left->set_xmax(s);
00610    right->set_xmin(s);
00611    cl->connect0(left,right);
00612    left->join0(right);
00613  } else if (v==1) {
00614    left = new CELL(*cl);
00615    right = new CELL(*cl);
00616    left->set_ymax(s);
00617    right->set_ymin(s);
00618    cl->connect1(left,right);
00619    left->join2(right);
00620  } else if (v==2) {
00621    left = new CELL(*cl);
00622    right = new CELL(*cl);
00623    left->set_zmax(s);
00624    right->set_zmin(s);
00625    cl->connect2(left,right);
00626    left->join2(right);
00627  }
00628 }
00629 //--------------------------------------------------------------------
00630 template<class CELL> bool 
00631 is_adjacent_cell3d(CELL* c1,CELL* c2){
00632   if(c1->xmax()<c2->xmin() || c2->xmax()<c1->xmin())
00633     return false;
00634   if(c1->ymax()<c2->ymin() || c2->ymax()<c1->ymin())
00635     return false;
00636   if(c1->zmax()<c2->zmin() || c2->zmax()<c1->zmin())
00637     return false;
00638   if((c1->xmax()==c2->xmin() || c2->xmax()==c1->xmin())) {
00639     if((c1->ymax()==c2->ymin() || c2->ymax()==c1->ymin()) ||
00640        (c1->zmax()==c2->zmin() || c2->zmax()==c1->zmin()) )
00641       return false;
00642   } else if((c1->ymax()==c2->ymin() || c2->ymax()==c1->ymin()) &&
00643             (c1->zmax()==c2->zmin() || c2->zmax()==c1->zmin()) )
00644     return false;
00645   return true;
00646 }
00647 
00648 //====================================================================
00649 } ; // namespace shape
00650 } ; // namespace mmx
00651 # undef TMPL
00652 # undef TMPL1
00653 # undef SELF
00654 # endif // shape_cell3d_hpp