// LinElas3d_WG_HexaQ03Q03RT03Q0_ModiLinSys1.cpp
// Solving linear elasticity in 3d by WG(Q0^3,Q0^3;RT[0]^3,Q0) on a hexa.mesh
// The single-matrix approach
// James Liu, Graham Harper, ColoState; 2014/07--2018/01

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

#include "HexaMesh.h"


int LinElas3d_WG_HexaQ03Q03RT03Q0_ModiLinSys1(
  SparseBlockMatrix &GlbMat, Vector &GlbRHS,
  const Vector &ABC, const HexaMesh &mesh, const Vector &flag,
  const Vector &GlbVecSource, const Vector &GlbVecDirichlet,
  const Vector &GlbVecNeumann)
{
  int labelElementA, labelElementB, labelFaceNeighbor[11];
  
  // Setup
  int Num0E = mesh.numberElements();
  int Num1C = mesh.numberFaces();
  int DOFs = 3*Num0E + 3*Num1C;
  GlbRHS.resize(DOFs);

  // Incorporating the Neumann boundary condition into the global lin.sys.
  GlbRHS = GlbVecSource - GlbVecNeumann;
  // GlbRHS.save2file("GlbRHS.dat");

  // Auxiliary
  FullMatrix eye3(3,3), zero33(3,3);
  eye3(1,1) = 1;
  eye3(2,2) = 1;
  eye3(3,3) = 1;

  // Approach "gentle": Maintaining symmetry
  // Adjusting the glb.disc.lin.sys. according to the Dirichlet bndry.cond.
  GlbRHS = GlbRHS - GlbMat * GlbVecDirichlet;
  for (int labelc=1; labelc<=Num1C; ++labelc) {
    int lc = labelc - 1;
    if (ABC(mesh.isBoundaryFace(labelc))==1) {  // For Dirichlet
      mesh.getFaceElement(labelc, labelElementA, labelElementB);
      mesh.getFaceFace(labelc,labelFaceNeighbor);
      GlbRHS(3*Num0E+3*lc+1) = GlbVecDirichlet(3*Num0E+3*lc+1);
      GlbRHS(3*Num0E+3*lc+2) = GlbVecDirichlet(3*Num0E+3*lc+2);
      GlbRHS(3*Num0E+3*lc+3) = GlbVecDirichlet(3*Num0E+3*lc+3);
      GlbMat.unitifyRow(Num0E+labelc,1);
      GlbMat.unitifyRow(Num0E+labelc,2);
      GlbMat.unitifyRow(Num0E+labelc,3);
      GlbMat.setBlock(labelElementA, Num0E+labelc, zero33);
      for (int j=0; j<6; ++j) {
        if (labelFaceNeighbor[j]==labelc)  continue;
        GlbMat.setBlock(Num0E+labelFaceNeighbor[j], Num0E+labelc, zero33);
      }
    }
  }
  // Adjusting the glb.disc.lin.sys. according to bndry.cond. sym. in x
  // GlbRHS: No change
  for (int labelc=1; labelc<=Num1C; ++labelc) {
    int lc = labelc - 1;
    if (ABC(mesh.isBoundaryFace(labelc))==4) {  // For symmetry in x-direction
      /*
      // JL20180126: ACTUALLY NOT NEEDED
      mesh.getFaceElement(labelc, labelElementA, labelElementB);
      mesh.getFaceFace(labelc,labelFaceNeighbor);
      */
      GlbRHS(3*Num0E+3*lc+1) = 0;
      GlbMat.unitifyRow(Num0E+labelc,1);
      GlbMat.unitifyColumn(Num0E+labelc,1);
    }
  }
  // Adjusting the glb.disc.lin.sys. according to bndry.cond. sym. in y
  // GlbRHS: No change
  for (int labelc=1; labelc<=Num1C; ++labelc) {
    int lc = labelc - 1;
    if (ABC(mesh.isBoundaryFace(labelc))==5) {  // For symmetry in y-direction
      /*
      // JL20180126: ACTUALLY NOT NEEDED
      mesh.getFaceElement(labelc, labelElementA, labelElementB);
      mesh.getFaceFace(labelc,labelFaceNeighbor);
      */
      GlbRHS(3*Num0E+3*lc+2) = 0;
      GlbMat.unitifyRow(Num0E+labelc,2);
      GlbMat.unitifyColumn(Num0E+labelc,2);
    }
  }
  // Adjusting the glb.disc.lin.sys. according to bndry.cond. sym. in z
  // GlbRHS: No change
  for (int labelc=1; labelc<=Num1C; ++labelc) {
    int lc = labelc - 1;
    if (ABC(mesh.isBoundaryFace(labelc))==6) {  // For symmetry in z-direction
      /*
      // JL20180126: ACTUALLY NOT NEEDED
      mesh.getFaceElement(labelc, labelElementA, labelElementB);
      mesh.getFaceFace(labelc,labelFaceNeighbor);
      */
      GlbRHS(3*Num0E+3*lc+3) = 0;
      GlbMat.unitifyRow(Num0E+labelc,3);
      GlbMat.unitifyColumn(Num0E+labelc,3);
    }
  }

  return(0);  // If successful
}

// LinElas3d_WG_HexaQ03Q03RT03Q0_ModiLinSys1.cpp
