00001
00002
00003
00004
00005
00006 #ifndef synaps_topology_sbd2d_H
00007 #define synaps_topology_sbd2d_H
00008
00009 #include <synaps/init.h>
00010 #include <synaps/arithm/sign.h>
00011 #include <synaps/upol/UPolDse.h>
00012 #include <synaps/upol/bezier/BEZIER.m>
00013 #include <synaps/linalg/VectDse.h>
00014 #include <synaps/linalg/MatrDse.h>
00015 #include <synaps/upol/bezier/rep1d.h>
00016 #include <synaps/usolve/bezier/sbd1d.h>
00017 #include <synaps/arithm/sign.h>
00018 #include <synaps/topology/cell2d.h>
00019 #include <synaps/topology/point_graph.h>
00020
00021 __BEGIN_NAMESPACE_SYNAPS
00022
00023 struct YX {
00024 template<class T>
00025 static bool less(const T& a, const T& b)
00026 {
00027 return (a[1]<b[1] || (a[1]==b[1] && a[0]< b[0]));
00028 }
00029 template<class T>
00030 bool operator()(const T& a, const T& b)
00031 {
00032 return less(a,b);
00033 }
00034 };
00035
00036
00037 struct XY {
00038 template<class T>
00039 static bool less(const T& a, const T& b)
00040 {
00041 return (a[0]<b[0] || (a[0]==b[0] && a[1]< b[1]));
00042 }
00043 template<class T>
00044 bool operator()(const T& a, const T& b)
00045 {
00046 return less(a,b);
00047 }
00048 };
00049
00050 namespace BEZIER {
00051
00052
00053 template<class MATR, class MPOL, class X>
00054 void monomial2bezier2d(MATR& M, const MPOL& p,
00055 const X& x0, const X& x1, const X& y0, const X& y1)
00056 {
00057 using let::assign;
00058 typedef typename MATR::coeff_t coeff_t;
00059 M.resize(degree(p,1)+1, degree(p,0)+1);
00060 for(typename MPOL::const_iterator it=p.begin(); it != p.end(); it++)
00061 {
00062 assign(M((*it)[1],(*it)[0]),it->coeff());
00063 }
00064 for(unsigned i=0;i<M.nbrow();i++)
00065 {
00066 VectDse<coeff_t> V = row(M,i), BZ(M.nbcol(),AsSize());
00067 BEZIER::monomials2bezier(BZ,V,V.size(),x0,x1);
00068
00069 for(unsigned j=0;j<M.nbcol();j++) M(i,j)=BZ[j];
00070 }
00071
00072 for(unsigned j=0;j<M.nbcol();j++)
00073 {
00074 VectDse<coeff_t> V= col(M,j), BZ(M.nbrow(),AsSize());
00075 BEZIER::monomials2bezier(BZ,V,V.size(),y0,y1);
00076
00077 for(unsigned i=0;i<M.nbrow();i++) M(i,j)=BZ[i];
00078 }
00079 }
00080
00081 template<class matr_t>
00082 void deCasteljau( matr_t& left, matr_t& right, ByRow M )
00083 {
00084 typedef typename matr_t::value_type coeff_t;
00085
00086 int col=left.nbcol(),row=left.nbrow();
00087 assert( col > 1 );
00088 coeff_t middle = coeff_t(0.5);
00089 for ( int i = 0 ; i < row ; ++i )
00090 {
00091 for ( int j = 1 ; j < col ; ++j )
00092 {
00093 left( i, j-1) = right( i, 0);
00094 for ( int k = 0 ; k < col - j ; ++k )
00095 {
00096 right( i, k ) = ( right( i, k ) + right( i, k+1) ) * middle;
00097 }
00098 }
00099 left( i, col-1 ) = right( i, 0);
00100 }
00101 }
00102
00103 template<class matr_t>
00104 void deCasteljau( matr_t & left, matr_t & right, ByCol M )
00105 {
00106
00107 typedef typename matr_t::value_type coeff_t;
00108
00109 int col=right.nbcol(),row=right.nbrow();
00110 assert( row > 1 );
00111 coeff_t middle = coeff_t(0.5);
00112 for ( int i = 0 ; i < col ; ++i )
00113 {
00114 for ( int j = 1 ; j < row ; ++j )
00115 {
00116 left( j-1, i) = right( 0, i);
00117 for ( int k = 0 ; k < row - j ; ++k )
00118 {
00119 right( k, i ) = ( right( k, i ) + right( k+1, i ) ) * middle;
00120 }
00121 }
00122 left( row-1, i) = right( 0, i);
00123 }
00124 }
00125
00126
00127 template<class cell_t>
00128 void split(std::vector<cell_t >& L, const cell_t& c)
00129 {
00130 typedef typename cell_t::matr_t matr_t;
00131
00132
00133
00134
00135 matr_t M1(c.bz_pol),M2(c.bz_pol);
00136 deCasteljau(M1,M2,ByRow());
00137 matr_t M11(M1),M21(M1);
00138 deCasteljau(M11,M21,ByCol());
00139
00140
00141 L.push_back(cell_t(c.x0,(c.x0+c.x1)/2,c.y0,(c.y0+c.y1)/2,M11));
00142 L.push_back(cell_t(c.x0,(c.x0+c.x1)/2,(c.y0+c.y1)/2,c.y1,M21));
00143 matr_t M12(M2),M22(M2);
00144 deCasteljau(M12,M22,ByCol());
00145
00146
00147 L.push_back(cell_t((c.x0+c.x1)/2,c.x1,c.y0,(c.y0+c.y1)/2,M12));
00148 L.push_back(cell_t((c.x0+c.x1)/2,c.x1,(c.y0+c.y1)/2,c.y1,M22));
00149 }
00150 }
00151
00152 namespace bezier
00153 {
00154 template<class X>
00155 struct bzslv: public Isole<X>
00156 {
00157
00158 X eps;
00159 bzslv():eps(1e-4){}
00160 template<class UPOL>
00161 Seq<X> solve(const UPOL&p, X x0, X x1)
00162 {
00163 typedef typename UPOL::value_type coeff_t;
00164 bezier::rep1d<coeff_t> bz(p.size());
00165 VECTOR::assign(bz,p);
00166 Seq<X> sol;
00167 this->Run(bz,eps*(x1-x0),false);
00168 this->GetMidPoints(sol);
00169 RConvert( sol, x0,x1);
00170 return sol;
00171 }
00172 };
00173
00174 template<class X, class C=X, class USLV= bezier::bzslv<C> >
00175 struct sbd2d
00176 {
00177 typedef Seq<X> u_solution;
00178 typedef topology::point<X> point_t;
00179 typedef topology::point_graph<X> pointgraph_t;
00180 typedef C coeff_t;
00181 typedef cell2d<X,MatrDse<C> > cell_t;
00182
00183 USLV slv;
00184 X eps, asr;
00185 sbd2d():eps(1e-3), asr(1) {}
00186 sbd2d(X e):eps(e), asr(1) {}
00187 sbd2d(X e, X a):eps(e), asr(a) {}
00188
00189 cell_t pop(std::vector<cell_t>& L)
00190 {
00191 cell_t c=L[L.size()-1];
00192 L.resize(L.size()-1);
00193 return c;
00194 }
00195
00196 template<class mpol_t>
00197 void assign(cell_t& C,
00198 const mpol_t& p,
00199 const X& x0, const X& x1, const X& y0, const X&y1)
00200 {
00201 typedef typename mpol_t::coeff_t coeff_t;
00202 MatrDse<coeff_t> fq;
00203 monomial2bezier2d(fq,p,x0,x1,y0,y1);
00204 C.bz_pol =MatrDse<coeff_t>(fq.nbrow(),fq.nbcol());
00205 MATRIX::assign(C.bz_pol,fq);
00206 C.x0=x0; C.x1=x1;C.y0=y0;C.y1=y1;
00207 }
00208
00209 void connect_xy(pointgraph_t & g, const cell_t & c, bool b);
00210
00211 int c_sign(const cell_t&);
00212 int dx_sign(const cell_t&);
00213 int dy_sign(const cell_t&);
00214 int regularity(cell_t & c);
00215
00216 template<typename mpol_t>
00217 void run(pointgraph_t& g, const mpol_t & p,
00218 const X& x0, const X& x1, const X& y0, const X& y1);
00219 };
00220
00221
00222 template<class X, class U, class S>
00223 void sbd2d<X,U,S>::connect_xy(pointgraph_t & g, const cell_t & c, bool b)
00224 {
00225 typedef typename pointgraph_t::point_t point_t;
00226 point_t pt(2);
00227 int n0, n1;
00228 Seq<point_t> l;
00229
00230 Seq<X> sol0 =slv.solve((typename cell_t::row_type)row(c.bz_pol,0),c.x0,c.x1);
00231 pt[1]=c.y0;
00232 for(unsigned i=0;i<sol0.size();i++)
00233 {
00234 pt[0]=sol0[i];
00235 l.push_back(pt);
00236 }
00237 Seq<X> sol1 =slv.solve((typename cell_t::row_type)row(c.bz_pol,c.bz_pol.nbrow()-1),
00238 c.x0,c.x1);
00239 pt[1]=c.y1;
00240 for(unsigned i=0;i<sol1.size();i++)
00241 {
00242 pt[0]=sol1[i]; l.push_back(pt);
00243 }
00244 Seq<X> sol2 =slv.solve((typename cell_t::col_type)col(c.bz_pol,0),c.y0,c.y1);
00245 pt[0]=c.x0;
00246 for(unsigned i=0;i<sol2.size();i++)
00247 {
00248 pt[1]=sol2[i]; l.push_back(pt);
00249 }
00250 Seq<X> sol3 =slv.solve((typename cell_t::col_type)col(c.bz_pol,c.bz_pol.nbcol()-1),
00251 c.y0,c.y1);
00252 pt[0]=c.x1;
00253 for(unsigned i=0;i<sol3.size();i++)
00254 {
00255 pt[1]=sol3[i]; l.push_back(pt);
00256 }
00257 if(b)
00258 sort(l.begin(),l.end(), YX());
00259 else
00260 sort(l.begin(),l.end(), XY());
00261 if(l.size()%2 !=0)
00262 {
00263 std::cout<<c.x0<<" "<<c.x1<<" "<<c.y0<<" "<<c.y1<<"\n ";
00264 VECTOR::print(std::cout,l)<<" => ";
00265 Seq<point_t> ll;
00266 for(unsigned i=0;i<l.size()-1;i++)
00267 if(fabs(l[i][0]-l[i+1][0])> slv.eps ||
00268 fabs(l[i][1]-l[i+1][1])> slv.eps)
00269 ll.push_back(l[i]);
00270 ll.push_back(l[l.size()-1]);
00271 l=ll;
00272 VECTOR::print(std::cout,l)<<std::endl;
00273 }
00274 for(int i=0;i<l.size(); i+=2)
00275 {
00276 if(i<l.size()-1)
00277 {
00278 n0=g.insert(l[i]);
00279 n1=g.insert(l[i+1]);
00280 g.insert(n0,n1);
00281 }
00282 }
00283 }
00284
00285
00286 template<class X, class U, class S>
00287 template<class mpol_t>
00288 void sbd2d<X,U,S>::run(topology::point_graph<X>& r, const mpol_t & p,
00289 const X& x0, const X& x1, const X& y0, const X& y1)
00290 {
00291 typedef typename mpol_t::coeff_t CF_T;
00292 typedef cell2d<coeff_t> cell_t;
00293
00294
00295 cell_t C;
00296 assign(C,p,x0,x1,y0,y1);
00297
00298 double sz = size(C)*asr;
00299 std::vector<cell_t> L, L_sing;
00300 L.push_back(C);
00301 while(L.size())
00302 {
00303 C=pop(L);
00304 int s=regularity(C);
00305 if(s==0)
00306 {
00307 if(size(C)<eps)
00308 {
00309 L_sing.push_back(C);
00310 addbox(r,C);
00311 addcenter(r,C);
00312 }
00313 else
00314 BEZIER::split(L,C);
00315 }
00316 else
00317 {
00318 if(s%4==0)
00319 {
00320 if(size(C)> sz)
00321 split(L,C);
00322 else
00323 {
00324 connect_xy(r,C, (C.reg/4) %4);
00325 }
00326 }
00327 }
00328 }
00329
00330 std::cout <<"Sing cell: "<<L_sing.size()<<std::endl;
00331 }
00332
00333
00334 template<class X, class U, class S>
00335 int sbd2d<X,U,S>::c_sign(const cell_t & M)
00336 {
00337 typedef typename cell_t::coeff_t coeff_t;
00338 int s=0;
00339 bool run=true;
00340 for(unsigned i=0;i<M.nbrow() && run;i++)
00341 for(unsigned j=0;j<M.nbcol()&& run;j++)
00342 if(M(i,j)>(0))
00343 {
00344 s= 1;
00345 run=false;
00346 }
00347 else if(M(i,j)<(0))
00348 {
00349 s= -1;
00350 run=false;
00351 }
00352
00353
00354 for(unsigned i=0;i<M.nbrow();i++)
00355 for(unsigned j=0;j<M.nbcol();j++)
00356 if(sign(M(i,j))==-s)
00357 return 0;
00358 return s;
00359 }
00360
00361 template<class X, class U, class S>
00362 int sbd2d<X,U,S>::dx_sign(const cell_t & M)
00363 {
00364 typedef typename cell_t::coeff_t coeff_t;
00365 int s=0; bool run=true;
00366 for(int i=0;i<M.nbrow() && run;i++)
00367 for(int j=0;j<M.nbcol()-1 && run;j++)
00368 if(M(i,j+1)>M(i,j))
00369 {
00370 s= 1;
00371 run=false;
00372 }
00373 else if(M(i,j+1)<M(i,j))
00374 {
00375 s= -1;
00376 run=false;
00377 }
00378 for(int i=0;i<M.nbrow();i++)
00379 for(int j=0;j<M.nbcol()-1;j++)
00380 {
00381 if(sign(M(i,j+1)-M(i,j))==-s)
00382 return 0;
00383 }
00384 return s;
00385 }
00386
00387 template<class X, class U, class S>
00388 int sbd2d<X,U,S>::dy_sign(const cell_t & M)
00389 {
00390 typedef typename cell_t::coeff_t coeff_t;
00391 int s=0; bool run=true;
00392 for(unsigned j=0;j<M.nbcol()&& run;j++)
00393 for(int i=0;i<M.nbrow()-1&& run;i++)
00394 if(M(i+1,j)>M(i,j))
00395 {
00396 s= 1;
00397 run=false;
00398 }
00399 else if(M(i+1,j)<M(i,j))
00400 {
00401 s= -1;
00402 run=false;
00403 }
00404
00405 for(unsigned j=0;j<M.nbcol();j++)
00406 for(int i=0;i<M.nbrow()-1;i++)
00407 {
00408 if(sign(M(i+1,j)-M(i,j))!= s)
00409 return 0;
00410 }
00411 return s;
00412 }
00413
00414 template<class X,class U,class S>
00415 int sbd2d<X,U,S>::regularity(cell_t& c)
00416 {
00417 int s0= c_sign(c);
00418 s0 += 4*dx_sign(c);
00419 s0 += 16*dy_sign(c);
00420 c.reg=s0;
00421 return s0;
00422 }
00423
00424 }
00425
00426 template<class cell_t, class X>
00427 void addbox(topology::point_graph<X>& g, const cell_t & c)
00428 {
00429 typedef typename topology::point_graph<X> pointgraph_t;
00430 typedef typename pointgraph_t::point_t point_t;
00431 point_t p1(2), p2(2);
00432 unsigned n1, n2;
00433
00434 X a0=c.x0,a1=c.x1,b0=c.y0,b1=c.y1;
00435
00436 p1[0]=a0;p1[1]=b0;
00437 p2[0]=a0;p2[1]=b1;
00438 n1=g.insert(p1);
00439 n2=g.insert(p2);
00440 g.insert(n1,n2);
00441
00442 p2[0]=a1;p2[1]=b0;
00443 n2=g.insert(p2);
00444 g.insert(n1,n2);
00445
00446 p1[0]=a0;p1[1]=b1;
00447 p2[0]=a1;p2[1]=b1;
00448 n1=g.insert(p1);
00449 n2=g.insert(p2);
00450 g.insert(n1,n2);
00451
00452 p1[0]=a1;p2[1]=b0;
00453 n1=g.insert(p2);
00454 g.insert(n1,n2);
00455 }
00456
00457 template<class cell_t, class X>
00458 void addplus(topology::point_graph<X>& g, const cell_t & c)
00459 {
00460 topology::point<X> p1(2),p2(2);
00461 p1[0]=(c.x0+c.x1)/2;
00462 p1[1]=(c.y0+3*c.y1)/4;
00463 int n1 = g.insert(p1);
00464 p2[0]=(c.x0+c.x1)/2;
00465 p2[1]=(3*c.y0+c.y1)/4;
00466 int n2= g.insert(p2);
00467 g.insert(n1,n2);
00468
00469 p1[0]=(3*c.x0+c.x1)/4;
00470 p1[1]=(c.y0+c.y1)/2;
00471 n1 = g.insert(p1);
00472 p2[0]=(c.x0+3*c.x1)/4;
00473 p2[1]=(c.y0+c.y1)/2;
00474 n2= g.insert(p2);
00475 g.insert(n1,n2);
00476
00477 }
00478
00479 template<class cell_t, class X>
00480 void addmoins(topology::point_graph<X>& g, const cell_t & c)
00481 {
00482 topology::point<X> p1(2),p2(2);
00483 p1[0]=(3*c.x0+c.x1)/4;
00484 p1[1]=(c.y0+c.y1)/2;
00485 int n1 = g.insert(p1);
00486 p2[0]=(c.x0+3*c.x1)/4;
00487 p2[1]=(c.y0+c.y1)/2;
00488 int n2= g.insert(p2);
00489 g.insert(n1,n2);
00490
00491 }
00492
00493 template<class cell_t, class X>
00494 void addcenter(topology::point_graph<X>& g, const cell_t & c)
00495 {
00496 topology::point<X> pt(2);
00497 pt[0]=(c.x0+c.x1)/2;
00498 pt[1]=(c.y0+c.y1)/2;
00499 g.insert(pt);
00500 }
00501
00502 __END_NAMESPACE_SYNAPS
00503
00504 #endif //synaps_topology_sbd2d_H