00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 #ifndef shape_implicit_solver_hpp
00013 #define shape_implicit_solver_hpp
00014
00015 # include <realroot/Interval.hpp>
00016 # include <realroot/ring_bernstein_tensor.hpp>
00017 # include <realroot/GMP.hpp>
00018 # include <realroot/solver_uv_sleeve.hpp>
00019 # include <realroot/solver_mv_bernstein.hpp>
00020 # include <shape/point.hpp>
00021 # include <shape/bounding_box.hpp>
00022 # include <shape/box_face.hpp>
00023 # include <shape/topological_degree.hpp>
00024
00025 # define TMPL template<class C, class V>
00026 # define SELF solver_implicit<C,V>
00027
00028
00029 namespace mmx { namespace shape {
00030
00031
00032 inline double lower(double x){return x;}
00033 inline double upper(double x){return x;}
00034
00035 template<class S1, class S2>
00036 point<double>* point_on_edge(const bounding_box<double>& bx, double u, const S1& s1, const S2& s2)
00037 {
00038 typedef point<double> Point;
00039 double p[3];
00040 p[S1::var] = S1::eval(bx);
00041 p[S2::var] = S2::eval(bx);
00042 int v= 3-S1::var-S2::var;
00043 p[v]= lower(bx,v)*(1-u) + upper(bx,v)*u;
00044
00045 return new Point(p[0],p[1],p[2]);
00046
00047 }
00048
00049 TMPL struct solver_implicit {
00050
00051 typedef C Scalar;
00052 typedef REF_OF(V) Ref;
00053 typedef vertex<C,3,V> Point;
00054 typedef bounding_box<C,Ref> BoundingBox;
00055 typedef box_face<C,Ref> BoxFace;
00056
00057 static BoxFace north_edge;
00058 static BoxFace south_edge;
00059 static BoxFace west_edge;
00060 static BoxFace east_edge;
00061 static BoxFace front_edge;
00062 static BoxFace back_edge;
00063
00064 static BoxFace north_back_edge;
00065 static BoxFace south_back_edge;
00066 static BoxFace west_back_edge;
00067 static BoxFace east_back_edge;
00068
00069 static BoxFace north_front_edge;
00070 static BoxFace south_front_edge;
00071 static BoxFace west_front_edge;
00072 static BoxFace east_front_edge;
00073
00074 static BoxFace north_west_edge;
00075 static BoxFace south_west_edge;
00076 static BoxFace north_east_edge;
00077 static BoxFace south_east_edge;
00078
00079 static BoxFace north_face;
00080 static BoxFace south_face;
00081 static BoxFace west_face;
00082 static BoxFace east_face;
00083 static BoxFace front_face;
00084 static BoxFace back_face;
00085
00086 static inline Point* new_point(double a, double b, double c) {
00087 return new Point(a,b,c);
00088 }
00089
00090
00091
00092
00093 template<class MultivariateDenseBernstein> static
00094 void edge_point(Seq<Point *>& lp,
00095 const MultivariateDenseBernstein & p,
00096 const BoxFace & E,
00097 const BoundingBox & bx) {
00098
00099 MultivariateDenseBernstein f;
00100
00101 tensor::face(f,p,E.cvar(0),E.side(0));
00102
00103 polynomial<double, with<Bernstein> > dw(1,f.size()-1), up(1,f.size()-1);
00104 for(unsigned i=0; i<f.size();i++) {
00105 up[i]= upper(f[i]);
00106 dw[i]= lower(f[i]);
00107 }
00108
00109 double
00110 U_b = E.lower(bx),
00111 V_b = E.upper(bx);
00112
00113 Seq<double> sol;
00114 solver< double, Sleeve<Approximate> >::solve(sol,dw,up,U_b,V_b);
00115
00116
00117
00118
00119
00120
00121
00122 Approximate approx;
00123 std::vector<double> pt(3);
00124 for(unsigned i=0;i<sol.size(); i++)
00125 if (i==0 || sol[i]-sol[i-1]> approx.eps){
00126 if(true) {
00127 lp<< E.new_point(bx,sol[i]);
00128 } else {
00129 std::cout<<"Not valid: ["<<(*E.new_point(bx,sol[i])) <<"] "
00130 <<bx<<" "
00131 <<E.m_v[0]<< " "<< E.m_s[0] <<" "
00132 <<E.m_v[1]<< " "<< E.m_s[1] << " "
00133 <<std::endl;
00134 lp<< E.new_point(bx,sol[i]);
00135 }
00136 }
00137 }
00138
00139 template<class MultivariateDenseBernstein> static
00140 int edge_sign_var(const MultivariateDenseBernstein & p,
00141 const BoxFace & E) {
00142 MultivariateDenseBernstein f;
00143 tensor::face(f,p,E.cvar(0),E.side(0));
00144 return sign_variation(f.begin(), f.end());
00145
00146 }
00147
00148 template<class MultivariateDenseBernstein> static
00149 void face_point(Seq<Point *>& sol,
00150 const Seq<MultivariateDenseBernstein>& eqs,
00151 const BoxFace& F,
00152 const BoundingBox& bx) {
00153
00154 typedef polynomial< double, with<Bernstein> > Polynomial;
00155
00156 Polynomial f, p;
00157 Seq<Polynomial> sys;
00158 for(unsigned i=0;i< eqs.size(); i++) {
00159 let::assign(p.rep(), eqs[i].rep());
00160 tensor::face(f.rep(), p.rep(), F.cvar(0), F.side(0));
00161 sys<<f;
00162 }
00163
00164 solver<double, ProjRd<SBD_RDRDL> > slv;
00165 std::vector<double> res;
00166 slv.solve_bernstein(res,sys);
00167 Approximate approx;
00168
00169 for(int i=0;i<(int)res.size()-1;i+=2) {
00170
00171
00172
00173
00174
00175
00176
00177 Point * p = F.new_point(bx, res[i], res[i+1]);
00178 if (F.is_valid(*p, bx)) sol<< p;
00179
00180 }
00181 }
00182
00183 template<class Polynomial>
00184 static void intersection(Seq<Point* >& sol,
00185 const Polynomial& f1,
00186 const Polynomial& f2, const BoundingBox & bx) {
00187
00188 typedef polynomial< double, with<Bernstein> > MultivariateDenseBernstein;
00189 MultivariateDenseBernstein p1, p2;
00190
00191
00192 let::assign(p1.rep(), f1.rep() );
00193 let::assign(p2.rep(), f2.rep() );
00194
00195 Seq<MultivariateDenseBernstein> sys;
00196 sys<<p1<<p2;
00197
00198 solver<double, ProjRd<SBD_RDRDL> > slv;
00199 std::vector<double> res;
00200 slv.solve_bernstein(res,sys);
00201
00202 double x_sc= (bx.xmax()-bx.xmin()), y_sc= (bx.ymax()-bx.ymin());
00203 for(unsigned i=0;i<res.size();i+=2)
00204 {
00205 sol << new_point(bx.xmin()+res[i]*x_sc, bx.ymin()+res[i+1]*y_sc, 0.0);
00206 }
00207 }
00208
00209 template<class Polynomial>
00210 static void extremal(Seq<Point* >& sol,
00211 const Polynomial& pol, const BoundingBox& b) {
00212
00213 typedef polynomial< double, with<Bernstein> > MultivariateDenseBernstein;
00214
00215
00216
00217
00218 typedef polynomial< GMP::rational, with<Bernstein> > MultivariateDenseBernsteinQ;
00219
00220 Seq<mmx::GMP::rational> bx;
00221 bx<<as<mmx::GMP::rational>(b.xmin());
00222 bx<<as<mmx::GMP::rational>(b.xmax());
00223 bx<<as<mmx::GMP::rational>(b.ymin());
00224 bx<<as<mmx::GMP::rational>(b.ymax());
00225
00226
00227
00228
00229
00230
00231 Polynomial dxpol = diff(pol,0), dypol = diff(pol,1);
00232 tensor::bernstein<mmx::GMP::rational>
00233 dxP(dxpol.rep(),bx),dyP(dypol.rep(),bx);
00234
00235 MultivariateDenseBernstein dxp,dyp;
00236 let::assign(dxp.rep(), dxP);
00237 let::assign(dyp.rep(), dyP);
00238
00239 Seq<MultivariateDenseBernstein> sys;
00240
00241 sys<<dxp<<dyp;
00242
00243 solver<double, ProjRd<SBD_RDRDL> > slv;
00244
00245 std::vector<double> res;
00246 slv.solve_bernstein(res,sys);
00247
00248 double x_sc= (b.xmax()-b.xmin()), y_sc= (b.ymax()-b.ymin());
00249 for(unsigned i=0;i<res.size();i+=2) {
00250 sol << new_point(b.xmin()+res[i]*x_sc, b.ymin()+res[i+1]*y_sc, 0.0);
00251 }
00252 }
00253
00254 template<class Polynomial>
00255 static void singular(Seq<Point* >& sol,
00256 const Polynomial& pol, const BoundingBox& b) {
00257
00258 typedef polynomial< double, with<Bernstein> > MultivariateDenseBernstein;
00259 typedef polynomial< GMP::rational, with<Bernstein> > MultivariateDenseBernsteinQ;
00260
00261 Seq<mmx::GMP::rational> bx;
00262 bx<<as<mmx::GMP::rational>(b.xmin());
00263 bx<<as<mmx::GMP::rational>(b.xmax());
00264 bx<<as<mmx::GMP::rational>(b.ymin());
00265 bx<<as<mmx::GMP::rational>(b.ymax());
00266
00267 Polynomial dxpol = diff(pol,0), dypol = diff(pol,1);
00268 tensor::bernstein<mmx::GMP::rational>
00269 P(pol.rep(),bx), dxP(dxpol.rep(),bx),dyP(dypol.rep(),bx);
00270
00271 MultivariateDenseBernstein p,dxp,dyp;
00272 let::assign(p.rep() , P);
00273 let::assign(dxp.rep(), dxP);
00274 let::assign(dyp.rep(), dyP);
00275
00276 Seq<MultivariateDenseBernstein> sys;
00277 sys<<p<<dxp<<dyp;
00278
00279 solver<double, ProjRd<SBD_RDRDL> > slv;
00280
00281
00282 std::vector<double> res;
00283 slv.solve_bernstein(res,sys);
00284
00285 double x_sc= (b.xmax()-b.xmin()), y_sc= (b.ymax()-b.ymin());
00286 for(unsigned i=0;i<res.size();i+=2) {
00287 sol << new_point(b.xmin()+res[i]*x_sc, b.ymin()+res[i+1]*y_sc, 0.0);
00288 }
00289 }
00290
00291
00292 template<class MultivariateDenseBernstein> static
00293 void common_edge_point(Seq<Point *>& lp,
00294 const MultivariateDenseBernstein & p,
00295 const BoxFace & E,
00296 const BoundingBox & bx1,
00297 const BoundingBox & bx2 ) {
00298 MultivariateDenseBernstein f;
00299
00300 tensor::face(f,p,E.cvar(0),E.side(0));
00301
00302 polynomial<double, with<Bernstein> > dw(1,f.size()-1), up(1,f.size()-1);
00303 for(unsigned i=0; i<f.size();i++) {
00304 up[i]= upper(f[i]);
00305 dw[i]= lower(f[i]);
00306
00307 }
00308
00309 double
00310 U_b = std::max( E.lower(bx1), E.lower(bx2) ) ,
00311 V_b = std::min( E.upper(bx1), E.upper(bx2) ) ;
00312
00313 Seq<double> sol;
00314
00315 solver< double, Sleeve<Approximate> >::solve(sol,dw,up,U_b,V_b);
00316
00317 sol.sort();
00318
00319 Approximate approx;
00320 for(unsigned i=0;i<sol.size(); i++)
00321 {
00322 lp<< E.new_point(bx1,sol[i]);
00323 }
00324 }
00325 };
00326
00327 # define REF REF_OF(V)
00328 TMPL box_face<C,REF> SELF::north_face(1,1);
00329 TMPL box_face<C,REF> SELF::south_face(1,0);
00330 TMPL box_face<C,REF> SELF::west_face (0,0);
00331 TMPL box_face<C,REF> SELF::east_face (0,1);
00332 TMPL box_face<C,REF> SELF::front_face(2,1);
00333 TMPL box_face<C,REF> SELF::back_face(2,0);
00334
00335
00336 TMPL box_face<C,REF> SELF::north_edge(1,1,2,0);
00337 TMPL box_face<C,REF> SELF::south_edge(1,0,2,0);
00338 TMPL box_face<C,REF> SELF::west_edge(0,0,2,0);
00339 TMPL box_face<C,REF> SELF::east_edge(0,1,2,0);
00340
00341 TMPL box_face<C,REF> SELF::front_edge(2,1,0,0);
00342 TMPL box_face<C,REF> SELF::back_edge (2,0,0,0);
00343
00344
00345 TMPL box_face<C,REF> SELF::north_back_edge (1,1,2,0);
00346 TMPL box_face<C,REF> SELF::south_back_edge (1,0,2,0);
00347 TMPL box_face<C,REF> SELF::west_back_edge (0,0,2,0);
00348 TMPL box_face<C,REF> SELF::east_back_edge (0,1,2,0);
00349
00350 TMPL box_face<C,REF> SELF::north_front_edge(1,1,2,1);
00351 TMPL box_face<C,REF> SELF::south_front_edge(1,0,2,1);
00352 TMPL box_face<C,REF> SELF::west_front_edge (0,0,2,1);
00353 TMPL box_face<C,REF> SELF::east_front_edge (0,1,2,1);
00354
00355 TMPL box_face<C,REF> SELF::north_west_edge (0,0,1,1);
00356 TMPL box_face<C,REF> SELF::south_west_edge (0,0,1,0);
00357 TMPL box_face<C,REF> SELF::north_east_edge (0,1,1,1);
00358 TMPL box_face<C,REF> SELF::south_east_edge (0,1,1,0);
00359
00360
00361 }
00362 }
00363
00364 # undef TMPL
00365 # undef SELF
00366 # undef REF
00367 # undef BoxFace
00368 # endif //shape_implicit_solver
