gsTraceCurve.hpp

Go to the documentation of this file.

 
#pragma once
 
#include <iostream>
#include <gsCore/gsLinearAlgebra.h>
#include <gsCore/gsFunction.h>
#include <gsCore/gsFunctionExpr.h>
#include <gsNurbs/gsBSplineBasis.h>
#include <gsUtils/gsPointGrid.h>
 
 
namespace gismo
{
 
template<class T>
void gsTraceLine( std::pair<gsFunction<T>*,gsFunction<T>*>  & map,
                       gsMatrix<T> const & x,gsMatrix<T> const & p, gsMatrix<T> & result)
{
    T h;
    int loops=10;
    int n_points = 10;
    T tolerance = 0.0001;
 
     gsMatrix<T> line(2,1), px, py;
     gsMatrix<T> next (2,1), delta;
     gsMatrix<T> tangent (2,1) ;
     gsMatrix<T> current (2,1);
     gsMatrix<T> Jmap (2,2);
 
 
 
     map.first->eval_into(x, px);
     map.second->eval_into(x, py);
 
 
 
  // matrix \line stores the evaluation of the \map at the current point
     line.row(0)= px;
     line.row(1)= py;
     next.row(0)= px;
     next.row(1)= py;
 
     tangent = p.col(0)-line.col(0);
 
      h = tangent.norm()/n_points ;
 
     tangent/=tangent.norm();
 
 
  // "current" stores the point we start with
     current = x.col(0);
 
   int count=0;
 
   for(int i=0; i!=n_points; i++)
   { next = next + h*tangent;
      loops=10;
 
  do
  {
      count++;
      //gsDebug<< "Point "<<count <<"= "<< current.transpose() <<"\n";
 
      Jmap.row(0) =  map.first->jacobian(current );
      Jmap.row(1) =  map.second->jacobian(current);
 
      if ( math::abs( Jmap.determinant() ) < 0.0001) //0.025 amoeba_hole,austria_hole
            {
               gsDebug<< "\n trace line: \n Jacobian vanished at : " << current.transpose() <<"\n"
                   "with Jacobian = "<< math::abs( Jmap.determinant() )<<"\n \n";
 
               loops =-1 ;
                break; // stop Newton iteration
            }
 
      delta = Jmap.inverse()*(next.col(0)-line.col(0));
      current += delta;
 
      map.first->eval_into(current,px);
      map.second->eval_into(current,py);
 
      line.row(0)=px;
      line.row(1)=py;
 
      // gsDebug<<"loops="<<loops<<"\n";
 
  }
  while((--loops > 0) && ( (line-next).norm() > tolerance  ) );
 
 
  result = current;
 
  // gsDebug<<"points traced : "<< result<< "\n";
 
 
 
    // example:
    // gsFunctionExpr<T> g1("x");
    // gsFunctionExpr<T> g2("y");
    // std::pair<gsFunction<T>*,gsFunction<T>*> map(g1,g2);
    // gsMatrix<T> res, p(2,2);
    // p<< 0,1,
    //     0,1 ;
    // gsTraceCurve(map, p, res);
 
 
 
}
 
}
 
 
template<class T>
void gsTraceCurve( std::pair<gsFunction<T>*,gsFunction<T>*>  & map,
                       gsMatrix<T> const & x, gsBSpline<T> * bs,gsMatrix<T> & result,const int n_points = 50,
                       const T tolerance = 0.0001)
{
 
    int loops=10;// usually are 100
 
 
    result.resize(2,n_points-1);
 
       gsMatrix<T> current (2,1);
       gsMatrix<T> next (2,1), delta;
       gsMatrix<T> line(2,1), px, py;
       gsMatrix<T> Jmap (2,2);
 
      // to do: change with bs->parameterRange()
      gsMatrix<T> pts =  gsPointGrid<T>(0,1,n_points);
 
      gsMatrix<T> ev =  bs->eval(pts);
 
      //next = ev.col(0); // ev.col(0)==0 point removed!
      next= ev.col(1);
 
      gsMatrix<T> rr; // rr stores the initial point P0 from which the tracing starts
      gsTraceLine( map, x, next, rr );
 
 
 
       current = rr.col(0);
 
 
       result.col(n_points-2)= current;
 
 
         next=ev.col(1);
 
     map.first->eval_into(current, px);
     map.second->eval_into(current, py);
 
  // matrix \line stores the evaluation of the \map at the current point
     line.row(0)= px;
     line.row(1)= py;
 
   int count=0;
 
   for(int i=2; i!=n_points; i++)
   {
      next = ev.col(i);
   loops=10;
  do
  {
 
  count++;
 
 
  // TO DO: transpose the result of  gsFunction::grad and add transpose()
  Jmap.row(0) =  map.first->jacobian(current);
  Jmap.row(1) =  map.second->jacobian(current);
 
 
   delta = Jmap.inverse()*(next.col(0)-line.col(0));
 
   current += delta;
 
   //   gsMatrix<T> * grad1 = map.first->deriv() ;
   //   Jmap.row(0) = *grad1;
   //   delete grad1;
 
   map.first->eval_into(current,px);
   map.second->eval_into(current,py);
 
   line.row(0)=px;
   line.row(1)=py;
 
   //gsDebug<<"\n loops = "<<loops<<"\n";
 
 
  }
  while((--loops > 0) && ( (line-next).norm() > tolerance  ) );
 
 
  result.col(n_points-i-1) = current;
// gsDebug<<"\n result.col("<<i<<") : \n" << result.col(i)<<"\n";
 
 
 
 
    // example:
    // gsFunctionExpr<T> g1("x");
    // gsFunctionExpr<T> g2("y");
    // std::pair<gsFunction<T>*,gsFunction<T>*> map(g1,g2);
    // gsMatrix<T> res, p(2,2);
    // p<< 0,1,
    //     0,1 ;
    // gsTraceCurve(map, p, res);
 
 
 
}
 
 
}
 
template<class T>
void gsTraceCurvePart(std::pair<gsFunction<T>*,gsFunction<T>*>  & map,
                      gsMatrix<T> const & x, gsBSpline<T> * bs, T t0, T t1,
                      gsMatrix<T> & result,const int n_points = 50,
                      const T tolerance = 0.0001)
{
    int loops=100;
 
    gsMatrix<T> current (2,1);
    current=x;
    gsMatrix<T> next (2,1), delta, rr;
    gsMatrix<T> line(2,1), px, py;
    gsMatrix<T> Jmap (2,2);
    gsMatrix<T> medium_point (1,1);
    medium_point<< 0.5 ;
 
    gsMatrix<T> discretI = gsPointGrid<T>(t0,t1, n_points );
    gsMatrix<T> ev = bs->eval(discretI);
    map.first->eval_into(current, px);
    map.second->eval_into(current, py);
    line.row(0)= px;
    line.row(1)= py;
 
    result.resize(2,ev.cols()+1);
    result.col(0)=x;
    bool cont = true;
    index_t i;
    for(i =1; i< ev.cols()&& cont; i++)
    {
        next= ev.col(i);
        loops=100;
        do
        {
            // TO DO: transpose the result of  gsFunction::grad and add transpose()
            Jmap.row(0) =  map.first->jacobian(current);
            Jmap.row(1) =  map.second->jacobian(current);
 
            //gsDebug<<"jacobian matrix is \n"<< Jmap<<"\n";
            if ( math::abs( Jmap.determinant() ) < 0.01) //0.025 amoeba_hole,austria_hole
            {
                gsDebug<< "Jacobian vanished at : " << current.transpose() <<"\n"
                "with Jacobian = "<< math::abs( Jmap.determinant() )<<"\n \n";
 
                cont = false;
                break; // stop Newton iteration
            }
 
            delta = Jmap.inverse()*(next.col(0)-line.col(0));
 
            current += delta;
 
            map.first->eval_into(current,px);
            map.second->eval_into(current,py);
 
            line.row(0)=px;
            line.row(1)=py;
 
 
 
        }
        while((--loops > 0) && ( (line-next).norm() > tolerance  ) );
 
        result.col(i) = current;
    }
 
  // gsDebug<< "Size expected:"<< result.cols() <<"\n";
    result.conservativeResize( gsEigen::NoChange, i);
    //gsDebug<< "Size finally :"<< result.cols()  <<"\n";
 
 
 
 
 
}
 
} // namespace gismo