// BrickMesh.cpp
// Class for brick meshes
// James Liu, ColoState; 2007/01--2017/02

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

#include "BrickMesh.h"
#include "cell3d.h"
#include "PtVec3d.h"


// BrickMesh: Default constructor

BrickMesh::BrickMesh() 
{
  NumNds = 0;  NumFcs = 0;  NumEms = 0;  NumBndryPcs = 0;
  BgnLblNd = 0;  BgnLblFc = 0;  BgnLblEm = 0;  BgnLblBndry = 0;
  IsBndryNd = 0;  NumNdsNd = 0;  LblNdNd = 0;  NumEmsNd = 0;  LblNdEm = 0;
  IsBndryFc = 0;  LblFcNd = 0;  LblFcFc = 0;  LblFcEm = 0;
  LblEmNd = 0;  LblEmFc = 0;
  nd = 0;  OutNmlBndry = 0;
  flag = 0;
  X = 0;  Y = 0;  Z = 0;
}


// BrickMesh: A (more meaningful) constructor:
// generating a uniform brick mesh on a brick domain

BrickMesh::BrickMesh(double xa, double xb, int nx,
                     double yc, double yd, int ny,
                     double ze, double zf, int nz)
{
  int i, j, k, l, lc, ld, le;
  // int labelc, labeld, labele;
  int labelElementA, labelElementB, labelVertex[8];
  int *CntEmsFc;
  double hx, hy, hz;

  NumNds = 0;  NumFcs = 0;  NumEms = 0;  NumBndryPcs = 0;
  BgnLblNd = 0;  BgnLblFc = 0;  BgnLblEm = 0;  BgnLblBndry = 0;
  IsBndryNd = 0;  NumNdsNd = 0;  LblNdNd = 0;  NumEmsNd = 0;  LblNdEm = 0;
  IsBndryFc = 0;  LblFcNd = 0;  LblFcFc = 0;  LblFcEm = 0;
  LblEmNd = 0;  LblEmFc = 0;
  nd = 0;  OutNmlBndry = 0;
  flag = 0;

  XA = xa;  XB = xb;  YC = yc;  YD = yd;  ZE = ze;  ZF = zf;
  NX = nx;  NY = ny;  NZ = nz;

  NumNds = (nx+1)*(ny+1)*(nz+1);
  NumFcs = nx*ny*(nz+1) + nz*nx*(ny+1) + ny*nz*(nx+1);
  NumEms = 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

  // Info of node coordinates
  // This also implies an ordering of the nodes
  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++;
      }
    }
  }

  // Info of element-vs-(8)nodes
  LblEmNd = new int[NumEms][8];
  le = 0;
  for (k=1; k<=nz; ++k) {
    for (j=1; j<=ny; ++j) {
      for (i=1; i<=nx; ++i) {
        labelVertex[0] = (k-1)*(ny+1)*(nx+1) + (j-1)*(nx+1) + i;
        labelVertex[1] = labelVertex[0] + 1;
        labelVertex[2] = labelVertex[1] + nx+1;
        labelVertex[3] = labelVertex[0] + nx+1;
        labelVertex[4] = labelVertex[0] + (ny+1)*(nx+1);
        labelVertex[5] = labelVertex[4] + 1;
        labelVertex[6] = labelVertex[5] + nx+1;
        labelVertex[7] = labelVertex[4] + nx+1;
        for (l=0; l<8; ++l)  LblEmNd[le][l] = labelVertex[l];
        le++;
      }
    }
  }
  
  // Info of face-vs-(4)nodes
  // LblFcNd = 0;
  LblFcNd = new int[NumFcs][4];

  // x-faces (faces perpendicular to x-axis)
  lc = 0;
  for (i=0; i<=nx; ++i) {
    for (k=1; k<=nz; ++k) {
      for (j=1; j<=ny; ++j) {
        labelVertex[0] = (k-1)*(ny+1)*(nx+1) + (j-1)*(nx+1) + i+1;
        labelVertex[1] = labelVertex[0] + (nx+1);
        labelVertex[2] = labelVertex[1] + (nx+1)*(ny+1);
        labelVertex[3] = labelVertex[0] + (nx+1)*(ny+1);
        for (l=0; l<4; ++l)  LblFcNd[lc][l] = labelVertex[l];
        lc++;
      }
    }
  }

  // y-faces (faces perpendicular to y-axis)
  lc = (nx+1)*ny*nz;
  for (j=0; j<=ny; ++j) {
    for (k=1; k<=nz; ++k) {
      for (i=1; i<=nx; ++i) {
        labelVertex[0] = (k-1)*(ny+1)*(nx+1) + j*(nx+1) + i;
        labelVertex[1] = labelVertex[0] + 1;
        labelVertex[2] = labelVertex[1] + (nx+1)*(ny+1);
        labelVertex[3] = labelVertex[0] + (nx+1)*(ny+1);
        for (l=0; l<4; ++l)  LblFcNd[lc][l] = labelVertex[l];
        lc++;
      }
    }
  }

  // z-faces (faces perpendicular to z-axis)
  lc = (nx+1)*ny*nz + nx*(ny+1)*nz;
  for (k=0; k<=nz; ++k) {
    for (j=1; j<=ny; ++j) {
      for (i=1; i<=nx; ++i) {
        labelVertex[0] = k*(ny+1)*(nx+1) + (j-1)*(nx+1) + i;
        labelVertex[1] = labelVertex[0] + 1;
        labelVertex[2] = labelVertex[1] + (nx+1);
        labelVertex[3] = labelVertex[0] + (nx+1);
        for (l=0; l<4; ++l)  LblFcNd[lc][l] = labelVertex[l];
        lc++;
      }
    }
  }

  // Info on which faces are on the boundaries
  // IsBndryFc = 0;
  IsBndryFc = new int[NumFcs];
  for (lc=0; lc<NumFcs; ++lc)  IsBndryFc[lc] = 0;
  // x-faces (faces perpendicular to x-axis)
  // i = 0;
  lc = 0;
  for (k=1; k<=nz; ++k) {
    for (j=1; j<=ny; ++j) {
      IsBndryFc[lc] = 1;
      lc++;
    }
  }
  // i = nx;
  lc = nx*ny*nz;
  for (k=1; k<=nz; ++k) {
    for (j=1; j<=ny; ++j) {
      IsBndryFc[lc] = 2;
      lc++;
    }
  }
  // y-faces (faces perpendicular to y-axis)
  // j = 0;
  lc = (nx+1)*ny*nz;
  for (k=1; k<=nz; ++k) {
    for (i=1; i<=nx; ++i) {
      IsBndryFc[lc] = 3;
      lc++;
    }
  }
  // j = ny;
  lc = (nx+1)*ny*nz + nx*ny*nz;
  for (k=1; k<=nz; ++k) {
    for (i=1; i<=nx; ++i) {
      IsBndryFc[lc] = 4;
      lc++;
    }
  }
  // z-faces (faces perpendicular to z-axis)
  // k = 0;
  lc = (nx+1)*ny*nz + nx*(ny+1)*nz;
  for (j=1; j<=ny; ++j) {
    for (i=1; i<=nx; ++i) {
      IsBndryFc[lc] = 5;
      lc++;
    }
  }
  // k = nz;
  lc = (nx+1)*ny*nz + nx*(ny+1)*nz + nz*nx*ny;
  for (j=1; j<=ny; ++j) {
    for (i=1; i<=nx; ++i) {
      IsBndryFc[lc] = 6;
      lc++;
    }
  }
  
  // Info of element-vs-(6)faces
  // LblEmFc = 0;
  LblEmFc = new int[NumEms][6];
  le = 0;
  for (k=1; k<=nz; ++k) {
    for (j=1; j<=ny; ++j) {
      for (i=1; i<=nx; ++i) {
        // x-faces: left, right
        LblEmFc[le][0] = (i-1)*ny*nz + (k-1)*ny + j;
        LblEmFc[le][1] =     i*ny*nz + (k-1)*ny + j;
        // y-faces: back, front
        LblEmFc[le][2] = (j-1)*nx*nz + (k-1)*nx + i + (nx+1)*ny*nz;
        LblEmFc[le][3] =     j*nx*nz + (k-1)*nx + i + (nx+1)*ny*nz;
        // z-faces: bottom, top
        LblEmFc[le][4] = (k-1)*nx*ny + (j-1)*nx + i + (nx+1)*ny*nz + nx*(ny+1)*nz;
        LblEmFc[le][5] =     k*nx*ny + (j-1)*nx + i + (nx+1)*ny*nz + nx*(ny+1)*nz;
        le++;
      }
    }
  }

  // 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);

  // Info of face-vs-(2)elements based on LblEmFc
  // LblFcEm = 0;
  LblFcEm = new int[NumFcs][2];
  CntEmsFc = new int[NumFcs];
  for (lc=0; lc<NumFcs; ++lc) {
    LblFcEm[lc][0] = 0;
    LblFcEm[lc][1] = 0;
    CntEmsFc[lc] = 0;
  }
  for (le=0; le<NumEms; ++le) {
    for (l=0; l<6; ++l) {
      lc = LblEmFc[le][l] - BgnLblFc;
      LblFcEm[lc][CntEmsFc[lc]] = le + BgnLblEm;
      CntEmsFc[lc] = CntEmsFc[lc] + 1;
    }
  }
  delete[] CntEmsFc;
 
  // Info of face-vs-(6/11)faces based on LblEmFc
  // LblFcFc = 0;
  LblFcFc = new int[NumFcs][11];
  // Setting to zero
  for (lc=0; lc<NumFcs; ++lc)
    for (l=0; l<11; ++l)
      LblFcFc[lc][l] = 0;
  for (lc=0; lc<NumFcs; ++lc) {
    labelElementA = LblFcEm[lc][0];
    labelElementB = LblFcEm[lc][1];
    for (l=0; l<6; ++l) {
      LblFcFc[lc][l] = LblEmFc[labelElementA-BgnLblEm][l];
    }
    if (labelElementB>0) {
      j = 0;
      for (l=0; l<6; ++l) {
        k = LblEmFc[labelElementB-BgnLblEm][l];
        if ((k-1)==lc)  continue;
        LblFcFc[lc][6+j] = k;
        j++;
      }
    }
  }

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


// BrickMesh: Constructor based on reading a mesh data file
/**********************************************************************
BrickMesh::BrickMesh(char *filename)
{
  FILE *fp;
  int ib, ic, ic1, ic2, id, ie, ie1, ie2, j;
  int k, k1, k2, label, le1, le2;
  double x, y, z;

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

  fscanf(fp, "%7d", &NumNds);
  fscanf(fp, "%7d", &NumFcs);
  fscanf(fp, "%7d", &NumEms);
  fscanf(fp, "%7d", &NumBndryPcs);
  fscanf(fp, "%2d", &BgnLblNd);
  fscanf(fp, "%2d", &BgnLblFc);
  fscanf(fp, "%2d", &BgnLblEm);
  fscanf(fp, "%2d", &BgnLblBndry);
  flag = 0;

  nd = new PtVec3d[NumNds];
  for (id=0; id<NumNds; ++id) {
    fscanf(fp, "%7d", &label);
    fscanf(fp, "%lf  %lf  %lf", &x, &y, &z);
    fscanf(fp, "%7d", &label);
    nd[id] = PtVec3d(x, y, z);
  }

  LblFcVrtx = new int [NumFcs][3];
  LblFcEm = new int [NumFcs][2];
  IsBndryFc = new int[NumFcs];
  for (ic=0; ic<NumFcs; ++ic) {
    fscanf(fp, "%7d", &label);
    fscanf(fp, "%7d", &LblFcVrtx[ic][0]);
    fscanf(fp, "%7d", &LblFcVrtx[ic][1]);
    fscanf(fp, "%7d", &LblFcVrtx[ic][2]);
    fscanf(fp, "%7d", &LblFcEm[ic][0]);
    fscanf(fp, "%7d", &LblFcEm[ic][1]);
    fscanf(fp, "%7d", &IsBndryFc[ic]);
  }

  LblEmVrtx = new int [NumEms][8];
  for (ie=0; ie<NumEms; ++ie) {
    fscanf(fp, "%7d", &label);
    fscanf(fp, "%7d", &LblEmVrtx[ie][0]);
    fscanf(fp, "%7d", &LblEmVrtx[ie][1]);
    fscanf(fp, "%7d", &LblEmVrtx[ie][2]);
    fscanf(fp, "%7d", &LblEmVrtx[ie][3]);
    fscanf(fp, "%7d", &LblEmVrtx[ie][4]);
    fscanf(fp, "%7d", &LblEmVrtx[ie][5]);
    fscanf(fp, "%7d", &LblEmVrtx[ie][6]);
    fscanf(fp, "%7d", &LblEmVrtx[ie][7]);
  }

  OutNmlBndry = new PtVec3d[NumBndryPcs];
  for (ib=0; ib<NumBndryPcs; ++ib) {
    fscanf(fp, "%2d", &label);
    fscanf(fp, "%lf  %lf  %lf", &x, &y, &z);
    OutNmlBndry[ib] = PtVec3d(x, y, z);
  }

  int (*LblEmNbr)[4] = new int [NumEms][4];
  // to be further used for generating face neighborhood info
  for (ie=0; ie<NumEms; ++ie) {
    fscanf(fp, "%7d", &label);
    fscanf(fp, "%7d", &LblEmNbr[ie][0]);
    fscanf(fp, "%7d", &LblEmNbr[ie][1]);
    fscanf(fp, "%7d", &LblEmNbr[ie][2]);
    fscanf(fp, "%7d", &LblEmNbr[ie][3]);
  }

  fclose(fp);

  // Generating element-vs-face info
  LblEmFc = new int [NumEms][4];
  int *NumFcsEm = new int[NumEms];
  for (ie=0; ie<NumEms; ++ie)  NumFcsEm[ie] = 0;
  for (ic=0; ic<NumFcs; ++ic) {
    le1 = LblFcEm[ic][0];
    le2 = LblFcEm[ic][1]; 
    if (le1>0) {
      ie1 = le1 - BgnLblEm;
      LblEmFc[ie1][NumFcsEm[ie1]] = ic + BgnLblFc;
      NumFcsEm[ie1] += 1;
    }
    if (le2>0) {
      ie2 = le2 - BgnLblEm;    
      LblEmFc[ie2][NumFcsEm[ie2]] = ic + BgnLblFc;
      NumFcsEm[ie2] += 1;
    }
  }

  // Generating face-vs-face neighborhood info
  LblFcNbr = new int [NumFcs][11];
  int *FlagAsFaceNbr = new int[NumFcs];
  for (ic=0; ic<NumFcs; ++ic) {
    for (j=0; j<25; ++j)  LblFcNbr[ic][j] = BgnLblFc - 1;
    for (j=0; j<NumFcs; ++j)  FlagAsFaceNbr[j] = 0;
    FlagAsFaceNbr[ic] = 1;
    le1 = LblFcEm[ic][0];
    le2 = LblFcEm[ic][1]; 
    if (le1>=BgnLblEm) { 
      ie1 = le1 - BgnLblEm;
      for (k=0; k<=3; ++k) {
        ic1 = LblEmFc[ie1][k] - BgnLblFc;
        FlagAsFaceNbr[ic1] = 1;
      }
      for (j=0; j<=3; ++j) { 
        k1 = LblEmNbr[ie1][j]-BgnLblEm;
        if (k1>=0) {
          for (k=0; k<=3; ++k) {
            ic1 = LblEmFc[k1][k] - BgnLblFc;
            FlagAsFaceNbr[ic1] = 1;
          }
        }
      }
    }
    if (le2>=BgnLblEm) { 
      ie2 = le2 - BgnLblEm;
      for (k=0; k<=3; ++k) {
        ic2 = LblEmFc[ie2][k] - BgnLblFc;
        FlagAsFaceNbr[ic2] = 1;
      }
      for (j=0; j<=3; ++j) { 
        k2 = LblEmNbr[ie2][j]-BgnLblEm;
        if (k2>=0) {
          for (k=0; k<=3; ++k) {
            ic2 = LblEmFc[k2][k] - BgnLblFc;
            FlagAsFaceNbr[ic2] = 1;
          }
        }
      }
    }
    k = 0;
    for (j=0; j<NumFcs; ++j) {
      if (FlagAsFaceNbr[j]==1) {
        LblFcNbr[ic][k] = j + BgnLblFc;
        k++;
      }
    }
  }
  delete[] FlagAsFaceNbr, NumFcsEm;
  delete[] LblEmNbr;
  // std::cout << "Reading/generating real mesh finished\n";
}
**********************************************************************/


// BrickMesh: Destructor

BrickMesh::~BrickMesh() 
{
  delete[] IsBndryNd, NumNdsNd, LblNdNd, NumEmsNd, LblNdEm;
  delete[] IsBndryFc, LblFcNd, LblFcFc, LblFcEm;
  delete[] LblEmNd, LblEmFc;
  delete[] nd, OutNmlBndry;
  NumNds = 0;  NumFcs = 0;  NumEms = 0;  NumBndryPcs = 0;
  BgnLblNd = 0;  BgnLblFc = 0;  BgnLblEm = 0;  BgnLblBndry = 0;
  IsBndryNd = 0;  NumNdsNd = 0;  LblNdNd = 0;  NumEmsNd = 0;  LblNdEm = 0;
  IsBndryFc = 0;  LblFcNd = 0;  LblFcFc = 0;  LblFcEm = 0;
  LblEmNd = 0;  LblEmFc = 0;
  nd = 0;  OutNmlBndry = 0;
  flag = 0;
  X = 0;  Y = 0;  Z = 0;
}


// BrickMesh: Enriching to TYpe-I for CG methods

void BrickMesh::enrich1()
{
  unsigned char *flagNdNbr;
  int i, id, ie, ii, j, jj, k, labele, labelVertex[8];
  int *Lpos, *Rpos;
  
  // For each node: Finding the number of neighboring elements & their labels
  NumEmsNd = new int[NumNds];
  LblNdEm = new int*[NumNds];
  //
  for (id=0; id<NumNds; ++id)  NumEmsNd[id] = 0;
  //
  for (labele=1; labele<=NumEms; ++labele) {
    for (j=0; j<8; ++j) {
      k = LblEmNd[labele-BgnLblEm][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 (labele=1; labele<=NumEms; ++labele) {
    for (j=0; j<8; ++j) {
      k = LblEmNd[labele-BgnLblEm][j] - BgnLblNd;
      LblNdEm[k][NumEmsNd[k]] = labele;
      NumEmsNd[k] = NumEmsNd[k] + 1;
    }
  }
  
  // For each node: Auxiliary stuff
  flagNdNbr = new unsigned char[NumNds];
  Lpos = new int[NumNds];
  Rpos = new int[NumNds];
  //
  for (id=0; id<NumNds; ++id) {
    Lpos[id] = id + BgnLblNd;
    Rpos[id] = id + BgnLblNd;
  }
  for (labele=1; labele<=NumEms; ++labele) {
    for (k=0; k<8; ++k) labelVertex[k] = LblEmNd[labele-BgnLblEm][k];
    for (i=0; i<8; ++i) {
      ii = labelVertex[i] - BgnLblNd;
      for (j=0; j<8; ++j) {
        if (labelVertex[j]<Lpos[ii])  Lpos[ii] = labelVertex[j];
        if (labelVertex[j]>Rpos[ii])  Rpos[ii] = labelVertex[j];
      }
    }
  }
  
  // For each node: Finding the number of its neighboring nodes and their labels
  NumNdsNd = new int[NumNds];
  LblNdNd = new int*[NumNds];
  for (ii=0; ii<NumNds; ++ii) {
    for (id=0; id<NumNds; ++id)  flagNdNbr[id] = 0;
    for (ie=0; ie<NumEmsNd[ii]; ++ie) {
      labele = LblNdEm[ii][ie];
      for (k=0; k<8; ++k) {
        labelVertex[k] = LblEmNd[labele-BgnLblEm][k];
        flagNdNbr[labelVertex[k]-BgnLblNd] = 1;
      }
    }
    NumNdsNd[ii] = 0;
    for (id=Lpos[ii]-1; id<=Rpos[ii]-1; ++id)  NumNdsNd[ii] += flagNdNbr[id];
    // Now for the labels of the neighboring nodes
    LblNdNd[ii] = new int[NumNdsNd[ii]];
    k = 0;
    for (id=Lpos[ii]-1; id<=Rpos[ii]-1; ++id) {
      if (flagNdNbr[id]==1) {
        LblNdNd[ii][k] = id + BgnLblNd;
        k++;
      }
    }
  }
  
  // Cleanup
  delete[] flagNdNbr, Lpos, Rpos;

  return;
}


// BrickMesh: Get an element

Brick BrickMesh::element(int labele) const
{
  int ie = labele - BgnLblEm;
  PtVec3d P0 = nd[LblEmNd[ie][0]-BgnLblNd];
  PtVec3d P1 = nd[LblEmNd[ie][6]-BgnLblNd];
  return Brick(P0,P1);
}


// BrickMesh: Get a face

Rect3d BrickMesh::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];
  PtVec3d D = nd[LblFcNd[ic][3]-BgnLblNd];
  return Rect3d(A,B,C,D);
}


// BrickMesh: Get node-vs-nodes info

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


// BrickMesh: Get face-vs-(4)nodes info

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


// BrickMesh: Get face-vs-(11)faces (neighborhood) info

void BrickMesh::getFaceFace(int labelc, int labelFace[11]) const
{
  int ic = labelc - BgnLblFc;
  for (int j=0; j<11; ++j)  labelFace[j] = LblFcFc[ic][j];
  return;
}


// BrickMesh: Get face-vs-(2)elements info

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


// BrickMesh: Get element-vs-(8)nodes info

void BrickMesh::getElementNode(int labele, int labelNode[8]) const
{
  int ie = labele - BgnLblEm;
  for (int j=0; j<8; ++j)  labelNode[j] = LblEmNd[ie][j];
  return;
}


// BrickMesh: Get element-vs-(6)faces info

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


// BrickMesh: max mesh size
/*
double BrickMesh::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) {
    }
  }

  return h;
}
*/


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

  h = 1E15;
  for (ie=0; ie<NumEms; ++ie) {
    for (j=0; j<8; ++j) {
    }
  }

  return h;
}
*/


// BrickMesh: Save to a data file

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

  if (NULL==(fp=fopen(filename, "w"))) {
    puts("Open data file failed.  Exit.\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, "%12.6f  ", nd[id].xCrd());
    fprintf(fp, "%12.6f  ", nd[id].yCrd());
    fprintf(fp, "%12.6f\n", nd[id].zCrd());
  }

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

  fprintf(fp, "Face-vs-face(neighbor) info\n");
  for (ic=0; ic<NumFcs; ++ic) {
    fprintf(fp, "%7d  ", ic+BgnLblFc);
    for (j=0; j<11; ++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<8; ++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<6; ++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;
}

// BrickMesh.cpp