// CharaTrkn2d.cpp
//
// Subroutines for characteristic tracking in 2-dim
//
// Jiangguo (James) Liu
// ColoState, Jan. 2007


#include <cmath>
#include <iostream>
#include <stdio.h>
#include "AuxMesh2d.h"
#include "chara2d.h"
#include "plygn.h"
#include "PtVec2d.h"
#include "RealMesh2d.h"
using namespace std;


// Chara2d: constructor

Chara2d::Chara2d(int maxNumPtsChara)
{
   maxNumPts = maxNumPtsChara;
   tm = new double[maxNumPts];
   pt = new PtVec2d[maxNumPts];
}


// Chara2d: resize

void Chara2d::resize(int maxNumPtsChara)
{
   delete[] tm,pt;
   maxNumPts = maxNumPtsChara;
   tm = new double[maxNumPts];
   pt = new PtVec2d[maxNumPts];
}


// Backward tracking, Euler method

void bkwdTrknEuler(Chara2d &cha, 
   const Plygn &realDom, const Plygn &innerDom, 
   const RealMesh2d &realMesh, const AuxMesh2d &auxMesh, 
   PtVec2d (*vel)(PtVec2d), double told, double tnew, 
   double thead, PtVec2d head) 
{
   int i;
   double alpha, deltat, dt;
   PtVec2d P0, P1, Pc, v;

   deltat = (tnew-told)/(cha.maxNumPts-1);
   for (i=0; i<cha.maxNumPts; ++i)  cha.tm[i] = told + i*deltat;

   cha.headOnBndry = 0;
   cha.footOnBndry = 0;

   cha.thead = thead;
   cha.head = head;

   if (innerDom.isInPlygn(head)) {
      cha.indHead = cha.maxNumPts-1;
      cha.tm[cha.indHead] = thead;
      cha.pt[cha.indHead] = head;
      P1 = head;
      for (i=cha.indHead; i>=1; --i) {
         v = vel(P1);
         P0 = P1 - deltat*v;
         cha.pt[i-1] = P0;
         P1 = P0;
      }
      cha.tfoot = told;
      cha.foot = P0;
      cha.indFoot = 0;
      cha.footOnBndry = 0;
   }

   if (!innerDom.isInPlygn(head)) {
      cha.indHead = (int)ceil((thead-told)/deltat);
      cha.tm[cha.indHead] = thead;
      cha.pt[cha.indHead] = head;
      P1 = head;
      for (i=cha.indHead; i>=1; --i) {
         v = vel(P1);
         dt = cha.tm[i]-cha.tm[i-1];
         P0 = P1 - dt*v;
         if (realDom.isInPlygn(P0)) {
            cha.pt[i-1] = P0;
            P1 = P0;
            cha.tfoot = cha.tm[i-1];
            cha.foot = P0;
            cha.footOnBndry = 0;
            cha.indFoot = i-1;
            continue;
         }
         else {
            Pc = realDom.bndryInterSxn(P1, P0);
            alpha = dist(Pc,P1)/dist(P0,P1);
            cha.tfoot = cha.tm[i-1]+alpha*dt;
            cha.foot = Pc;
            cha.indFoot = i-1;
            cha.footOnBndry = 1;
            cha.tm[cha.indFoot] = cha.tfoot;
            cha.pt[cha.indFoot] = cha.foot;
            break;
         }
      }
      cha.Deltat = cha.tm[cha.indHead]-cha.tm[cha.indFoot];
      cha.Jacobian = 1.0;
   }

   cha.Deltat = cha.tm[cha.indHead]-cha.tm[cha.indFoot];
   cha.Jacobian = 1.0;
   return;
}


// Backward tracking, 2nd-order Runge-Kutta method

void bkwdTrknRK2(Chara2d &cha, 
   const Plygn &realDom, const Plygn &innerDom, 
   const RealMesh2d &realMesh, const AuxMesh2d &auxMesh, 
   PtVec2d (*vel)(PtVec2d P), double told, double tnew, 
   double thead, PtVec2d head)
{
   int i;
   double deltat;
   PtVec2d P0, P1, F0, F1;

   deltat = (tnew-told)/(cha.maxNumPts-1);
   for (i=0; i<cha.maxNumPts; ++i)  cha.tm[i] = told + i*deltat;

   cha.headOnBndry = 0;
   cha.footOnBndry = 0;

   cha.thead = thead;
   cha.head = head;

//   if (innerDom.isInPlygn(head)) {
      cha.indHead = cha.maxNumPts-1;
      cha.tm[cha.indHead] = thead;
      cha.pt[cha.indHead] = head;
      P1 = head;
      for (i=cha.indHead; i>=1; --i) {
         F1 = deltat*vel(P1);
         P0 = P1 - F1;
         F0 = deltat*vel(P0);
         P0 = P1 - 0.5*(F1 + F0);
         cha.pt[i-1] = P0;
         P1 = P0;
      }
      cha.tfoot = told;
      cha.foot = P0;
      cha.indFoot = 0;
      cha.footOnBndry = 0;
      cha.Deltat = cha.thead - cha.tfoot;
      cha.Jacobian = 1.0;
//   }

   return;
}