// TetraMesh.cpp
// Class for tetrahedral meshes
// James Liu, Graham Harper, ColoState; 2007/01--2018/05

#include <iostream>
#include <stdio.h>
#include <stdlib.h>

#include "cell3d.h"
#include "PtVec3d.h"
#include "TetraMesh.h"


// TetraMesh: Default constructor

TetraMesh::TetraMesh() 
{
  NumNds = 0;  NumEms = 0;  NumFcs = 0;  NumBndryPcs = 0;
  BgnLblNd = 1;  BgnLblEm = 1;  BgnLblFc = 1;  BgnLblBndry = 1;  // by default
  IsBndryNd = 0;  NumNdsNd = 0;  NumEmsNd = 0;  NumFcsNd = 0;
  LblNdNd = 0;  LblNdEm = 0;  LblNdFc = 0;
  IsBndryFc = 0;  LblFcNd = 0;  LblFcEm = 0;  LblFcFc = 0;
  LblEmNd = 0;  LblEmFc = 0;  SignEmFc = 0;
  nd = 0;  OutNmlBndry = 0;
  flag = 0;
}


// TetraMesh: A (more meaningful) constructor:
// generating a uniform tetrahedral mesh with provided parameters

TetraMesh::TetraMesh(double xa, double xb, int nx,
                     double yc, double yd, int ny,
                     double ze, double zf, int nz)
{
  int i, j, k, l, labelc, labele, lc, ld, le, lea, leb, m, n;
  int labelElementA, labelElementB;
  // JL20150122: Note lableCorner[0] and labelVrtx[0] are not used
  int labelCorner[5], labelVrtx[9], labelLast[7], labelThis[7];
  int *NumFcsEm;
  double hx, hy, hz;
  double *x, *y, *z;

  NumNds = 0;  NumEms = 0;  NumFcs = 0;  NumBndryPcs = 0;
  BgnLblNd = 1;  BgnLblEm = 1;  BgnLblFc = 1;  BgnLblBndry = 1;  // by default
  IsBndryNd = 0;  NumNdsNd = 0;  NumEmsNd = 0;  NumFcsNd = 0;
  LblNdNd = 0;  LblNdEm = 0;  LblNdFc = 0;
  IsBndryFc = 0;  LblFcNd = 0;  LblFcEm = 0;  LblFcFc = 0;
  LblEmNd = 0;  LblEmFc = 0;  SignEmFc = 0;
  nd = 0;  OutNmlBndry = 0;
  flag = 0;

  NumNds = (nx+1)*(ny+1)*(nz+1);
  NumFcs = (nx+1)*ny*nz*2 + nx*(ny+1)*nz*2 + nx*ny*(nz+1)*2 + nx*ny*nz*6;
  NumEms = 6*nx*ny*nz;
  NumBndryPcs = 6;

  // BgnLblNd = 1;  BgnLblFc = 1;  BgnLblEm = 1;  BgnLblBndry = 1;

  x = new double[nx+1];
  hx = (xb-xa)/nx;
  x[0] = xa;
  for (i=1; i<=nx-1; ++i)  x[i] = x[i-1] + hx;
  x[nx] = xb;  // Avoiding rounding errors
  
  y = new double[ny+1];
  hy = (yd-yc)/ny;
  y[0] = yc;
  for (j=1; j<=ny-1; ++j)  y[j] = y[j-1] + hy;
  y[ny] = yd;  // Avoiding rounding errors
  
  z = new double[nz+1];
  hz = (zf-ze)/nz;
  z[0] = ze;
  for (k=1; k<=nz-1; ++k)  z[k] = z[k-1] + hz;
  z[nz] = zf;  // Avoiding rounding errors

  // Generating nodal coordinates
  nd = new PtVec3d[NumNds];
  ld = 0;
  for (k=0; k<=nz; ++k) {
    for (j=0; j<=ny; ++j) {
      for (i=0; i<=nx; ++i) {
        nd[ld] = PtVec3d(x[i], y[j], z[k]);
        ld++;
      }
    }
  }
  delete[] x;
  delete[] y;
  delete[] z;
  
  // Setting boundary node flags/marks
  IsBndryNd = new int[NumNds];
  // For z-direction
  for (j=0; j<=ny; ++j) {
    for (i=0; i<=nx; ++i) {
      int ji = j*(nx+1) + i;
      IsBndryNd[ji] = 1;
      ld = nz*(ny+1)*(nx+1) + ji;
      IsBndryNd[ld] = 1;
    }
  }
  // For y-direction
  for (k=0; k<=nz; ++k) {
    for (i=0; i<=nx; ++i) {
      int ki = k*(ny+1)*(nx+1) + i;
      IsBndryNd[ki] = 1;
      ld = ny*(nx+1) + ki;
      IsBndryNd[ld] = 1;
    }
  }
  // For x-direction
  for (k=0; k<=nz; ++k) {
    for (j=0; j<=ny; ++j) {
      int kj = k*(ny+1)*(nx+1) + j*(nx+1);
      IsBndryNd[kj] = 1;
      ld = nx + kj;
      IsBndryNd[ld] = 1;
    }
  }
  
  // Allocation for face info
  LblFcNd = new int[NumFcs][3];
  LblFcEm = new int[NumFcs][2];

  // Generating info for the block faces perpendicular to z-axis
  lc = 0;
  for (k=1; k<=nz+1; ++k) {
    for (j=1; j<=ny; ++j) {
      for (i=1; i<=nx; ++i) {
        //
        n = (((k-1)-1)*ny*nx+(j-1)*nx+i-1)*6;
        for (m=1; m<=6; ++m)  labelLast[m] = m+n;
        if (k==1)  {for (m=1; m<=6; ++m)  labelLast[m] = 0;}
        //
        if (k==nz+1) {
          for (m=1; m<=6; ++m)  labelThis[m] = 0;
        } else {
          n = ((k-1)*ny*nx+(j-1)*nx+i-1)*6;
          for (m=1; m<=6; ++m)  labelThis[m] = m+n;
        }

        labelCorner[1] = glbLblNd(i, j, k, nx, ny);
        labelCorner[2] = glbLblNd(i+1, j, k, nx, ny);
        labelCorner[3] = glbLblNd(i+1, j+1, k, nx, ny);
        labelCorner[4] = glbLblNd(i, j+1, k, nx, ny);

        LblFcNd[lc][0] = labelCorner[1];
        LblFcNd[lc][1] = labelCorner[3];
        LblFcNd[lc][2] = labelCorner[4];

        LblFcEm[lc][0] = labelLast[3];
        LblFcEm[lc][1] = labelThis[1];

        LblFcNd[lc+1][0] = labelCorner[1];
        LblFcNd[lc+1][1] = labelCorner[2];
        LblFcNd[lc+1][2] = labelCorner[3];

        LblFcEm[lc+1][0] = labelLast[6];
        LblFcEm[lc+1][1] = labelThis[4];

        lc += 2;
      }
    }
  }

  // Generating info for the block faces perpendicular to y-axis
  lc = 2*nx*ny*(nz+1);
  for (j=1; j<=ny+1; ++j) {
    for (k=1; k<=nz; ++k) {
      for (i=1; i<=nx; ++i) {
        //
        if (j==1) {
          for (m=1; m<=6; ++m)  labelLast[m] = 0;
        } else {
          n = ((k-1)*ny*nx+((j-1)-1)*nx+i-1)*6;
          for (m=1; m<=6; ++m)  labelLast[m] = m+n;
        }
        //
        if (j==(ny+1)) {
          for (m=1; m<=6; ++m)  labelThis[m] = 0;
        } else {
          n = ((k-1)*ny*nx+(j-1)*nx+i-1)*6;
          for (m=1; m<=6; ++m)  labelThis[m] = m+n;
        }
        
        labelCorner[1] = glbLblNd(i, j, k, nx, ny);
        labelCorner[2] = glbLblNd(i+1, j, k, nx, ny);
        labelCorner[3] = glbLblNd(i+1, j, k+1, nx, ny);
        labelCorner[4] = glbLblNd(i, j, k+1, nx, ny);

        LblFcNd[lc][0] = labelCorner[1];
        LblFcNd[lc][1] = labelCorner[3];
        LblFcNd[lc][2] = labelCorner[4];

        LblFcEm[lc][0] = labelLast[3];
        LblFcEm[lc][1] = labelThis[5];

        LblFcNd[lc+1][0] = labelCorner[1];
        LblFcNd[lc+1][1] = labelCorner[2];
        LblFcNd[lc+1][2] = labelCorner[3];
        
        LblFcEm[lc+1][0] = labelLast[2];
        LblFcEm[lc+1][1] = labelThis[4];

        lc += 2;
      }
    }
  }

  // Generating info for the block faces perpendicualr to x-axis
  lc = 2*nx*ny*(nz+1) + 2*nx*nz*(ny+1);
  for (i=1; i<=nx+1; ++i) {
    for (k=1; k<=nz; ++k) {
      for (j=1; j<=ny; ++j) {
        //
        if (i==1) {
          for (m=1; m<=6; ++m)  labelLast[m] = 0;
        } else {
          n = ((k-1)*ny*nx+(j-1)*nx+(i-1)-1)*6;
          for (m=1; m<=6; ++m)  labelLast[m] = m+n;
        }
        //
        if (i==(nx+1)) {
          for (m=1; m<=6; ++m)  labelThis[m] = 0;
        } else {
          n = ((k-1)*ny*nx+(j-1)*nx+i-1)*6;
          for (m=1; m<=6; ++m)  labelThis[m] = m+n;
        }
        //
        labelCorner[1] = glbLblNd(i, j, k, nx, ny);
        labelCorner[2] = glbLblNd(i, j+1, k, nx, ny);
        labelCorner[3] = glbLblNd(i, j+1, k+1, nx, ny);
        labelCorner[4] = glbLblNd(i, j, k+1, nx, ny);
        //
        LblFcNd[lc][0] = labelCorner[1];
        LblFcNd[lc][1] = labelCorner[2];
        LblFcNd[lc][2] = labelCorner[4];
        //
        LblFcEm[lc][0] = labelLast[4];
        LblFcEm[lc][1] = labelThis[1];
        //
        LblFcNd[lc+1][0] = labelCorner[2];
        LblFcNd[lc+1][1] = labelCorner[3];
        LblFcNd[lc+1][2] = labelCorner[4];
        //
        LblFcEm[lc+1][0] = labelLast[6];
        LblFcEm[lc+1][1] = labelThis[3];
        lc += 2;
      }
    }
  }

  // Setting boundary face flags/marks
  IsBndryFc = new int[NumFcs];
  for (lc=0; lc<NumFcs; ++lc)  IsBndryFc[lc] = 0;
  //
  // For the block faces at the 2 ends of z-axis
  for (j=1; j<=ny; ++j) {
    for (i=1; i<=nx; ++i) {
      // low end ;
      lc = ((j-1)*nx+(i-1))*2;
      IsBndryFc[lc] = 5;
      IsBndryFc[lc+1] = 5;
      // high end
      lc += nz*(ny*nx*2);
      IsBndryFc[lc] = 6;
      IsBndryFc[lc+1] = 6;
    }
  }
  // For the block faces at the 2 ends of y-axis
  for (k=1; k<=nz; ++k) {
    for (i=1; i<=nx; ++i) {
      // low end
      lc = 2*nx*ny*(nz+1) + ((k-1)*nx+(i-1))*2;
      IsBndryFc[lc] = 3;
      IsBndryFc[lc+1] = 3;
      // high end
      lc += ny*(nz*nx*2);
      IsBndryFc[lc] = 4;
      IsBndryFc[lc+1] = 4;
    }
  }
  // For the block faces at the 2 ends of x-axis
  for (k=1; k<=nz; ++k) {
    for (j=1; j<=ny; ++j) {
      // low end
      lc = 2*nx*ny*(nz+1) + 2*nx*nz*(ny+1) + ((k-1)*ny+(j-1))*2;
      IsBndryFc[lc] = 1;
      IsBndryFc[lc+1] = 1;
      // high end
      lc += nx*(nz*ny*2);
      IsBndryFc[lc] = 2;
      IsBndryFc[lc+1] = 2;
    }
  }

  // Generating info for the faces inside the boxes
  lc = 2*nx*ny*(nz+1) + 2*nx*nz*(ny+1) + 2*ny*nz*(nx+1);
  for (k=1; k<=nz; ++k) {
    for (j=1; j<=ny; ++j) {
      for (i=1; i<=nx; ++i) {
        // cmptGlbLblsVrtxBlk(lblVrtx=labelVrtx, i, j, k, nx, ny);
        labelVrtx[1] = (k-1)*(ny+1)*(nx+1) + (j-1)*(nx+1) + i;
        labelVrtx[2] = labelVrtx[1] + 1;
        labelVrtx[3] = labelVrtx[1] + nx+2;
        labelVrtx[4] = labelVrtx[1] + nx+1;
        labelVrtx[5] = labelVrtx[1] + (ny+1)*(nx+1);
        labelVrtx[6] = labelVrtx[5] + 1;
        labelVrtx[7] = labelVrtx[5] + nx+2;
        labelVrtx[8] = labelVrtx[5] + nx+1;

        n = ((k-1)*ny*nx+(j-1)*nx+i-1)*6;
        for (m=1; m<=6; ++m)  labelThis[m] = m+n;

        LblFcNd[lc][0] = labelVrtx[1];
        LblFcNd[lc][1] = labelVrtx[3];
        LblFcNd[lc][2] = labelVrtx[5];

        LblFcEm[lc][0] = labelThis[1];
        LblFcEm[lc][1] = labelThis[5];

        LblFcNd[lc+1][0] = labelVrtx[3];
        LblFcNd[lc+1][1] = labelVrtx[5];
        LblFcNd[lc+1][2] = labelVrtx[7];

        LblFcEm[lc+1][0] = labelThis[2];
        LblFcEm[lc+1][1] = labelThis[6];

        LblFcNd[lc+2][0] = labelVrtx[3];
        LblFcNd[lc+2][1] = labelVrtx[4];
        LblFcNd[lc+2][2] = labelVrtx[5];

        LblFcEm[lc+2][0] = labelThis[1];
        LblFcEm[lc+2][1] = labelThis[2];

        LblFcNd[lc+3][0] = labelVrtx[4];
        LblFcNd[lc+3][1] = labelVrtx[5];
        LblFcNd[lc+3][2] = labelVrtx[7];

        LblFcEm[lc+3][0] = labelThis[2];
        LblFcEm[lc+3][1] = labelThis[3];

        LblFcNd[lc+4][0] = labelVrtx[1];
        LblFcNd[lc+4][1] = labelVrtx[3];
        LblFcNd[lc+4][2] = labelVrtx[6];

        LblFcEm[lc+4][0] = labelThis[4];
        LblFcEm[lc+4][1] = labelThis[5];

        LblFcNd[lc+5][0] = labelVrtx[3];
        LblFcNd[lc+5][1] = labelVrtx[5];
        LblFcNd[lc+5][2] = labelVrtx[6];

        LblFcEm[lc+5][0] = labelThis[5];
        LblFcEm[lc+5][1] = labelThis[6];

        lc += 6;
      }
    }
  }
  
  // JL20150525: CRITICAL: "SWAPPING" for boundary faces
  for (lc=0; lc<NumFcs; ++lc) {
    // if (IsBndryFc[lc]==0)  continue;
    // if ((LblFcEm[lc][0]<=0) && (LblFcEm[lc][1]>0)) {
    if (LblFcEm[lc][0]<=0) {
      LblFcEm[lc][0] = LblFcEm[lc][1];
      LblFcEm[lc][1] = 0;
    }
  }

  // Generating element-vs-(4)nodes info
  LblEmNd = new int[NumEms][4];
  le = 0;
  for (k=1; k<=nz; ++k) {
    for (j=1; j<=ny; ++j) {
      for (i=1; i<=nx; ++i) {
        // cmptGlbLblsVrtxBlk(lblVrtx, i, j, k, nx, ny);
        labelVrtx[1] = (k-1)*(ny+1)*(nx+1) + (j-1)*(nx+1) + i;
        labelVrtx[2] = labelVrtx[1] + 1;
        labelVrtx[3] = labelVrtx[1] + nx+2;
        labelVrtx[4] = labelVrtx[1] + nx+1;
        labelVrtx[5] = labelVrtx[1] + (ny+1)*(nx+1);
        labelVrtx[6] = labelVrtx[5] + 1;
        labelVrtx[7] = labelVrtx[5] + nx+2;
        labelVrtx[8] = labelVrtx[5] + nx+1;
  
        LblEmNd[le][0] = labelVrtx[1];
        LblEmNd[le][1] = labelVrtx[3];
        LblEmNd[le][2] = labelVrtx[4];
        LblEmNd[le][3] = labelVrtx[5];
  
        LblEmNd[le+1][0] = labelVrtx[3];
        LblEmNd[le+1][1] = labelVrtx[4];
        LblEmNd[le+1][2] = labelVrtx[5];
        LblEmNd[le+1][3] = labelVrtx[7];

        LblEmNd[le+2][0] = labelVrtx[4];
        LblEmNd[le+2][1] = labelVrtx[5];
        LblEmNd[le+2][2] = labelVrtx[7];
        LblEmNd[le+2][3] = labelVrtx[8];

        LblEmNd[le+3][0] = labelVrtx[1];
        LblEmNd[le+3][1] = labelVrtx[2];
        LblEmNd[le+3][2] = labelVrtx[3];
        LblEmNd[le+3][3] = labelVrtx[6];

        LblEmNd[le+4][0] = labelVrtx[1];
        LblEmNd[le+4][1] = labelVrtx[3];
        LblEmNd[le+4][2] = labelVrtx[5];
        LblEmNd[le+4][3] = labelVrtx[6];

        LblEmNd[le+5][0] = labelVrtx[3];
        LblEmNd[le+5][1] = labelVrtx[5];
        LblEmNd[le+5][2] = labelVrtx[6];
        LblEmNd[le+5][3] = labelVrtx[7];

        le += 6;
      }
    }
  }

  // Generating info of element-vs-(4)faces based on LblFcEm
  LblEmFc = new int[NumEms][4];
  NumFcsEm = new int[NumEms];
  for (le=0; le<NumEms; ++le)  NumFcsEm[le] = 0;
  for (lc=0; lc<NumFcs; ++lc) {
    labelElementA = LblFcEm[lc][0];
    labelElementB = LblFcEm[lc][1];
    if (labelElementA>0) {
      lea = labelElementA - BgnLblEm;
      LblEmFc[lea][NumFcsEm[lea]] = lc + BgnLblFc;
      NumFcsEm[lea] += 1;
    }
    if (labelElementB>0) {
      leb = labelElementB - BgnLblEm;
      LblEmFc[leb][NumFcsEm[leb]] = lc + BgnLblFc;
      NumFcsEm[leb] += 1;
    }
  }
  delete[] NumFcsEm;

  // Generating info face-vs-(4/7)faces based on LblEmFc
  LblFcFc = new int[NumFcs][7];
  // Initialization
  for (lc=0; lc<NumFcs; ++lc)
    for (l=0; l<7; ++l)
      LblFcFc[lc][l] = 0;
  //
  // Allocating auxiliary ...
  int *Lpos = new int[NumFcs];
  int *Rpos = new int[NumFcs];
  unsigned char *flagFcNbr = new unsigned char[NumFcs];
  // Initialization
  for (lc=0; lc<NumFcs; ++lc) {
    Lpos[lc] = lc + BgnLblFc;
    Rpos[lc] = lc + BgnLblFc;
  }
  //
  for (lc=0; lc<NumFcs; ++lc) {
    //
    for (i=0; i<NumFcs; ++i)  flagFcNbr[i] = 0;
    //
    labelElementA = LblFcEm[lc][0];
    labelElementB = LblFcEm[lc][1];
    lea = labelElementA - BgnLblEm;
    leb = labelElementB - BgnLblEm;
    //
    if (labelElementA>0) {
      for (l=0; l<4; ++l) {
        k = LblEmFc[lea][l];
        flagFcNbr[k-BgnLblFc] = 1;
        if (k<Lpos[lc])  Lpos[lc] = k;
        if (k>Rpos[lc])  Rpos[lc] = k;
      }
    }
    //
    if (labelElementB>0) {
      for (l=0; l<4; ++l) {
        k = LblEmFc[leb][l];
        flagFcNbr[k-BgnLblFc] = 1;
        if (k<Lpos[lc])  Lpos[lc] = k;
        if (k>Rpos[lc])  Rpos[lc] = k;
      }
    }
    j = 0;
    for (i=Lpos[lc]; i<=Rpos[lc]; ++i) {
      if (flagFcNbr[i-BgnLblFc]>0) {
        LblFcFc[lc][j] = i;
        j++;
      }
    }
  }
  // Releasing auxiliary ...
  delete[] Lpos;
  delete[] Rpos;
  delete[] flagFcNbr;

  // Info of outward normal vectors on boundaries
  OutNmlBndry = new PtVec3d[NumBndryPcs];
  OutNmlBndry[0] = PtVec3d(-1,0,0);
  OutNmlBndry[1] = PtVec3d( 1,0,0);
  OutNmlBndry[2] = PtVec3d(0,-1,0);
  OutNmlBndry[3] = PtVec3d(0, 1,0);
  OutNmlBndry[4] = PtVec3d(0,0,-1);
  OutNmlBndry[5] = PtVec3d(0,0, 1);

  // Level-1 dynamic mesh
  flag = 1;
}


// TetraMesh: Destructor

TetraMesh::~TetraMesh()
{
  if (flag==1) {
    for (int ld=0; ld<=NumNds-1; ++ld) {
      delete LblNdNd[ld];
      delete LblNdEm[ld];
    }
  }
  delete[] IsBndryNd, NumNdsNd, LblNdNd, NumEmsNd, LblNdEm;
  delete[] IsBndryFc, LblFcNd, LblFcFc, LblFcEm;
  delete[] LblEmNd, LblEmFc, SignEmFc;
  delete[] nd, OutNmlBndry;
  NumNds = 0;  NumEms = 0;  NumFcs = 0;  NumBndryPcs = 0;
  BgnLblNd = 1;  BgnLblEm = 1;  BgnLblFc = 1;  BgnLblBndry = 1;  // by default
  IsBndryNd = 0;  NumNdsNd = 0;  NumEmsNd = 0;  NumFcsNd = 0;
  LblNdNd = 0;  LblNdEm = 0;  LblNdFc = 0;
  IsBndryFc = 0;  LblFcNd = 0;  LblFcEm = 0;  LblFcFc = 0;
  LblEmNd = 0;  LblEmFc = 0;  SignEmFc = 0;
  nd = 0;  OutNmlBndry = 0;
  flag = 0;
}


void TetraMesh::fillNodeInfo(int numberNodes, double (*crd)[3],
                             int *boundaryNodeMark)
{
  int id, j;

  //
  NumNds = numberNodes;
  nd = new PtVec3d[NumNds];
  IsBndryNd = new int[NumNds];
  
  //
  for (id=0; id<NumNds; ++id)
    nd[id] = PtVec3d(crd[id][0], crd[id][1], crd[id][2]);
  
  //
  for (id=0; id<NumNds; ++id)
    IsBndryNd[id] = boundaryNodeMark[id];
  
  return;
}


void TetraMesh::fillFaceInfo(int numberFaces, int (*faceNode)[3],
                             int (*faceElement)[2], int *boundaryFaceMark)
{
  int ic, j;
  
  //
  NumFcs = numberFaces;
  LblFcNd = new int[NumFcs][3];
  LblFcEm = new int[NumFcs][2];
  IsBndryFc = new int[NumFcs];
  
  //
  for (ic=0; ic<NumFcs; ++ic)
    for (j=0; j<3; ++j)
      LblFcNd[ic][j] = faceNode[ic][j];
  
  //
  for (ic=0; ic<NumFcs; ++ic)
    for (j=0; j<2; ++j)
      LblFcEm[ic][j] = faceElement[ic][j];

  //
  for (ic=0; ic<NumFcs; ++ic)
    IsBndryFc[ic] = boundaryFaceMark[ic];
  
  // JL20150525: CRITICAL: "SWAPPING" for boundary faces
  for (ic=0; ic<NumFcs; ++ic) {
    if (IsBndryFc[ic]==0)  continue;
    if ((LblFcEm[ic][0]<=0) && (LblFcEm[ic][1]>0)) {
      LblFcEm[ic][0] = LblFcEm[ic][1];
      LblFcEm[ic][1] = 0;
    }
  }
  
  return;
}


void TetraMesh::fillElementInfo(int numberElements, int (*elementNode)[4])
{
  int ie, j;
  
  //
  NumEms = numberElements;
  LblEmNd = new int[NumEms][4];
  
  //
  for (ie=0; ie<NumEms; ++ie)
    for (j=0; j<4; ++j)
      LblEmNd[ie][j] = elementNode[ie][j];

  return;
}


// TetraMesh: Enriching to Type I mesh for CG methods
// Generating node-nodes info

void TetraMesh::enrich1()
{
  int i, id, ie, j, jd, k, labele, labelVertex[4];
  unsigned char *flagNdNbr;
  int *Lpos, *Rpos;
  
  // For all nodes: Finding the number of neighboring elements & their labels
  // using LblEmNd
  NumEmsNd = new int[NumNds];
  LblNdEm = new int*[NumNds];
  //
  for (id=0; id<NumNds; ++id)  NumEmsNd[id] = 0;
  //
  for (ie=0; ie<NumEms; ++ie) {
    for (j=0; j<4; ++j) {
      k = LblEmNd[ie][j] - BgnLblNd;
      NumEmsNd[k] = NumEmsNd[k] + 1;
    }
  }
  //
  for (id=0; id<NumNds; ++id)
    LblNdEm[id] = new int[NumEmsNd[id]];
  // A bit unusal but working well
  for (id=0; id<NumNds; ++id)  NumEmsNd[id] = 0;
  for (ie=0; ie<NumEms; ++ie) {
    for (j=0; j<4; ++j) {
      k = LblEmNd[ie][j] - BgnLblNd;
      LblNdEm[k][NumEmsNd[k]] = ie + BgnLblEm;
      NumEmsNd[k] = NumEmsNd[k] + 1;
    }
  }

  // For all nodes: Allocation for auxiliary stuff
  Lpos = new int[NumNds];
  Rpos = new int[NumNds];

  // For all nodes: Calculations of aux. stuff: leftmost-, rightmost- positions
  for (id=0; id<NumNds; ++id) {
    Lpos[id] = id + BgnLblNd;
    Rpos[id] = id + BgnLblNd;
  }
  for (ie=0; ie<NumEms; ++ie) {
    for (i=0; i<4; ++i)
      labelVertex[i] = LblEmNd[ie][i];
    for (i=0; i<4; ++i) {
      id = labelVertex[i] - BgnLblNd;
      for (j=0; j<4; ++j) {
        if (labelVertex[j]<Lpos[id])  Lpos[id] = labelVertex[j];
        if (labelVertex[j]>Rpos[id])  Rpos[id] = labelVertex[j];
      }
    }
  }

  // For all nodes: Finding the number of its neighboring nodes & their labels
  // using LblNdEm
  flagNdNbr = new unsigned char[NumNds];
  NumNdsNd = new int[NumNds];
  LblNdNd = new int*[NumNds];
  //
  for (id=0; id<NumNds; ++id) {
    for (jd=0; jd<NumNds; ++jd)  flagNdNbr[jd] = 0;
    for (ie=0; ie<NumEmsNd[id]; ++ie) {
      labele = LblNdEm[id][ie];
      for (j=0; j<4; ++j) {
        labelVertex[j] = LblEmNd[labele-BgnLblEm][j];
        flagNdNbr[labelVertex[j]-BgnLblNd] = 1;
      }
    }
    //
    NumNdsNd[id] = 0;
    for (jd=Lpos[id]-1; jd<=Rpos[id]-1; ++jd)
      NumNdsNd[id] += flagNdNbr[jd];
    //
    LblNdNd[id] = new int[NumNdsNd[id]];
    k = 0;
    for (jd=Lpos[id]-1; jd<=Rpos[id]-1; ++jd) {
      if (flagNdNbr[jd]==1) {
        LblNdNd[id][k] = jd + BgnLblNd;
        k++;
      }
    }
  }

  // Releasing/Cleanup
  delete[] flagNdNbr, Lpos, Rpos;

  flag = 1;
  return;
}


// TetraMesh: Enriching to Type II mesh for DG methods
// Generating element-elements, face-elements info

void TetraMesh::enrich2()
{
  flag = 2;
  return;
}


// TetraMesh: Enriching to Type III mesh for WG methods
// Generating element-faces, face-elements, face-faces info

void TetraMesh::enrich3()
{
  int i, j, k, l, lc, le, lea, leb;
  int labelElementA, labelElementB;
  int *NumFcsEm;
  int *Lpos;
  int *Rpos;
  int *flagFcNbr;

  // std::cout << "NumEms=" << NumEms << "\n";
  // std::cout << "NumFcs=" << NumFcs << "\n";
  
  // Generating info of element-vs-(4)faces based on LblFcEm
  LblEmFc = new int[NumEms][4];
  NumFcsEm = new int[NumEms];
  for (le=0; le<NumEms; ++le)
    NumFcsEm[le] = 0;
  for (lc=0; lc<NumFcs; ++lc) {
    labelElementA = LblFcEm[lc][0];
    labelElementB = LblFcEm[lc][1];
    if (labelElementA>0) {
      lea = labelElementA - BgnLblEm;
      LblEmFc[lea][NumFcsEm[lea]] = lc + BgnLblFc;
      NumFcsEm[lea] += 1;
    }
    if (labelElementB>0) {
      leb = labelElementB - BgnLblEm;
      LblEmFc[leb][NumFcsEm[leb]] = lc + BgnLblFc;
      NumFcsEm[leb] += 1;
    }
  }
  delete[] NumFcsEm;

  // Generating info face-vs-(4/7)faces based on LblEmFc
  LblFcFc = new int[NumFcs][7];
  // Initialization
  for (lc=0; lc<NumFcs; ++lc) {
    for (l=0; l<7; ++l) {
      LblFcFc[lc][l] = 0;
    }
  }
  //
  // Allocating auxiliary ...
  // std::cout << "NumEms=" << NumEms << "\n";
  // std::cout << "NumFcs=" << NumFcs << "\n";
  Lpos = new int[NumFcs];
  Rpos = new int[NumFcs];
  flagFcNbr = new int[NumFcs];
  // Initialization
  for (lc=0; lc<NumFcs; ++lc) {
    Lpos[lc] = lc + BgnLblFc;
    Rpos[lc] = lc + BgnLblFc;
  }
  //
  for (lc=0; lc<NumFcs; ++lc) {
    //
    for (i=0; i<NumFcs; ++i)  flagFcNbr[i] = 0;
    //
    labelElementA = LblFcEm[lc][0];
    labelElementB = LblFcEm[lc][1];
    lea = labelElementA - BgnLblEm;
    leb = labelElementB - BgnLblEm;
    //
    if (labelElementA>0) {
      for (l=0; l<4; ++l) {
        k = LblEmFc[lea][l];
        flagFcNbr[k-BgnLblFc] = 1;
        if (k<Lpos[lc])  Lpos[lc] = k;
        if (k>Rpos[lc])  Rpos[lc] = k;
      }
    }
    //
    if (labelElementB>0) {
      for (l=0; l<4; ++l) {
        k = LblEmFc[leb][l];
        flagFcNbr[k-BgnLblFc] = 1;
        if (k<Lpos[lc])  Lpos[lc] = k;
        if (k>Rpos[lc])  Rpos[lc] = k;
      }
    }
    j = 0;
    for (i=Lpos[lc]; i<=Rpos[lc]; ++i) {
      if (flagFcNbr[i-BgnLblFc]>0) {
        LblFcFc[lc][j] = i;
        j++;
      }
    }
  }
  // Releasing auxiliary ...
  delete[] Lpos;
  delete[] Rpos;
  delete[] flagFcNbr;

  flag = 3;
  return;
}


// TetraMesh: Enriching Category IV mesh info
// Generating mesh info on node-faces, SignEmFc

void TetraMesh::enrich4()
{
  // Setting up NumFcsNd, LblNdFc
  NumFcsNd = new int[NumNds];
  LblNdFc = new int*[NumNds];
  // 1st-round sweeping: Initialization
  for (int ld=0; ld<=NumNds-1; ld++) {
    NumFcsNd[ld] = 0;
  }
  // 1st-round sweeping: Counting
  for (int ic=0; ic<=NumFcs-1; ic++) {
    for (int j=0; j<=2; ++j) {
      int k = LblFcNd[ic][j] - BgnLblNd;
      NumFcsNd[k] += 1;
    }
  }
  // Allocation
  for (int ld=0; ld<=NumNds-1; ld++) {
    LblNdFc[ld] = new int[NumFcsNd[ld]];
  }
  // 2nd-round sweeping: Resetting
  for (int ld=0; ld<=NumNds-1; ld++) {
    NumFcsNd[ld] = 0;
  }
  // 2nd-round sweeping: Recording
  for (int ic=0; ic<=NumFcs-1; ic++) {
    for (int j=0; j<=2; ++j) {
      int k = LblFcNd[ic][j] - BgnLblNd;
      LblNdFc[k][NumFcsNd[k]] = ic + BgnLblFc;
      NumFcsNd[k] += 1;
    }
  }

  // Auxiliary SignFcEm
  int (*SignFcEm)[2] = new int[NumFcs][2];
  // Default values
  for (int lc=0; lc<=NumFcs-1; ++lc) {
    SignFcEm[lc][0] =  1;
    SignFcEm[lc][1] = -1;
  }

  // JL20180217: Temporarily as placeholder
  // Setting up SignEmFc
  SignEmFc = new int[NumEms][4];
  for (int ie=0; ie<=NumEms-1; ie++) {
    for (int j=0; j<=3; ++j) {
      SignEmFc[ie][j] = 1;
    }
  }
  
  delete[] SignFcEm;
  
  flag = 4;
  return;
}


// TetraMesh: get an element

Tetra TetraMesh::element(int labele) const
{
  int ie = labele - BgnLblEm;
  PtVec3d A = nd[LblEmNd[ie][0]-BgnLblNd];
  PtVec3d B = nd[LblEmNd[ie][1]-BgnLblNd];
  PtVec3d C = nd[LblEmNd[ie][2]-BgnLblNd];
  PtVec3d D = nd[LblEmNd[ie][3]-BgnLblNd];
  return Tetra(A,B,C,D);
}


// TetraMesh: get a face

Tri3d TetraMesh::face(int labelc) const
{
  int ic = labelc - BgnLblFc;
  PtVec3d A = nd[LblFcNd[ic][0]-BgnLblNd];
  PtVec3d B = nd[LblFcNd[ic][1]-BgnLblNd];
  PtVec3d C = nd[LblFcNd[ic][2]-BgnLblNd];
  return Tri3d(A,B,C);
}


// TetraMesh: get neighboring nodes for a given node

void TetraMesh::getNodeNode(int labeld, int *labelNode) const
{
  int id = labeld - BgnLblNd;
  for (int i=0; i<NumNdsNd[id]; ++i)
    labelNode[i] = LblNdNd[id][i];
  return;
}


// TetraMesh: get neighboring elements for a given node

void TetraMesh::getNodeElement(int labeld, int *labelElement) const
{
  int id = labeld - BgnLblNd;
  for (int i=0; i<NumEmsNd[id]; ++i)
    labelElement[i] = LblNdEm[id][i];
  return;
}


// TetraMesh: get neighboring faces for a given node

void TetraMesh::getNodeFace(int labeld, int *labelFace) const
{
  int id = labeld - BgnLblNd;
  for (int i=0; i<NumFcsNd[id]; ++i)
    labelFace[i] = LblNdFc[id][i];
  return;
}


// TetraMesh: get face-vs-(3)nodes (or vertcies) info

void TetraMesh::getFaceNode(int labelc, int labelVertex[3]) const
{
  int ic = labelc - BgnLblFc;
  for (int j=0; j<3; ++j)  labelVertex[j] = LblFcNd[ic][j];
  return;
}


// TetraMesh: get face-vs-(4/7)faces (neighborhood info

void TetraMesh::getFaceFace(int labelc, int labelFaceNbr[7]) const
{
  int ic = labelc - BgnLblFc;
  if (IsBndryFc[ic]>0) {
    for (int j=0; j<4; ++j)  labelFaceNbr[j] = LblFcFc[ic][j];
  }
  if (IsBndryFc[ic]==0) {
    for (int j=0; j<7; ++j)  labelFaceNbr[j] = LblFcFc[ic][j];
  }
  return;
}


// TetraMesh: get face-vs-(1/2)element(s) info

void TetraMesh::getFaceElement(int labelc,
                               int &labelElementA, int &labelElementB) const
{
  int lc = labelc - BgnLblFc;
  labelElementA = LblFcEm[lc][0];
  labelElementB = LblFcEm[lc][1];
  return;
}


// TetraMesh: get element-vs-(4)nodes info

void TetraMesh::getElementNode(int labele, int labelVertex[4]) const
{
  int le = labele - BgnLblEm;
  for (int j=0; j<4; ++j)  labelVertex[j] = LblEmNd[le][j];
  return;
}


// TetraMesh: get element-vs-(4)faces info

void TetraMesh::getElementFace(int labele, int labelFace[4]) const
{
  int le = labele - BgnLblEm;
  for (int j=0; j<4; ++j)  labelFace[j] = LblEmFc[le][j];
  return;
}


// TetraMesh: get element-vs-(4)faces signs

void TetraMesh::getElementFaceSign(int labele, int sign[4]) const
{
  int le = labele - BgnLblEm;
  for (int j=0; j<4; ++j)  sign[j] = SignEmFc[le][j];
  return;
}


// TetraMesh: max mesh size
/*
double TetraMesh::maxSize() const 
{
  int ie, j;
  double h, maxEdge, minEdge;
  PtVec3d vrtx[4];

  h = 0;
  for (ie=0; ie<NumEms; ++ie) {
    for (j=0; j<4; ++j) 
      vrtx[j] = nd[LblEmVrtx[ie][j]-BgnLblNd];
    Tetra Te(vrtx);
    Te.getMaxMinEdges(maxEdge, minEdge);
    if (maxEdge>h)  h=maxEdge;
  }

  return h;
}
*/


// TetraMesh: min mesh size
/*
double TetraMesh::minSize() const 
{
  int ie, j;
  double h, maxEdge, minEdge;
  PtVec3d vrtx[4];

  h = 1E15;
  for (ie=0; ie<NumEms; ++ie) {
    for (j=0; j<4; ++j) 
      vrtx[j] = nd[LblEmVrtx[ie][j]-BgnLblNd];
    Tetra Te(vrtx);
    Te.getMaxMinEdges(maxEdge, minEdge);
    if (minEdge<h)  h=minEdge;
  }

  return h;
}
*/


// TetraMesh: save to a data file

void TetraMesh::save2file(char *filename) const
{
  FILE *fp;
  int ib, ic, id, ie, j;

  if (NULL==(fp=fopen(filename, "w"))) {
    std::cout << "Open data file failed." << "\n";
    exit(-1);
  }

  fprintf(fp, "%7d  %7d  %7d  %7d\n", NumNds, NumFcs, NumEms, NumBndryPcs);
  fprintf(fp, "%2d  %2d  %2d  %2d\n", BgnLblNd, BgnLblFc, BgnLblEm, BgnLblBndry);

  fprintf(fp, "Node info\n");
  for (id=0; id<NumNds; ++id) {
    fprintf(fp, "%7d  ", id+BgnLblNd);
    fprintf(fp, "%2d  ", IsBndryNd[id]);
    fprintf(fp, "%12.6f  ", nd[id].xCrd());
    fprintf(fp, "%12.6f  ", nd[id].yCrd());
    fprintf(fp, "%12.6f\n", nd[id].zCrd());
  }

  fprintf(fp, "Node-vs-nodes info\n");
  for (id=0; id<NumNds; ++id) {
    fprintf(fp, "%7d  ", id+BgnLblNd);
    fprintf(fp, "%3d  ", NumNdsNd[id]);
    for (j=0; j<NumNdsNd[id]; ++j)
      fprintf(fp, "%6d  ", LblNdNd[id][j]);
    fprintf(fp, "\n");
  }

  fprintf(fp, "Face-vs node, element, boundary info\n");
  for (ic=0; ic<NumFcs; ++ic) {
    fprintf(fp, "%7d  ", ic+BgnLblFc);
    fprintf(fp, "%2d  ", IsBndryFc[ic]);
    fprintf(fp, "%7d  ", LblFcNd[ic][0]);
    fprintf(fp, "%7d  ", LblFcNd[ic][1]);
    fprintf(fp, "%7d  ", LblFcNd[ic][2]);
    fprintf(fp, "%7d  ", LblFcEm[ic][0]);
    fprintf(fp, "%7d\n", LblFcEm[ic][1]);
  }

  fprintf(fp, "Face-vs-face(neighbor) info\n");
  for (ic=0; ic<NumFcs; ++ic) {
    fprintf(fp, "%7d  ", ic+BgnLblFc);
    for (j=0; j<7; ++j)
      fprintf(fp, "%6d  ", LblFcFc[ic][j]);
    fprintf(fp, "\n");
  }

  fprintf(fp, "Element-vs-node info\n");
  for (ie=0; ie<NumEms; ++ie)  {
    fprintf(fp, "%7d  ", ie+BgnLblEm);
    for (j=0; j<4; ++j)  fprintf(fp, "%7d  ", LblEmNd[ie][j]);
    fprintf(fp, "\n");
  }

  fprintf(fp, "Element-vs-face info\n");
  for (ie=0; ie<NumEms; ++ie)  {
    fprintf(fp, "%7d  ", ie+BgnLblEm);
    for (j=0; j<4; ++j)  fprintf(fp, "%7d  ", LblEmFc[ie][j]);
    fprintf(fp, "\n");
  }

  fprintf(fp, "Boundary info\n");
  for (ib=0; ib<NumBndryPcs; ++ib) {
    fprintf(fp, "%2d  ", ib+BgnLblBndry);
    fprintf(fp, "%6.3f  ", OutNmlBndry[ib].xCrd());
    fprintf(fp, "%6.3f  ", OutNmlBndry[ib].yCrd());
    fprintf(fp, "%6.3f\n", OutNmlBndry[ib].zCrd());
  }

  fclose(fp);

  return;
}


// TetraMesh: (private) the global index of a node

inline int TetraMesh::glbLblNd(int i, int j, int k, int nx, int ny)
{
  return (k-1)*(nx+1)*(ny+1) + (j-1)*(nx+1) + i;
}

// TetraMesh.cpp
