// Darcy3d_WG_TetraP0P0RT0_AsmGlbMat1.cpp
// Darcy: Weak Galerkin (P0,P0;RT0) on a tetrahedral mesh
// The single-matrix approach
// James Liu, ColoState; 2014/07--2017/01

#include "matrix.h"
#include "vector.h"

#include "cell3d.h"
#include "GaussQuad.h"
#include "mat3.h"
#include "PtVec3d.h"
#include "TetraMesh.h"
#include "WG3d_Tetra.h"

int Darcy3d_WG_TetraP0P0RT0_AsmGlbMat1(SparseMatrix &GlbMat,
                                       const TetraMesh &mesh, Mat3 *PermK,
                                       const GaussQuad &GQTe)
{
  int i, ii, j, jj, lc, le;
  int labelc, labele, labelElementA, labelElementB;
  int labelFace[4], labelFaceNeighbor[7], sign[4];
  int DOFs, Num0E, Num1C, NumCols, NumRows;
  int *NumEntrsInRow, **ColIdx;
  PtVec3d EmCntr, FcCntr, nml;
  Tri3d STF[4], tri;
  Tetra tetra;
  Mat3 MatK;
  FullMatrix EGGMK(5,5);

  // Setup
  Num0E = mesh.numberElements();
  Num1C = mesh.numberFaces();
  DOFs = Num0E + Num1C;
  NumRows = DOFs;
  NumCols = NumRows;
  NumEntrsInRow = new int[NumRows];
  ColIdx = new int*[NumRows];

  // Setting up the global matrix: element--(1)element+(4)faces interaction
  for (le=0; le<Num0E; ++le) {
    labele = le + mesh.beginLabelElement();
    mesh.getElementFace(labele, labelFace);
    NumEntrsInRow[le] = 5;
    ColIdx[le] = new int[NumEntrsInRow[le]];
    ColIdx[le][0] = labele;
    for (j=0; j<4; ++j)  ColIdx[le][j+1] = Num0E + labelFace[j];
  }

  // Setting up the global matrix: face--(1/2)elements+(4/7)faces interaction
  for (lc=0; lc<Num1C; ++lc) {
    labelc = lc + mesh.beginLabelFace();
    mesh.getFaceElement(labelc, labelElementA, labelElementB);
    mesh.getFaceFace(labelc, labelFaceNeighbor);
    if (mesh.isBoundaryFace(labelc)==0) {
      NumEntrsInRow[Num0E+lc] = 9;  // 2+7
      ColIdx[Num0E+lc] = new int[NumEntrsInRow[Num0E+lc]];
      ColIdx[Num0E+lc][0] = labelElementA;
      ColIdx[Num0E+lc][1] = labelElementB;
      for (j=0; j<7; ++j)
        ColIdx[Num0E+lc][j+2] = Num0E + labelFaceNeighbor[j];
    }
    if (mesh.isBoundaryFace(labelc)>0) {
      NumEntrsInRow[Num0E+lc] = 5;  // 1+4
      ColIdx[Num0E+lc] = new int[NumEntrsInRow[Num0E+lc]];
      ColIdx[Num0E+lc][0] = labelElementA;
      for (j=0; j<4; ++j)
        ColIdx[Num0E+lc][j+1] = Num0E + labelFaceNeighbor[j];
    }
  }

  // Setting up GlbMat: resize & release
  GlbMat.resize(NumRows,NumCols,NumEntrsInRow,ColIdx);
  for (ii=0; ii<NumRows; ++ii)  delete[] ColIdx[ii];
  delete[] ColIdx;
  delete[] NumEntrsInRow;
  
  // Assembling the global matrix
  for (labele=1; labele<=mesh.numberElements(); ++labele) {
    le = labele - mesh.beginLabelElement();
    mesh.getElementFace(labele, labelFace);
    tetra = mesh.element(labele);
    MatK = PermK[le];

    // std::cout << "labele= " << labele << "\n";
    // std::cout << MatK << "\n";

    for (j=0; j<4; ++j)  STF[j] = mesh.face(labelFace[j]);

    // Determining signs by using centers and normals
    EmCntr = tetra.center();
    for (j=0; j<4; ++j) {
      sign[j] = 1;
      tri = STF[j];
      FcCntr = tri.center();
      nml = tri.normal();
      if (dotProduct(nml,FcCntr-EmCntr)<0)  sign[j] = -1;
    }

    WG3d_TetraP0P0RT0_EltGradGradMatK(EGGMK, tetra, STF, sign, MatK, GQTe);
    // std::cout << EGGMK << "\n";

    //
    i = 1;  ii = labele;
    j = 1;  jj = labele;
    GlbMat.addEntry(ii, jj, EGGMK(1,1));
    for (j=1; j<=4; ++j) {
      jj = Num0E + labelFace[j-1];
      GlbMat.addEntry(ii, jj, EGGMK(1,j+1));
    }
    //
    for (i=1; i<=4; ++i) {
      ii = Num0E + labelFace[i-1];
      j = 1;  jj = labele;
      GlbMat.addEntry(ii, jj, EGGMK(i+1,1));
      for (j=1; j<=4; ++j) {
        jj = Num0E + labelFace[j-1];
        GlbMat.addEntry(ii, jj, EGGMK(i+1,j+1));
      }
    }
  }

  return(0);  // if successful
}