// main.cpp
//
// J.Liu, R.Cali
// ColoState, Jan.-Jun. 2007

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

#ifndef M_PI
#define M_PI 3.141592653589793
#endif  // M_PI

int main() 
{
   FILE *fp;
   char fn[10], str[4];
   int eltType, i, numHitsBndry=0, numVrtz;
   int *lblVrtx=0;
   double ratio, tfoot, thead;
   PtVec2d foot, head, P, V;
   PtVec2d *vrtx=0;

   RealMesh2d realMesh("TrigMesh.dat");
   cout << "Real mesh max size= " << realMesh.maxMeshSize() << "\n";
   cout << "Real mesh min size= " << realMesh.minMeshSize() << "\n";
   cout << "\n";
   realMesh.saveMesh();

   AuxMesh2d auxMesh(realMesh);
   cout << "avg. num. trigs= " << auxMesh.avgNumEltsPerRect() << "\n";
   auxMesh.saveRectVsEltsData();
   cout << "\n";

   Plygn realDom("RealDom.dat");
   Plygn innerDom(realDom);

   double T = 1.0;
   int NT = 20;
   double Deltat = T/(double)NT;
   double vmax = sqrt(2.0)/2;

   Chara2d cha(16);
   for (int n=14; n<15; ++n) {
      double told = n*Deltat;
      double tnew = (n+1)*Deltat;
      updateInnerDom(innerDom, realDom, vmax, told, tnew);

      cout << "n=" << n << "\n";
      strcpy(fn, "Ut");
      sprintf(str, "%d", n); 
      strcat(fn, str);
      strcat(fn, ".dat");

      if (NULL==(fp=fopen(fn, "w"))) {
         puts("Open Utxxx.dat file failed.  Exit.");
         exit(-1);
      }

      fprintf(fp, "n=%2d\n", n);

      for (int ie=realMesh.bgnLblElt(); ie<=realMesh.endLblElt(); ++ie) {
         lblVrtx = realMesh.getEltInfo(eltType, ie);
         thead = tnew;

         if (eltType==trig) {  // triangle
            numVrtz = 3;
            vrtx = new PtVec2d[numVrtz];
            for (i=0; i<numVrtz; ++i)  vrtx[i] = realMesh.node(lblVrtx[i]);
            Trig T(vrtx[0],vrtx[1],vrtx[2]);
            delete[] vrtx;
            head = T.center();
         }

         bkwdTrknEuler(cha, realDom, innerDom, realMesh, auxMesh, 
            vel2dExRotating, told, tnew, thead, head);
         cha.getFoot(tfoot, foot);

        if (cha.isFootOnBndry()) {
           cout << "ie=" << ie << ", ";
           cout << "head=" << head;
           cout << "foot=" << foot << "\n";
        }

         numHitsBndry += cha.isFootOnBndry();
         int loc = locInRealMesh(foot, realMesh, auxMesh);
      }

      ratio = (double)numHitsBndry/(double)realMesh.numElts();
      cout << "numHitsBndry= " << numHitsBndry << "\n";
      cout << "ratio= " << ratio << "\n";

      fclose(fp);
   }

   cout << "\nEvery thing is fine!\n";
   return 0;
}