// JL20180518: TO BE REVISED FOR USING Array for Boundary Conditions (ABC)
// JL20150409: FOR NOW ASSUMING ALL DIRICHLET BOUNDARY CONDITIONS
// Darcy3d_WG_HexaQ0Q0RT0_AsmBndryConds.cpp
// The Schur-complement approach
// Using ABC: Array of Boundary Conditions
// James Liu, Graham Harper, ColoState; 2014/07--2018/05

#include <cmath>

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

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


int Darcy3d_WG_HexaQ0Q0RT0_AsmBndryConds(Vector &GlbVecD, Vector &GlbVecN,
                                         double (*fxnpD)(PtVec3d),
                                         double (*fxnuN)(PtVec3d),
                                         // const Vector &ABC,
                                         const HexaMesh &mesh,
                                         const GaussQuad &GQQ)
{
  PtVec3d cntr, qp;

  // Setup
  int Num1C = mesh.numberFaces();

  // JL20180518: TO BE REVISED FOR USING Array for Boundary Conditions (ABC)
  // Sorting out Dirichlet boundary faces
  /*
  int NumDiricheltBndryFcs=0, NumNeumannBndryFcs=0, NumRobinBndryFcs=0;
  int NumSymXBndryFcs=0, NumSymYBndryFcs=0, NumSymZBndryFcs=0;
  for (int labelc=1; labelc<=Num1C; ++labelc) {
    if (mesh.isBoundaryFace(labelc)==0) continue;  // Interior faces
    if (ABC(mesh.isBoundaryFace(labelc))==1) {  // Dirichlet boundary faces
      NumDiricheltBndryFcs++;
    }
    if (ABC(mesh.isBoundaryFace(labelc))==2) {  // Neumann boundary faces
      NumNeumannBndryFcs++;
    }
    if (ABC(mesh.isBoundaryFace(labelc))==3) {  // Robin boundary faces
      NumRobinBndryFcs++;
    }
    if (ABC(mesh.isBoundaryFace(labelc))==4) {  // Symmetric boundary faces
      NumSymXBndryFcs++;
    }
    if (ABC(mesh.isBoundaryFace(labelc))==4) {  // Symmetric boundary faces
      NumSymYBndryFcs++;
    }
    if (ABC(mesh.isBoundaryFace(labelc))==6) {  // Symmetric boundary faces
      NumSymZBndryFcs++;
    }
  }
  std::cout << "NumDiricheltBndryFcs=" << NumDiricheltBndryFcs << "\n";
  */
  
  // JL20150409: FOR NOW ASSUMING ALL DIRICHLET BOUNDARY CONDITIONS
  // Processing Dirichlet boundary conditions to get GlbVecD
  GlbVecD.resize(Num1C);
  for (int labelc=1; labelc<=Num1C; ++labelc) {
    if (mesh.isBoundaryFace(labelc)>0) {  // Dirichlet face
      Quadri3d quadri = mesh.face(labelc);
      /*
      for (j=0; j<4; ++j)
        std::cout << quadri.vertex(j) << "\n";
      std::cout << "\n";
      */
      double pDavg = 0, 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);  // JL20150425: REALLY NEEDED?
        PtVec3d qp(0,0,0);
        for (int j=0; j<GQQ.numberVertices(); ++j) {
          qp = qp + GQQ.baryCoordinate(k,j)*quadri.vertex(j);
        }
        double pval = fxnpD(qp);
        pDavg += GQQ.weight(k) * jac * pval;
        area += GQQ.weight(k) * jac;
      }
      pDavg /= area;
      GlbVecD(labelc) = pDavg;
    }
  }

  // Processing Neumann boundary conditions to get GlbVecN
  GlbVecN.resize(Num1C);

  return(0);  // If successful
}
