00001
00002
00003
00004
00005 #ifndef synaps_topology_TopSbdBdg2d_h_
00006 #define synaps_topology_TopSbdBdg2d_h_
00007
00008 #include <synaps/init.h>
00009 #include <synaps/arithm/gmp.h>
00010 #include <synaps/mpol.h>
00011 #include <synaps/topology/sbd2d.h>
00012 #include <synaps/arithm/Interval.h>
00013
00014 __BEGIN_NAMESPACE_SYNAPS
00015
00016 template<class O, class X, class I> struct TopSbdBdg2d;
00017
00018 template<class X, int N>
00019 void print_mx(const MatrDse<Interval<X,N> >& M)
00020 {
00021
00022 double dmx=0, imx=0,imn=1;
00023 for(unsigned i=0; i< M.nbrow();i++)
00024 for(unsigned j=0; j< M.nbcol();j++)
00025 {
00026 imx= std::max(imx,M(i,j).sup()-M(i,j).inf());
00027 imn= std::min(imn,M(i,j).sup()-M(i,j).inf());
00028 if( M(i,j).sup()<0)
00029 dmx = std::max(dmx,-M(i,j).sup());
00030 else
00031 dmx = std::max(dmx,M(i,j).sup());
00032 }
00033 std::cout<<"DMax "<<dmx<<std::endl;
00034 std::cout<<"IMax "<<imx<<std::endl;
00035 std::cout<<"IMin "<<imn<<std::endl;
00036 }
00037
00038 template<class X, int N>
00039 void print_sign(const MatrDse<Interval<X,N> >& M)
00040 {
00041 for(unsigned i=0;i<M.nbrow();i++)
00042 {
00043 std::cout<<"[";
00044 for(unsigned j=0; j<M.nbcol();j++)
00045 {
00046 if(M(i,j).sup()<0)
00047 std::cout<<"-";
00048 else if(M(i,j).inf()>=0)
00049 std::cout<<"+";
00050 else if(M(i,j).inf()<0 && M(i,j).sup()>0)
00051 std::cout<<"?";
00052 else
00053 std::cout<<" ";
00054
00055 }
00056 std::cout<<"]"<<std::endl;
00057 }
00058 }
00059 template<class Int>
00060 bool IsSmall (const Int& I)
00061 {
00062 return (I.inf()<= 0 && I.sup() >= 0 && (I.inf()!=I.sup()));
00063 }
00064
00065 template<class X, int N>
00066 int sign(const Interval<X,N>& I)
00067 {
00068 if(I.sup()< (0))
00069 return -1;
00070 else if(I.inf()>(0))
00071 return 1;
00072 else if(I.sup()==0 && I.inf()==0)
00073 return 0;
00074 else
00075 return 2;
00076 }
00077
00078 namespace bezier
00079 {
00080 template<class X>
00081 struct bzislv: public Isole<X>
00082 {
00083 X eps;
00084 bzislv():eps(1e-4){}
00085
00086 template<class UPOL>
00087 Seq<X> solve(const UPOL& p, X x0, X x1)
00088 {
00089 typedef typename UPOL::value_type coeff_t;
00090 bezier::rep1d<X> ubp(p.size()), dbp(p.size());
00091 for(unsigned i=0;i<p.size();i++)
00092 {
00093 ubp[i]=p[i].sup();
00094 dbp[i]=p[i].inf();
00095 assert(dbp[i]<= ubp[i]);
00096 }
00097 this->Run(dbp,ubp,eps,false);
00098 Seq<X> sol; sol.resize(0);
00099 this->GetMidPoints(sol);
00100 for (unsigned i = 0; i < sol.size(); i ++ )
00101 {
00102 if ( this->m_res[i].t )
00103 std::cout<<" xxx"<<std::endl;
00104 }
00105 RConvert(sol,x0,x1);
00106 return sol;
00107 }
00108 };
00109
00110 }
00111
00112 template<class O, class X=double, class I= Interval<X,3> >
00113 struct TopSbdBdg2d : public bezier::sbd2d<X,I,bezier::bzislv<X> >
00114 {
00115 typedef topology::point<X> point_t;
00116 typedef topology::point_graph<X> pointgraph_t;
00117 typedef I interval_t;
00118 typedef interval_t coeff_t;
00119 typedef cell2d<X,MatrDse<coeff_t> > cell_t;
00120 bool draw_box;
00121
00122 TopSbdBdg2d<O,X,I>(bool bx=true):bezier::sbd2d<X,I,bezier::bzislv<X> >(1.e-3,0.1), draw_box(bx) {}
00123 TopSbdBdg2d<O,X,I>(const X& e, const X& a, bool bx=true):bezier::sbd2d<X,I,bezier::bzislv<X> >(e,a), draw_box(bx) {}
00124
00125
00126 int c_sign(const cell_t&);
00127 int dx_sign(const cell_t&);
00128 int dy_sign(const cell_t&);
00129 int regularity(cell_t & c)
00130 {
00131 int s0=c_sign(c), s1=4*dx_sign(c);
00132 s1+=16*dy_sign(c);
00133 if(s1)
00134 {
00135 if(s0==1 || s0 == -1)
00136 s0+=s1;
00137 else
00138 s0=s1;
00139 }
00140 c.reg=s0;
00141
00142 return s0;
00143 }
00144
00145 template<class mpol_t>
00146 void assign(cell_t& C,
00147 const mpol_t& p,
00148 const X& x0, const X& x1, const X& y0, const X&y1)
00149 {
00150 typedef typename mpol_t::coeff_t Coeff_t;
00151 MatrDse<Coeff_t> fq;
00152 BEZIER::monomial2bezier2d(fq,p,Coeff_t(x0),Coeff_t(x1),Coeff_t(y0),Coeff_t(y1));
00153 Coeff_t mx=0;
00154 for(unsigned i=0; i< fq.nbrow();i++)
00155 for(unsigned j=0; j< fq.nbcol();j++)
00156 if( fq(i,j)<Coeff_t(0))
00157 mx = std::max(mx,Coeff_t(-fq(i,j)));
00158 else
00159 mx = std::max(mx,fq(i,j));
00160 std::cout<<"Max "<<mx<<std::endl;
00161 fq/=mx;
00162 mx=Coeff_t(0);
00163 C.bz_pol =MatrDse<coeff_t>(fq.nbrow(),fq.nbcol());
00164 MATRIX::assign(C.bz_pol,fq);
00165 C.x0=x0; C.x1=x1;C.y0=y0;C.y1=y1;
00166
00167 }
00168
00169 template<typename mpol_t>
00170 void run(topology::point_graph<double>& r, const mpol_t & p,
00171 const X& x0, const X& x1, const X& y0, const X& y1);
00172 };
00173
00174
00175 template<class O, class X, class I>
00176 template<class mpol_t>
00177 void TopSbdBdg2d<O,X,I>::run(topology::point_graph<double>& r, const mpol_t & p,
00178 const X& x0, const X& x1, const X& y0, const X& y1)
00179 {
00180 cell_t C;
00181 assign(C,p,x0,x1,y0,y1);
00182
00183 double min_sz=size(C)*this->eps,
00184 max_sz = size(C)*this->asr;
00185 std::cout<<"Min size "<<min_sz<<std::endl;
00186 std::cout<<"Max size "<<max_sz<<std::endl;
00187 if(draw_box) std::cout<<"Drawing singular boxes"<<std::endl;
00188 std::vector<cell_t> L, L_sing;
00189 L.push_back(C);
00190 while(L.size())
00191 {
00192 C=this->pop(L);
00193 int s=regularity(C);
00194
00195 if(s==0)
00196 {
00197
00198 if(size(C)<min_sz)
00199 {
00200 L_sing.push_back(C);
00201 if(draw_box)
00202 {
00203 addbox(r,C);
00204 addcenter(r,C);
00205 }
00206 }
00207 else
00208 BEZIER::split(L,C);
00209 }
00210 else
00211 {
00212 if(s==2)
00213 {
00214 std::cout<<"Zoom [ "<<C.x0<<" "<<C.x1<<" "<<C.y0<<" "<<C.y1<<" ]"<<std::endl;
00215 assign(C,p,C.x0,C.x1,C.y0,C.y1);
00216 L.push_back(C);
00217 if(draw_box)
00218 {
00219 addbox(r,C);
00220
00221 }
00222 }
00223 else if(s%4==0)
00224 {
00225 if(size(C)> max_sz)
00226 BEZIER::split(L,C);
00227 else
00228 {
00229 connect_xy(r,C,(C.reg/4)%4);
00230
00231 }
00232 }
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242 }
00243 }
00244 std::cout <<"Sing cell: "<<L_sing.size()<<std::endl;
00245 point_t pt(2);
00246 }
00247
00248 template<class O, class X, class I>
00249 int TopSbdBdg2d<O,X,I>::c_sign(const cell_t& M)
00250 {
00251 typedef typename cell_t::coeff_t coeff_t;
00252 int s=0;
00253 bool loop=true;
00254 unsigned i,j;
00255
00256 for(i=0;i<M.nbrow() && IsSmall(M(i,0)); i++);
00257 if(i==M.nbrow()) return 2;
00258 for(i=0;i<M.nbrow() && IsSmall(M(i,M.nbcol()-1)); i++);
00259 if(i==M.nbrow()) return 2;
00260
00261 for(j=0;j<M.nbcol() && IsSmall(M(0,j)); j++);
00262 if(j==M.nbcol()) return 2;
00263 for(j=0;j<M.nbcol() && IsSmall(M(M.nbrow()-1,j)); j++);
00264 if(j==M.nbcol()) return 2;
00265
00266 for(i=0;i<M.nbrow() && loop;i++)
00267 for(j=0;j<M.nbcol()&& loop;j++)
00268 if(M(i,j).inf()>(0))
00269 {
00270 s= 1;
00271 loop=false;
00272 }
00273 else if(M(i,j).sup()<(0))
00274 {
00275 s= -1;
00276 loop=false;
00277 }
00278 for(i=0;i<M.nbrow();i++)
00279 for(j=0;j<M.nbcol();j++)
00280 if(s==1 && sign(M(i,j))==-1)
00281 return 0;
00282 else if(s==-1 && sign(M(i,j))==1)
00283 return 0;
00284 return s;
00285 }
00286
00287
00288 template<class O, class X,class I>
00289 int TopSbdBdg2d<O,X,I>::dx_sign(const cell_t & M)
00290 {
00291 typedef typename cell_t::coeff_t coeff_t;
00292 int s=0; bool loop=true;
00293 for(int i=0;i<M.nbrow() && loop;i++)
00294 for(int j=0;j<M.nbcol()-1 && loop;j++)
00295 if(M(i,j+1).inf()>M(i,j).sup())
00296 {
00297 s= 1;
00298 loop=false;
00299 }
00300 else if(M(i,j+1).sup()<M(i,j).inf())
00301 {
00302 s= -1;
00303 loop=false;
00304 }
00305 else if(IsSmall(interval_t(M(i,j+1)-M(i,j))))
00306 {
00307 s=0;
00308 loop=false;
00309 }
00310 for(int i=0;i<M.nbrow();i++)
00311 for(int j=0;j<M.nbcol()-1;j++)
00312 {
00313 interval_t D=M(i,j+1)-M(i,j);
00314 if(sign(D)!=s)
00315 return 0;
00316 }
00317 return s;
00318 }
00319
00320 template<class O,class X, class I>
00321 int TopSbdBdg2d<O,X,I>::dy_sign(const cell_t & M)
00322 {
00323 typedef typename cell_t::coeff_t coeff_t;
00324 int s=0; bool loop=true;
00325 for(unsigned j=0;j<M.nbcol()&& loop;j++)
00326 for(int i=0;i<M.nbrow()-1&& loop;i++)
00327 if(M(i+1,j).inf()>M(i,j).sup())
00328 {
00329 s= 1;
00330 loop=false;
00331 }
00332 else if(M(i+1,j).sup()<M(i,j).sup())
00333 {
00334 s= -1;
00335 loop=false;
00336 }
00337 else if(IsSmall(interval_t(M(i+1,j)-M(i,j))))
00338 {
00339 s=0;
00340 loop=false;
00341 }
00342 for(unsigned j=0;j<M.nbcol();j++)
00343 for(int i=0;i<M.nbrow()-1;i++)
00344 {
00345 interval_t D=M(i+1,j)-M(i,j);
00346 if(sign(D)!= s)
00347 return 0;
00348 }
00349 return s;
00350 }
00351
00352
00353 namespace topology
00354 {
00355 template<class C, class U>
00356 void assign(point_graph<C>& r, const MPol<U>& s,
00357 TopSbdBdg2d<C> mth,
00358 C x0, C x1, C y0, C y1,
00359 unsigned Z=1)
00360 {
00361 typedef QQ coeff_t;
00362 typedef MPol<coeff_t> mpol_t;
00363 mpol_t p; let::assign(p,s);
00364 double eps=mth.eps,a0=x0,b0=x1,a1=y0,b1=y1;
00365 C d0=(b0-a0)/Z, d1=(b1-a1)/Z;
00366 for(unsigned i=0;i<Z;i++)
00367 for(unsigned j=0;j<Z;j++)
00368 {
00369 mth.run(r,p,a0+d0*i,a0+(i+1)*d0,a1+d1*j,a1+d1*(j+1));
00370 }
00371 }
00372 }
00373 __END_NAMESPACE_SYNAPS
00374
00375 #endif //synaps_topology_TopSbdBdg2d_h_