// Darcy3d_WG_TetraP0P0RT0_AsmBndryConds1.cpp
// Darcy: Weak Galerkin (P0,P0;RT0) on a tetrahedral mesh
// The single-matrix approach
// James Liu, ColoState; 2014/07--2017/01

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

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

int Darcy3d_WG_TetraP0P0RT0_AsmBndryConds1(Vector &GlbVecDirichlet,
                                           Vector &GlbVecNeumann,
                                           double (*fxnpD)(PtVec3d),
                                           double (*fxnuN)(PtVec3d),
                                           const TetraMesh &mesh,
                                           const GaussQuad &GQT)
{
  int i, j;
  double pDavg, pDval;
  PtVec3d qp;
  Tri3d tri;
  
  // Setup
  int Num0E = mesh.numberElements();
  int Num1C = mesh.numberFaces();
  int DOFs = Num0E + Num1C;
  GlbVecDirichlet.resize(DOFs);
  GlbVecNeumann.resize(DOFs);
  
  // Processing Dirichlet boundary conditions
  for (int labelc=1; labelc<=Num1C; ++labelc) {
    if (mesh.isBoundaryFace(labelc)>0) {
      tri = mesh.face(labelc);
      pDavg = 0;
      for (i=0; i<GQT.numberQuadraturePoints(); ++i) {
        qp = PtVec3d(0,0,0);
        for (j=0; j<GQT.numberVertices(); ++j) {
          qp = qp + GQT.baryCoordinate(i,j) * tri.vertex(j);
        }
        pDval = fxnpD(qp);
        pDavg += GQT.weight(i) * pDval;
      }
      GlbVecDirichlet(Num0E+labelc) = pDavg;
    }
  }
  
  // Processing Neumann boundary conditions
  // None for now
  
  return(0);  // if successful
}