// JL20150409: FOR NOW: ASSUMING ALL DIRICHLET BOUNDARY CONDITIONS
// Darcy3d_WG_TetraP0P0RT0_AsmBndryConds.cpp
// The Schur-complement approach
// James Liu, Graham Harper, ColoState; 2014/07--2017/03

#include <cmath>

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

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


int Darcy3d_WG_TetraP0P0RT0_AsmBndryConds(Vector &GlbVecD, Vector &GlbVecN,
                                          double (*fxnpD)(PtVec3d),
                                          double (*fxnuN)(PtVec3d),
                                          const TetraMesh &mesh,
                                          const GaussQuad &GQT)
{
  int j, k;
  double pDavg, pDval;
  PtVec3d cntr, qp;
  Tri3d tri;

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

  // Processing Dirichlet boundary conditions to get GlbVecD
  GlbVecD.resize(Num1C);
  for (int labelc=1; labelc<=Num1C; ++labelc) {
    if (mesh.isBoundaryFace(labelc)>0) {  // Dirichlet face
      tri = mesh.face(labelc);
      pDavg = 0;
      for (k=0; k<GQT.numberQuadraturePoints(); ++k) {
        qp = PtVec3d(0,0,0);
        for (j=0; j<GQT.numberVertices(); ++j) {
          qp = qp + GQT.baryCoordinate(k,j) * tri.vertex(j);
        }
        pDval = fxnpD(qp);
        pDavg += GQT.weight(k) * pDval;
      }
      GlbVecD(labelc) = pDavg;
    }
	}

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

  return(0);  // If successful
}