// Darcy3d_WG_TetraP0P0RT0_AsmGlbMats.cpp
// The Schur-complement approach: Forming 4 global matrices
// James Liu, Graham Harper, ColoState; 2014/07--2017/03

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

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


int Darcy3d_WG_TetraP0P0RT0_AsmGlbMats(BlockDiagMatrix &GlbMatEE,
                                       SparseBlockMatrix &GlbMatEC,
                                       SparseBlockMatrix &GlbMatCE,
                                       SparseBlockMatrix &GlbMatCC,
                                       const TetraMesh &mesh, Mat3 *PermK,
                                       const GaussQuad &GQTe,
                                       const GaussQuad &GQT)
{
  int labelElementA, labelElementB, labelFace[4], labelFaceNeighbor[7], sign[4];
  int NumRows, NumColumns, numRowBands, numColumnBands;
  int *dimRowBand, *dimColumnBand, *numBlocksInRowBand, **columnBandIndex;

  PtVec3d EmCntr, FcCntr, nml;
  Tri3d STF[4], tri;
  Tetra tetra;
  Mat3 MatK;

  // Setup
  int Num0E = mesh.numberElements();
  int Num1C = mesh.numberFaces();
  
  //////////////////////////////////////////////////////////////////////////////
  
  // For GlbMatEE: Element itself interior interaction
  int numBlks = Num0E;
  int *dimBlk = new int[numBlks];
  for (int ie=0; ie<Num0E; ++ie)  dimBlk[ie] = 1;
  GlbMatEE.resize(numBlks, dimBlk);
  delete[] dimBlk;
  
  //////////////////////////////////////////////////////////////////////////////
  
  // For GlbMatEC: Element-Face interaction
  NumRows = Num0E;
  NumColumns = Num1C;
  //
  numRowBands = Num0E;
  numColumnBands = Num1C;
  //
  dimRowBand = new int[numRowBands];
  for (int ie=0; ie<Num0E; ++ie)  dimRowBand[ie] = 1;
  //
  dimColumnBand = new int[numColumnBands];
  for (int jc=0; jc<Num1C; ++jc)  dimColumnBand[jc] = 1;
  //
  numBlocksInRowBand = new int[numRowBands];
  columnBandIndex = new int*[numRowBands];
  //
  // Setting up GlbMatEC: the sparsity structure
  for (int le=0; le<Num0E; ++le) {
    int labele = le + mesh.beginLabelElement();
    mesh.getElementFace(labele, labelFace);
    numBlocksInRowBand[le] = 4;
    columnBandIndex[le] = new int[numBlocksInRowBand[le]];
    for (int j=0; j<4; ++j)
      columnBandIndex[le][j] = labelFace[j];
  }
  //
  // Setting up GlbMatEC: resize
  GlbMatEC.resize(numRowBands, numColumnBands, dimRowBand, dimColumnBand,
                  numBlocksInRowBand, columnBandIndex);
  //
  // Releasing
  for (int i=0; i<numRowBands; ++i)  delete[] columnBandIndex[i];
  delete[] columnBandIndex;
  delete[] numBlocksInRowBand;
  delete[] dimRowBand;
  delete[] dimColumnBand;
  
  //////////////////////////////////////////////////////////////////////////////
  
  // For GlbMatCE: Face-Element(interior) interaction
  NumRows = Num1C;
  NumColumns = Num0E;
  //
  numRowBands = Num1C;
  numColumnBands = Num0E;
  //
  dimRowBand = new int[numRowBands];
  for (int i=0; i<numRowBands; ++i)  dimRowBand[i] = 1;
  //
  dimColumnBand = new int[numColumnBands];
  for (int j=0; j<numColumnBands; ++j)  dimColumnBand[j] = 1;
  //
  numBlocksInRowBand = new int[numRowBands];
  columnBandIndex = new int*[numRowBands];
  //
  // Setting up GlbMatCE: the sparsity structure
  for (int lc=0; lc<Num1C; ++lc) {
    int labelc = lc + mesh.beginLabelFace();
    mesh.getFaceElement(labelc, labelElementA, labelElementB);
    if (mesh.isBoundaryFace(labelc)==0) {  // An interior face
      numBlocksInRowBand[lc] = 2;  // 2 adjacent elements
      columnBandIndex[lc] = new int[numBlocksInRowBand[lc]];
      columnBandIndex[lc][0] = labelElementA;
      columnBandIndex[lc][1] = labelElementB;
    }
    if (mesh.isBoundaryFace(labelc)>0) {  // A boundary face
      numBlocksInRowBand[lc] = 1;  // 1 adjacent element
      columnBandIndex[lc] = new int[numBlocksInRowBand[lc]];
      columnBandIndex[lc][0] = labelElementA;
    }
  }
  // Setting up GlbMatCE: resize
  GlbMatCE.resize(numRowBands, numColumnBands, dimRowBand, dimColumnBand,
                  numBlocksInRowBand, columnBandIndex);
  //
  // Releasing
  for (int i=0; i<numRowBands; ++i)  delete[] columnBandIndex[i];
  delete[] columnBandIndex;
  delete[] numBlocksInRowBand;
  delete[] dimRowBand;
  delete[] dimColumnBand;
  
  //////////////////////////////////////////////////////////////////////////////
  
  // For GlbMatCC: Face-Face interaction
  NumRows = Num1C;
  NumColumns = NumRows;
  //
  numRowBands = Num1C;
  numColumnBands = numRowBands;
  //
  dimRowBand = new int[numRowBands];
  for (int ic=0; ic<Num1C; ++ic)  dimRowBand[ic] = 1;
  //
  dimColumnBand = new int[numColumnBands];
  for (int jc=0; jc<Num1C; ++jc)  dimColumnBand[jc] = 1;
  //
  numBlocksInRowBand = new int[numRowBands];
  columnBandIndex = new int*[numRowBands];
  //
  // Setting up GlbMatCC: the sparsity structure
  for (int lc=0; lc<Num1C; ++lc) {
    int labelc = lc + mesh.beginLabelFace();
    mesh.getFaceFace(labelc, labelFaceNeighbor);
    if (mesh.isBoundaryFace(labelc)==0) {
      numBlocksInRowBand[lc] = 7;  // 7 adjacent faces
      columnBandIndex[lc] = new int[numBlocksInRowBand[lc]];
      for (int j=0; j<7; ++j)
        columnBandIndex[lc][j] = labelFaceNeighbor[j];
    }
    if (mesh.isBoundaryFace(labelc)>0) {
      numBlocksInRowBand[lc] = 4;  // 4 adjacent faces
      columnBandIndex[lc] = new int[numBlocksInRowBand[lc]];
      for (int j=0; j<4; ++j)
        columnBandIndex[lc][j] = labelFaceNeighbor[j];
    }
  }
  //
  // Setting up GlbMatCC: resize
  GlbMatCC.resize(numRowBands, numColumnBands, dimRowBand, dimColumnBand,
                  numBlocksInRowBand, columnBandIndex);
  //
  // Releasing
  for (int i=0; i<numRowBands; ++i)  delete[] columnBandIndex[i];
  delete[] columnBandIndex;
  delete[] numBlocksInRowBand;
  delete[] dimRowBand;
  delete[] dimColumnBand;
  
  //////////////////////////////////////////////////////////////////////////////
  
  // Preparation for assembly
  FullMatrix EGGMK(5,5), EltMat(5,5), SubMat;

  // Assembling the global matrix
  std::cout << "le/1000= ";
  for (int le=0; le<Num0E; ++le) {
    if (le%1000==0)  std::cout << le/1000 << "  " << std::flush;
    int labele = le + mesh.beginLabelElement();
    mesh.getElementFace(labele, labelFace);
    tetra = mesh.element(labele);
    // double diam = tetra.diameter();
    
    for (int j=0; j<4; ++j)  STF[j] = mesh.face(labelFace[j]);
    
    // Determining signs by using centers and normals
    EmCntr = tetra.center();
    for (int 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;
    }
    
    // JL20150502: THIS IS THE WAY TO GO
    MatK = PermK[le];
    
    WG3d_TetraP0P0RT0_EltGradGradMatK(EGGMK, tetra, STF, sign, MatK, GQTe);
    // std::cout << EGGMK << "\n";
    EltMat = EGGMK;

    // For element interior interaction
    SubMat = submatrix(EltMat,1,1,1,1);
    GlbMatEE.setBlock(labele, SubMat);
    
    // For element(interior)-face interaction
    for (int j=1; j<=4; ++j) {
      int jb = labelFace[j-1];
      SubMat = submatrix(EltMat, 1, 1+j, 1, 1);
      GlbMatEC.addBlock(labele, jb, SubMat);
    }
    
    // For face-element(interior) interaction
    for (int i=1; i<=4; ++i) {
      int ib = labelFace[i-1];
      SubMat = submatrix(EltMat, 1+i, 1, 1, 1);
      GlbMatCE.addBlock(ib, labele, SubMat);
    }
    
    // For face-face interaction
    for (int i=1; i<=4; ++i) {
      int ib = labelFace[i-1];
      for (int j=1; j<=4; ++j) {
        int jb = labelFace[j-1];
        SubMat = submatrix(EltMat, 1+i, 1+j, 1, 1);
        GlbMatCC.addBlock(ib, jb, SubMat);
        // GlbMatCC.addBlock(jb, ib, SubMat);
      }
    }
  }
  // std::cout << "\n";

  return(0);  // If successful
}