// LinElas3d_WG_HexaQ03Q03RT03Q0_AsmBndryConds1.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--2017/11

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

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


int LinElas3d_WG_HexaQ03Q03RT03Q0_AsmBndryConds1(
  Vector &GlbVecDirichlet, Vector &GlbVecNeumann, Vector &flag,
  PtVec3d (*fxnuD)(PtVec3d), PtVec3d (*fxntN)(PtVec3d), const Vector &ABC,
  const HexaMesh &mesh, const GaussQuad &GQQ)
{
  double area;
  PtVec3d qp, tNcond, tNval, uDavg, uDval;

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

  // 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";
  std::cout << "NumBndryFcsNeumann=" << NumBndryFcsNeumann << "\n";
  int *LabelFaceDirichlet = new int[NumBndryFcsDirichlet];
  int *LabelFaceNeumann = new int[NumBndryFcsNeumann];
  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))==2) {  // Boundary faces: Neumann
      LabelFaceNeumann[NumBndryFcsNeumann] = labelc;
      NumBndryFcsNeumann++;
    }
    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++;
    }
  }
  std::cout << "NumBndryFcsNeumann=" << NumBndryFcsNeumann << "\n";
  std::cout << "NumBndryFcsDirichlet=" << NumBndryFcsDirichlet << "\n";

  // Processing Dirichlet boundary conditions
  GlbVecDirichlet.resize(DOFs);
  for (int ic=0; ic<NumBndryFcsDirichlet; ++ic) {
    int labelc = LabelFaceDirichlet[ic];
    int lc = LabelFaceDirichlet[ic] - mesh.beginLabelFace();
    Quadri3d quadri = mesh.face(labelc);
    uDavg = PtVec3d(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);
      }
      uDval = fxnuD(qp);
      uDavg = uDavg + GQQ.weight(k) * jac * uDval;
      area += GQQ.weight(k) * jac;
    }
    uDavg = (1/area)*uDavg;
    GlbVecDirichlet(3*Num0E+3*lc+1) = uDavg.xCrd();
    GlbVecDirichlet(3*Num0E+3*lc+2) = uDavg.yCrd();
    GlbVecDirichlet(3*Num0E+3*lc+3) = uDavg.zCrd();
  }

  // Processing Neumann boundary conditions
  GlbVecNeumann.resize(DOFs);
  for (int ic=0; ic<NumBndryFcsNeumann; ++ic) {
    int labelc = LabelFaceNeumann[ic];
    int lc = LabelFaceNeumann[ic] - mesh.beginLabelFace();
    Quadri3d quadri = mesh.face(labelc);
    tNcond =  PtVec3d(0,0,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);
      }
      tNval = fxntN(qp);
      // std::cout << tNval << "\n";
      tNcond = tNcond + GQQ.weight(k) * jac * tNval;
    }
    GlbVecNeumann(3*Num0E+3*lc+1) = tNcond.xCrd();
    GlbVecNeumann(3*Num0E+3*lc+2) = tNcond.yCrd();
    GlbVecNeumann(3*Num0E+3*lc+3) = tNcond.zCrd();
  }

  // Clean up
  delete[] LabelFaceDirichlet, LabelFaceNeumann;
  delete[] LabelFaceSymX, LabelFaceSymY, LabelFaceSymZ;

  return(0);  // If successful
}

// LinElas3d_WG_HexaQ03Q03RT03Q0_AsmBndryConds1.cpp
