// JL20171105: BEING REVISED FOR GENERAL TREATMENT OF BOUNDARY CONDITIONS
// LinElas3d_WG_HexaQ03Q03RT03Q0_AsmBndryConds.cpp
// Solving linear elasticity in 3d by WG(Q0^3,Q0^3;RT[0]^3,Q0) on a hexa.mesh
// The Schur-complement approach
// James Liu, Graham Harper, ColoState; 2014/07--2017/11

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

#include "GaussQuad.h"
#include "PtVec3d.h"
#include "HexaMesh.h"


int LinElas3d_WG_HexaQ03Q03RT03Q0_AsmBndryConds(
  Vector &GlbVecD, Vector &GlbVecN, Vector &flag,
  PtVec3d (*fxnuD)(PtVec3d), PtVec3d (*fxntN)(PtVec3d),
  const Vector &ABC, const HexaMesh &mesh, const GaussQuad &GQQ)
{
  PtVec3d qp, uDavg;
  
  // Setup
  int Num1C = mesh.numberFaces();
  int DOFsFC = 3*Num1C;
  flag.resize(DOFsFC);

  // Sorting out boundary faces
  int NumBndryFcsDirichlet=0, NumBndryFcsNeumann=0, NumBndryFcsRobin=0;
  int NumBndryFcsSymX=0, NumBndryFcsSymY=0, NumBndryFcsSymZ=0;
  for (int labelc=1; labelc<=Num1C; ++labelc) {
    if (mesh.isBoundaryFace(labelc)==0) continue;  // Interior faces
    if (ABC(mesh.isBoundaryFace(labelc))==1) NumBndryFcsDirichlet++;
    if (ABC(mesh.isBoundaryFace(labelc))==2) NumBndryFcsNeumann++;
    if (ABC(mesh.isBoundaryFace(labelc))==3) NumBndryFcsRobin++;
    if (ABC(mesh.isBoundaryFace(labelc))==4) NumBndryFcsSymX++;
    if (ABC(mesh.isBoundaryFace(labelc))==5) NumBndryFcsSymY++;
    if (ABC(mesh.isBoundaryFace(labelc))==6) NumBndryFcsSymZ++;
  }
  std::cout << "NumBndryFcsDirichlet=" << NumBndryFcsDirichlet << "\n";
  int *LabelFaceDirichlet = new int[NumBndryFcsDirichlet];
  int *LabelFaceSymX = new int[NumBndryFcsSymX];
  int *LabelFaceSymY = new int[NumBndryFcsSymY];
  int *LabelFaceSymZ = new int[NumBndryFcsSymZ];
  NumBndryFcsDirichlet=0;  NumBndryFcsNeumann=0;  NumBndryFcsRobin=0;
  NumBndryFcsSymX=0, NumBndryFcsSymY=0, NumBndryFcsSymZ=0;
  for (int labelc=1; labelc<=Num1C; ++labelc) {
    if (mesh.isBoundaryFace(labelc)==0) continue;  // Interior faces
    if (ABC(mesh.isBoundaryFace(labelc))==1) {  // Boundary faces: Dirichlet
      LabelFaceDirichlet[NumBndryFcsDirichlet] = labelc;
      NumBndryFcsDirichlet++;
    }
    if (ABC(mesh.isBoundaryFace(labelc))==4) {  // Boundary faces: Symmetry in x-direction
      LabelFaceSymX[NumBndryFcsSymX] = labelc;
      NumBndryFcsSymX++;
    }
    if (ABC(mesh.isBoundaryFace(labelc))==5) {  // Boundary faces: Symmetry in y-direction
      LabelFaceSymY[NumBndryFcsSymY] = labelc;
      NumBndryFcsSymY++;
    }
    if (ABC(mesh.isBoundaryFace(labelc))==6) {  // Boundary faces: Symmetry in z-direction
      LabelFaceSymZ[NumBndryFcsSymZ] = labelc;
      NumBndryFcsSymZ++;
    }
  }

  // Processing Dirichlet boundary conditions
  GlbVecD.resize(DOFsFC);
  for (int ic=0; ic<NumBndryFcsDirichlet; ++ic) {
    int labelc = LabelFaceDirichlet[ic];
    int lc = LabelFaceDirichlet[ic] - mesh.beginLabelFace();
    Quadri3d quadri = mesh.face(labelc);
    PtVec3d uDavg(0,0,0);
    double area = 0;
    for (int k=0; k<GQQ.numberQuadraturePoints(); ++k) {
      double xhat = GQQ.CartesianCoordinate(k,0);
      double yhat = GQQ.CartesianCoordinate(k,1);
      double jac = quadri.JacobianDeterminant(xhat, yhat);
      jac = fabs(jac);  // Necessary?
      qp = PtVec3d(0,0,0);
      for (int j=0; j<GQQ.numberVertices(); ++j) {
        qp = qp + GQQ.baryCoordinate(k,j)*quadri.vertex(j);
      }
      PtVec3d uval = fxnuD(qp);
      uDavg = uDavg + GQQ.weight(k) * jac * uval;
      area += GQQ.weight(k) * jac;
    }
    uDavg = (1/area)*uDavg;
    GlbVecD(3*lc+1) = uDavg.xCrd();
    GlbVecD(3*lc+2) = uDavg.yCrd();
    GlbVecD(3*lc+3) = uDavg.zCrd();
    flag(3*lc+1) = 1;
    flag(3*lc+2) = 1;
    flag(3*lc+3) = 1;
  }

  // Processing Neumann boundary conditions to get GlbVecN
  // None for now
  GlbVecN.resize(DOFsFC);

  // Cleanup
  delete[] LabelFaceDirichlet;
  delete[] LabelFaceSymX, LabelFaceSymY, LabelFaceSymZ;

  return(0);  // If successful
}

// LinElas3d_WG_HexaQ03Q03RT03Q0_AsmBndryConds.cpp
