// Darcy3d_WG_HexaQ0Q0RT0_ModiLinSys.cpp
// The Schur-complement approach
// Modifying the partitioned linear system
// James Liu, Graham Harper, ColoState; 2014/07--2017/03

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

#include "HexaMesh.h"


int Darcy3d_WG_HexaQ0Q0RT0_ModiLinSys(SparseBlockMatrix &GlbMatEF,
                                      SparseBlockMatrix &GlbMatFE,
                                      SparseBlockMatrix &GlbMatFF,
                                      Vector &GlbVecE, Vector &GlbVecF,
                                      Vector &GlbVecD, Vector &GlbVecN,
                                      const HexaMesh &mesh)
{
  int labelElementA, labelElementB, labelFaceNeighbor[11];

  // Auxilliary stuff
  FullMatrix zero1(1,1), eye1(1,1);
  eye1(1,1) = 1;

  // Setup
  int Num0E = mesh.numberElements();
  int Num1F = mesh.numberFaces();
  
  // JL20170311: CHECK CORRECTNESS !!
  // Modifying GlbVecE: RHS "upper"/element-interior part
  GlbVecE = GlbVecE - GlbMatEF * GlbVecD;

  // JL20170311: CHECK CORRECTNESS !!
  // Updating GlbVecF: RHS "lower"/face part
  GlbVecF = GlbVecN - GlbMatFF * GlbVecD;
  for (int labelc=1; labelc<=Num1F; ++labelc) {
    int lc = labelc - 1;
    if (mesh.isBoundaryFace(labelc)>0) {  // Dirichlet face
      GlbVecF[lc] = GlbVecD[lc];
    }
  }
  
  // Modifying GlbMatEF, GlbMatFE simultaneously
  for (int labelc=1; labelc<=Num1F; ++labelc) {
    if (mesh.isBoundaryFace(labelc)>0) {  // Dirichlet face
      mesh.getFaceElement(labelc, labelElementA, labelElementB);
      GlbMatEF.setBlock(labelElementA, labelc, zero1);
      GlbMatFE.setBlock(labelc, labelElementA, zero1);
    }
  }

  // Modifying GlbMatFF
  for (int labelc=1; labelc<=Num1F; ++labelc) {
    if (mesh.isBoundaryFace(labelc)>0) {  // Dirichlet face
      mesh.getFaceFace(labelc, labelFaceNeighbor);
      GlbMatFF.setBlock(labelc, labelc, eye1);
      for (int j=0; j<6; ++j) {
        if (labelFaceNeighbor[j]==labelc)  continue;
        GlbMatFF.setBlock(labelc, labelFaceNeighbor[j], zero1);
        GlbMatFF.setBlock(labelFaceNeighbor[j], labelc, zero1);
      }
    }
  }

  return(0);  // If successful
}