// plygn.cpp
//
// 2-dim convex polygonal domain
// 
// Jiangguo (James) Liu
// ColoState, Jan.-Jun. 2007


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


// Plygn: copy constructor

Plygn::Plygn(const Plygn &plygn)
{
   numVrtz = plygn.numVrtz;
   vrtx = new PtVec2d[numVrtz+1];  // +1, new memory is allocated
   for (int i=0; i<=numVrtz; ++i)  vrtx[i] = plygn.vrtx[i];
}


// Plygn: constructor based on a given data file

Plygn::Plygn(char *filename)
{
   FILE *fp;  // File pointer
   double x, y;

   if (NULL==(fp=fopen(filename, "r"))) {
      puts("Open data file failed.  Exit.");
      exit(-1);
   }

   fscanf(fp, "%2d", &numVrtz);

   vrtx = new PtVec2d[numVrtz+1];  // Don't forget +1

   for (int i=0; i<numVrtz; ++i) {
      fscanf(fp, "%lf  %lf", &x, &y);
      vrtx[i] = PtVec2d(x,y);
   }

   fclose(fp);

   vrtx[numVrtz] = vrtx[0];
}


// Plygn: the center of a polygon

PtVec2d Plygn::center() const 
{
   PtVec2d P(0,0);
   for (int i=0; i<numVrtz; ++i)  P = P + vrtx[i];
   return (1.0/numVrtz)*P;
}


// Plygn: the area of a polygon

double Plygn::area() const 
{
   double a = 0.0;
   for (int i=0; i<numVrtz; ++i) {
      Trig T(center(), vrtx[i], vrtx[i+1]);
      a += T.area();
   }
   return a;
}


// Plygn: is a given point inside the polygon
// 0:  outside
// +:  the label of the decomposed triangle the point is in

int Plygn::isInPlygn(PtVec2d P) const 
{
   int j = 0;
   for (int i=0; i<numVrtz; ++i) {
      Trig T(center(), vrtx[i], vrtx[i+1]);
      if (T.isInTrig(P)) {
         j = i+1;
         break;
      }
   }
   return j;
}


// Plygn: boundary intersection 
// Assume Po outside, Pi inside 
// Return Pc on the boundary

PtVec2d Plygn::bndryInterSxn(const PtVec2d &Pi, const PtVec2d &Po) const 
{
   int j = isInPlygn(Pi);
   Trig T(center(), vrtx[j-1], vrtx[j]);
   return T.hypInterSxn(Pi,Po);
}


// Update the inner domain during each time step

void updateInnerDom(Plygn &innerDom, const Plygn &realDom, 
   double vmax, double told, double tnew)
{
   double distance = vmax*(tnew-told);
   for (int i=0; i<innerDom.numVrtz; ++i) {
      PtVec2d Pi = realDom.vrtx[i];
      PtVec2d Pi1 = realDom.vrtx[i+1];
      double side = dist(Pi, Pi1);
      Trig T(realDom.center(), realDom.vrtx[i], realDom.vrtx[i+1]);
      double height = 2*T.area()/side;
      double alpha = distance/height;
      innerDom.vrtx[i] = alpha*realDom.center() + (1-alpha)*realDom.vrtx[i];
   }
   innerDom.vrtx[innerDom.numVrtz] = innerDom.vrtx[0];
   return;
}