solver_fatarcs.hpp

Go to the documentation of this file.

# include <realroot/fatarcs.hpp>
# include <realroot/polynomial_bernstein.hpp>
# include <realroot/polynomial_tensor.hpp>

#pragma once

using namespace mmx;


// solver_fatarcs struct

template<class C>
struct solver_mv_fatarcs
{
 
  typedef std::vector<C>                   vec_t;
  typedef domain<C>                        box_t;
  typedef polynomial< C ,with<Bernstein> > bernstein_t;
  typedef box_rep< C >                     box_rep_t;
  typedef arc_rep< C >                     arc_rep_t;
  typedef std::stack< box_t >              stbox_t;
  typedef Seq< box_t >                     seqbox_t;

  C epsilon;
  bernstein_t  poly1,  poly2;
  
  //CONSTRUCTOR
  solver_mv_fatarcs(bernstein_t p1, bernstein_t p2, C eps){poly1=p1; poly2=p2; epsilon=eps; };


  //FCTS

  box_t box_gen(box_rep_t m1, box_rep_t m2, int MTH)
  {

    arc_rep_t mc1,mc2;
    Seq<vec_t> l1=m1.event_list();
    if(l1.size()==4){mc1=median(m1, l1, MTH);}else{mc1=arc_rep_t();};
    Seq<vec_t> l2=m2.event_list();
    if(l2.size()==4){mc2=median(m2, l2, MTH);}else{mc2=arc_rep_t();};
   
    box_t  bbox=m1.box;
  
    if(is_arc(mc1) && is_arc(mc2)){ 

     C c1=m1.min_grad();//std::cout << c1 <<std::endl;
     C c2=m2.min_grad();//std::cout << c2 <<std::endl;
     if(c1>C(0) && c2>C(0) ){
    C r1=m1.max_eval(mc1)/sqrt(c1);//std::cout << r1 <<std::endl;
        C r2=m2.max_eval(mc2)/sqrt(c2);//std::cout << r2 <<std::endl;
        Seq<vec_t> extrema;
        extrema=extpts(mc1,r1,mc2,r2);

        bbox=minmax_box(extrema, m1.box);
      }else{/*std::cout << "no grad "<<std::endl;*/} 
    }else{/*std::cout << "no med "<<std::endl;*/}

    return(bbox);
 };/*box around fatarc intersection*/



void prepro(bernstein_t & pol1, bernstein_t & pol2, vec_t midpt)
  { 
    bernstein_t p1=pol1,p2=pol2;
    C c1x,c1y,c2x,c2y;
  tensor::eval(c1x, diff(p1,0).rep(),midpt);
  tensor::eval(c2x, diff(p2,0).rep(),midpt);
  tensor::eval(c1y, diff(p1,1).rep(),midpt);
  tensor::eval(c2y, diff(p2,1).rep(),midpt);
  if(c1x!=c1y && c2x!=c2y){
  pol1=c1x*p1+c2x*p2;
  pol2=c1y*p1+c2y*p2; }
};


  seqbox_t solver(int MTH){
  stbox_t list, mo;
  seqbox_t moseq;
  
   int subdiv=0;
   int arcgen=0;
   int it=0;

   box_t unit(int(2));
   box_t p[4]; 
    list.push(unit);
 
    while(!list.empty()){it++;
     
      box_t rec=list.top();  
    
      list.pop();

      box_rep_t m1(poly1, rec); 
      box_rep_t m2(poly2, rec);
     
      if(m1.not_empty() && m2.not_empty()){

      // prepro( m1.poly, m2.poly, (rec.llc()+rec.urc())/((coeff_t)(2))); 
      // std::cout <<"b"<<std::endl;
      box_t newrec=box_gen(m1, m2, MTH); arcgen++; 

       
      if(newrec.dim()==0){ //std::cout <<"no fatarc intersection"<<std::endl;
      }else{   //std::cout <<"new"<<std::endl;newrec.print(5);
          if(newrec.diam() <( rec.diam()/ C(2) )  ){;

        if(newrec.diam()<epsilon){mo.push(newrec);//std::cout <<"small1"<<std::endl; newrec.print(6);
        }
        else{list.push(newrec);//std::cout <<"list new"<<std::endl;
        }

      }
      else{
        if(rec.diam()<epsilon){mo.push(rec);}
        else{rec.split(p);subdiv++;/*std::cout <<"sub1"<<std::endl;*/
           for(int i=0; i<4; i++){list.push(p[3-i]);}
        }
        
      };

      };
         
      }else{//std::cout <<"empty"<<std::endl;
      };
    }
   
 
       
std::cout << " * * * * * "<<std::endl;
std::cout << "iteration: "<<it<<std::endl;
std::cout << "subdivisoins: "<<subdiv<<std::endl;
std::cout << "try fatarc generation: "<<arcgen<<std::endl;
//std::cout << "solution: "<<mo.size()<<std::endl;
   

 while(!mo.empty()){

 box_t rb=mo.top(); 
 moseq<< rb;
 mo.pop();
 };

 return moseq;

 
 }/*solver*/


};//end of structure solver_mv_fatarcs

//                            TEST


template<class C>
  Seq< domain<C> >  solve_bv_fatarcs( polynomial< C, with<Bernstein> > p1, 
                      polynomial< C, with<Bernstein> > p2, C eps  ) 
{
  int MTH=1;
  solver_mv_fatarcs<C> sys( p1, p2, eps);
  return sys.solver(MTH);
}