#include <iostream>
#include <fstream>
#include <cstdio>
#include <cstdlib>

#include <map>

#include "global.h"
#include "mesh3d.h"
#include "n2f_dd.h"

void makeTriArray(int &numAbcTri, tri3d ***triArray, char *meshName )
{
  int numNodes;

  std::map<edge3d, edge3d*> edgemap;

  char triName[180];
  sprintf(triName, "%s.tri", meshName);

  ifstream foo(triName, ios::in);
  if (!foo) {
    cout << "could not find file " << triName << endl;
    exit(0);
  }

  foo >> unitConv;
  // read in nodes
  foo >> numNodes;
  //  cout << "numNodes = " << numNodes << endl;
  node3d **nodeArray = new node3d*[numNodes];
  int i;
  double x, y, z;
  for (i = 0; i < numNodes; i ++) {
    foo >> x >> y >> z;

    x *= unitConv;
    y *= unitConv;
    z *= unitConv;

	node3d *nd = new node3d();
    nd->SetCoord(x, y, z);
	nd->SetId(i);
    nodeArray[i] = nd;
  }

  // read in the triangles
  foo >> numAbcTri;
  cout << "numAbcTri = " << numAbcTri << endl;
  (*triArray) = new tri3d*[numAbcTri];
  node3d *node[3];
  tri3d tri, *triPtr;
  edge3d *egPtr[3];
  int j, id;
  for (i = 0; i < numAbcTri; i ++) {
    // read in the three nodes for the triangle
    for (j = 0; j < 3; j ++) {
      foo >> id;
      node[j] = nodeArray[id];
    }
    // insert the three edges into the tree
    for (j = 0; j < 3; j ++) {
      int j1, j2;
      j1 = (j + 1) % 3;
      j2 = (j + 2) % 3;

	  edge3d eg;
      eg.SetEdge_no_arrange(node[j1], node[j2]);
	  std::map<edge3d,edge3d*>::iterator f = edgemap.find(eg);
	  if (f != edgemap.end())
		  egPtr[j] = f->second;
	  else
		  egPtr[j] = edgemap[eg] = new edge3d(eg);
    }
    // cross-link between triangle and edges
    (*triArray)[i] = new tri3d;
    (*triArray)[i]->SetTri(egPtr[0], egPtr[1], egPtr[2]);
    (*triArray)[i]->setNum( i );
  }

  delete[] nodeArray;
  foo.close();
}

void setRegister(regcvtr **regCur, int numAbcTri,
		 tri3d **abcTriArray, char *curName, int cnt)
{
  ifstream foo(curName, ios::in);
  if (!foo) {
    cout << "could not find file " << curName << endl;
    exit(0);
  }

  int i, j, k;
  int m, n;
  tri3d *fc;
  double Xcoord, Ycoord, Zcoord;                   
  vtr coord[6];                                   
  cVtr JorM;
  Complex JorM_x, JorM_y, JorM_z;
  double re, im;

  *regCur = new regcvtr[numAbcTri];

  for (i = 0; i < numAbcTri; i ++) {
    int tetID, faceIndex;
    vtr normal;
    double nx, ny, nz;

    (*regCur)[i].init(cnt);

    fc = abcTriArray[i];
    for (j = 0; j < cnt; j ++) {
      foo >> re >> im;
      JorM_x = Complex(re, im);
      foo >> re >> im;
      JorM_y = Complex(re, im);
      foo >> re >> im;
      JorM_z = Complex(re, im);
      JorM.setcvtr(JorM_x, JorM_y, JorM_z);
      (*regCur)[i].SetField(j, JorM);
    }
  }
}
void setRegister1(regcvtr **regCur, int numAbcTri,
		  tri3d **abcTriArray, char *curName, vtr k, cVtr HorEinc, 
		  int opt )
{
  ifstream foo(curName, ios::in);
  k = k * unitConv;
  int i, j;
  int m, n;
  tri3d *fc;
  double Xcoord, Ycoord, Zcoord;
  double nx, ny, nz;
  vtr coord[6], normal;
  cVtr JorM;
  cVtr JorMinc;
  Complex JorM_x, JorM_y, JorM_z;
  double re, im;
  double const static ZERO_7 = 1.e-7;

  char null = 0;
  if (*regCur == NULL) null = 1;
  
  if (null) 
    *regCur = new regcvtr[numAbcTri];
  for (i = 0; i < numAbcTri; i ++) {
    if (null) (*regCur)[i].init(3);

    if(foo.is_open()==false)continue;

    fc = abcTriArray[i];
    fc->GetNormal(&nx, &ny, &nz);

    normal.setvtr( nx, ny, nz );
    for (j = 0; j < 3; j++) {
      Xcoord = (fc->GetNodePtr(j))->GetX();
      Ycoord = (fc->GetNodePtr(j))->GetY();
      Zcoord = (fc->GetNodePtr(j))->GetZ();
      coord[j] = vtr(Xcoord, Ycoord, Zcoord);
    }
    coord[3] = (coord[1] + coord[2]) * 0.5;
    coord[4] = (coord[0] + coord[2]) * 0.5;
    coord[5] = (coord[0] + coord[1]) * 0.5;

    for (j = 0; j < 3; j ++) {
      foo >> re >> im;
      JorM_x = Complex(re, im);
      foo >> re >> im;
      JorM_y = Complex(re, im);
      foo >> re >> im;
      JorM_z = Complex(re, im);

      JorM.setcvtr(JorM_x, JorM_y, JorM_z);
 //      JorMinc = HorEinc * exp(Complex(0.0, -dotP(k,coord[j]))) * normal;
      JorMinc = HorEinc * exp(Complex(0.0, -dotP(k,coord[j])));
      if (opt == 1)
	JorM = JorM - JorMinc;
      (*regCur)[i].SetField(j, JorM);
    }
  }
}

void ComputeTriangleNormals(int triCnt, tri3d **triArrayPtr)
{
  // first simply the normal for every triangle and reset the done flag
  int i;
  vtr ndVtr;
  for (i = 0; i < triCnt; i ++) {
    triArrayPtr[i]->ComputeNormal();
    triArrayPtr[i]->ResetDone();
  }

  triArrayPtr[0]->SetDone();
  tri3d *neigTri;
  for (i = 0; i < 3; i ++) {
    neigTri = GetNeighborTriPtr(triArrayPtr[0], i);
    
    if (neigTri == 0) continue;
    if (neigTri->IsDone() != 1) CheckNormal(neigTri);
  }
}


void CheckNormal(tri3d *tt)
{
  if (tt == 0) return;
  if (tt->IsDone()) return;
  
  int id;
  tri3d *doneNeighbor;
  doneNeighbor = FindDoneNeighborTriPtr(tt, &id);
  if (doneNeighbor == 0) return;

  if (Consistent(tt, doneNeighbor, id) == 0) tt->NormalReverse();
  tt->SetDone();

  int i;
  tri3d *neigTri;
  for (i = 0; i < 3; i ++) {
    neigTri = GetNeighborTriPtr(tt, i);
    
    if (neigTri == 0) continue;
    if (neigTri->IsDone() != 1) CheckNormal(neigTri);
  }
}


tri3d* FindDoneNeighborTriPtr(tri3d *tt, int *id)
{
  int i;
  for (i = 0; i < 3; i ++) {
    tri3d *neigTri = GetNeighborTriPtr(tt, i);
    if (neigTri == 0) continue;

    if (neigTri->IsDone()) {
      (*id) = i;
      return neigTri;
    }
  }

  return 0;
}


int Consistent(tri3d *t1, tri3d *t2, int id)
{
  edge3d *eg;
  int pol1, pol2;

  eg = t1->GetEdgePtr(id);
  pol1 = Polarity(t1, eg);
  pol2 = Polarity(t2, eg);
  if (pol1 == pol2) return 0;

  return 1;
}


int Polarity(tri3d *tt, edge3d *eg)
{
  node3d *nd, *nnd, *nnd1;

  nd = FindOppositeNode(tt, eg);
  if (nd == 0) return 0;

  nnd = eg->GetNodePtr(0);
  nnd1 = eg->GetNodePtr(1);
  vtr avtr, bvtr, nvtr;
  double vx, vy, vz;
  tt->GetNormal(&vx, &vy, &vz);
  nvtr.setvtr(vx, vy, vz);
  avtr.setvtr(nnd1->GetX() - nnd->GetX(), nnd1->GetY()-nnd->GetY(),
	      nnd1->GetZ() - nnd->GetZ());
  bvtr.setvtr(nd->GetX()-nnd->GetX(), nd->GetY()-nnd->GetY(),
	      nd->GetZ()-nnd->GetZ());
  if (dotP((nvtr * avtr), bvtr) > 0.0) return 1;

  return -1;
}


node3d *FindOppositeNode(tri3d *tt, edge3d *eg)
{
  int i;

  i = 0;
  do {
    if (tt->GetEdgePtr(i) == eg) return tt->GetNodePtr(i);

    i ++;
  } while (i < 3);

  return 0;
}