// AuxMesh2d.cpp
//
// Jiangguo (James) Liu
// ColoState, Jan.-Jun. 2007


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


// AuxMesh2d: a more meaningful constructor
// Creating the auxiliary mesh from the real mesh

AuxMesh2d::AuxMesh2d(const RealMesh2d &realMesh)
{
   int i, id, ie, ir, j, k, kk, loc;
   int eltType;
   int numVrtz, *lblVrtx=0;
   double xmax, xmin, xtmp;
   double ymax, ymin, ytmp;
   PtVec2d *vrtx=0;

// Finding the minimal covering rectangle 

   DomA = 0.0;
   DomB = 0.0;
   DomC = 0.0;
   DomD = 0.0;

   for (id=realMesh.bgnLblNd(); id<=realMesh.endLblNd(); ++id) {
      xtmp = (realMesh.node(id)).xCrd();
      ytmp = (realMesh.node(id)).yCrd();
      if (xtmp<DomA)  DomA = xtmp;
      if (xtmp>DomB)  DomB = xtmp;
      if (ytmp<DomC)  DomC = ytmp;
      if (ytmp>DomD)  DomD = ytmp;
   }
/*
   cout << "DomA= " << DomA << "\n";
   cout << "DomB= " << DomB << "\n";
   cout << "DomC= " << DomC << "\n";
   cout << "DomD= " << DomD << "\n";
*/

// Creating the auxiliary mesh

   nx = (int)ceil((DomB-DomA)/realMesh.minMeshSize());
   ny = (int)ceil((DomD-DomC)/realMesh.minMeshSize());
   nr = nx*ny;

   Deltax = (DomB-DomA)/nx;
   Deltay = (DomD-DomC)/ny;

   x = new double[nx+1];
   y = new double[ny+1];

   x[0] = DomA;
   y[0] = DomC;

   for (i=1; i<=nx; ++i)  x[i] = x[i-1] + Deltax;
   for (j=1; j<=ny; ++j)  y[j] = y[j-1] + Deltay;

   x[nx] = DomB;
   y[ny] = DomD;

   // cout << "Auxiliary mesh created\n";
   // cout << "nx=" << nx << "\n";
   // cout << "ny=" << ny << "\n";
   // printf("x[nx]=%48.45f\n", x[nx]);

// Creating the connection table: element versus rectangular boxes

   int ne = realMesh.numElts();

   int *nre = new int[ne];  // number of rect. boxes per element
   int **er = new int *[ne];  // element vs. rectangular box

   for (ie=0; ie<ne; ++ie)  nre[ie] = 0;  // Number 0
   for (ie=0; ie<ne; ++ie)  er[ie] = 0;  // Null pointer

   ner = new int[nr];
   for (ir=0; ir<nr; ++ir)  ner[ir] = 0;  // Number 0

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

      if (eltType==rect2d) {  // rectangle in 2-dim
         numVrtz = 2;
         vrtx = new PtVec2d[numVrtz];
         vrtx[0] = realMesh.node(lblVrtx[0]);
         vrtx[1] = realMesh.node(lblVrtx[1]);
         Rect2d R(vrtx[0], vrtx[1]);
         delete[] vrtx;
         R.getMaxMinCrds(xmax, xmin, ymax, ymin);
      }

      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;
         T.getMaxMinCrds(xmax, xmin, ymax, ymin);
      }
/*
      if (eltType==quadri) {  // quadrilateral
         numVrtz = 4;
         vrtx = new PtVec2d[numVrtz];
         for (i=0; i<numVrtz; ++i)  vrtx[i] = realMesh.node(lblVrtx[i]);
         Quadri Q(vrtx[0], vrtx[1], vrtx[2], vrtx[3]);
         delete[] vrtx;
         Q.getMaxMinCrds(xmax, xmin, ymax, ymin);
      }
*/
      if (xmin<DomA)  xmin = DomA;
      if (xmax>DomB)  xmax = DomB;
      if (ymin<DomC)  ymin = DomC;
      if (ymax>DomD)  ymax = DomD;

      int Imin = (int)ceil((xmin-DomA)/Deltax);
      int Imax = (int)ceil((xmax-DomA)/Deltax);
      int Jmin = (int)ceil((ymin-DomC)/Deltay);
      int Jmax = (int)ceil((ymax-DomC)/Deltay);

      if (Imin<1)  Imin=1;
      if (Imax>nx)  Imax = nx;
      if (Jmin<1)  Jmin=1;
      if (Jmax>ny)  Jmax = ny;

      loc = ie-realMesh.bgnLblElt();
      nre[loc] = (Imax-Imin+1)*(Jmax-Jmin+1);
      er[loc] = new int[nre[loc]];

      k = 0;
      for (i=Imin; i<=Imax; ++i) {
         for (j=Jmin; j<=Jmax; ++j) {
            kk = (i-1)*ny + j-1;
            er[loc][k] = kk;
            ner[kk] += 1;
            k += 1;
         }
      }
   }

// Inverting the connection table to create the look-up table: 
// rectangular box versus elements

   re = new int *[nr];
   for (ir=0; ir<nr; ++ir)  re[ir] = new int[ner[ir]];

   for (ir=0; ir<nr; ++ir)  ner[ir] = 0;  // a little bit weird

   for (ie=0; ie<ne; ++ie) {
      for (j=0; j<nre[ie]; ++j) {
         kk = er[ie][j];
         re[kk][ner[kk]] = ie + realMesh.bgnLblElt();
         ner[kk] += 1;
      }
   }

// Cleanning up

   for (ie=0; ie<ne; ++ie)  delete[] er[ie];
   delete[] er, nre;

   return;
}


// AuxMesh2d: destructor

AuxMesh2d::~AuxMesh2d() 
{
   delete[] x, y;
   delete[] ner;
   for (int i=0; i<nr; ++i)  delete[] re[i];
   delete[] re;
}


// AuxMesh2d: Average number of elements per rectangle

float AuxMesh2d::avgNumEltsPerRect() const 
{
   int m = 0;
   for (int i=0; i<nr; ++i)  m += ner[i];
   return (float)m/(float)nr;
}


// AuxMesh2d: Numeber of elements intersecting with rectangle box ind

int AuxMesh2d::numEltsWithRect(Index2 ind) const 
{
   Index2 MaxInd(nx,ny);
   int lind = lindex(ind, MaxInd);
   return ner[lind-1];
}


// AuxMesh2d: Look-up table: 
// rectangular box ind versus j-th element 
// which intersects with the recangular box

int AuxMesh2d::rectVsElt(Index2 ind, int j) const 
{
   Index2 MaxInd(nx,ny);
   int lind = lindex(ind, MaxInd);
   // cout << "lind,j=" << lind << " , " << j << "\n";
   return re[lind-1][j-1];
}


// AuxMesh2d: 
// Saving the look-up table to a data file
// Fixed filename: "RectVsElts.dat"

void AuxMesh2d::saveRectVsEltsData() const 
{
   FILE *fp;

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

   // fprintf(fp, "DomABCD = %7.4f  %7.4f  %7.4f  %7.4f\n", DomA, DomB, DomC, DomD);
   // fprintf(fp, "Rectangular boxes = %7d  %7d  %10d\n", nx, ny, nr);
   fprintf(fp, "%10d\n", nr);

   for (int ix=1; ix<=nx; ++ix) {
      for (int iy=1; iy<=ny; ++iy) {
         int i = (ix-1)*ny + iy;
         // fprintf(fp, "(%7d, %7d)  %4d:  ", ix, iy, ner[i-1]);
         fprintf(fp, "%7d  %4d  ", i, ner[i-1]);
         for (int j=1; j<=ner[i-1]; ++j) {
            fprintf(fp, "%7d", re[i-1][j-1]);
         }
         fprintf(fp, "\n");
      }
   }

   fclose(fp);
   return;
}


// AuxMesh2d: Locating a point in the auxiliary rectangular mesh

Index2 AuxMesh2d::locInAuxMesh(PtVec2d P) const 
{
   double Px = P.xCrd();
   double Py = P.yCrd();

   int ix=0;
   for (int i=1; i<=nx; ++i) {
      if (x[i-1]<=Px && Px<=x[i]) {
         ix = i;
         break;
      }
   }

   int iy=0;
   for (int j=1; j<=ny; ++j) {
      if (y[j-1]<=Py && Py<=y[j]) {
         iy = j;
         break;
      }
   }

   return Index2(ix,iy);
}


// Locating a point in the real mesh 
// via the look-up table in the aux. mesh

int locInRealMesh(PtVec2d P, const RealMesh2d &realMesh, 
   const AuxMesh2d &auxMesh)
{
   int loc;
   int eltType;
   int i, j, k, numVrtz;
   int *lblVrtx=0;
   PtVec2d *vrtx=0;

   Index2 ind = auxMesh.locInAuxMesh(P);

   loc = -1;
   for (j=1; j<=auxMesh.numEltsWithRect(ind); ++j) {
      k = auxMesh.rectVsElt(ind, j);
      lblVrtx = realMesh.getEltInfo(eltType, k);

      if (eltType==rect2d) {  // rectangle in 2-dim
         numVrtz = 2;
         vrtx = new PtVec2d[numVrtz];
         vrtx[0] = realMesh.node(lblVrtx[0]);
         vrtx[1] = realMesh.node(lblVrtx[1]);
         Rect2d R(vrtx[0], vrtx[1]);
         delete[] vrtx;
         if (R.isInRect2d(P))  loc = k;
      }

      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;
         if (T.isInTrig(P))  loc = k;
      }
   }

   return loc;
}