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

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

#include "TetraMesh.h"


int Darcy3d_WG_TetraP0P0RT0_ModiLinSys1(SparseMatrix &GlbMat, Vector &GlbRHS,
                                        const Vector &GlbVecSource,
                                        const Vector &GlbVecDirichlet,
                                        const Vector &GlbVecNeumann,
                                        const TetraMesh &mesh)
{
  int labelElementA, labelElementB, labelFaceNeighbor[7];

  // Setup
  int Num0E = mesh.numberElements();
  int Num1C = mesh.numberFaces();
  int DOFs = Num0E + Num1C;
  GlbRHS.resize(DOFs);

  // Incorporating the Neumann boundary condition into the disc.lin.sys.
  GlbRHS = GlbVecSource - GlbVecNeumann;
  
  // JL20170131: This approach is not recommended hence not implemented here
  // Approach "harsh": Easy
  // Adjusting the disc.lin.sys. according to the Dirichlet bndry.cond.
  
  // Approach "gentle": Maintaining symmetry
  // Adjusting the disc.lin.sys. according to the Dirichlet bndry.cond.
  GlbRHS = GlbRHS - GlbMat * GlbVecDirichlet;
  for (int labelc=1; labelc<=Num1C; ++labelc) {
    if (mesh.isBoundaryFace(labelc)>0) {
      mesh.getFaceElement(labelc, labelElementA, labelElementB);
      mesh.getFaceFace(labelc, labelFaceNeighbor);
      GlbMat.unitfyRow(Num0E+labelc);
      GlbRHS(Num0E+labelc) = GlbVecDirichlet(Num0E+labelc);
      GlbMat.setEntry(labelElementA, Num0E+labelc, 0);
      for (int j=0; j<4; ++j) {  // Actually only 3 neighboring faces
        if (labelFaceNeighbor[j]==labelc)  continue;
        GlbMat.setEntry(Num0E+labelFaceNeighbor[j], Num0E+labelc, 0);
      }
    }
  }

  return(0);  // If successful
}
