// Darcy3d_WG_TetraP0P0RT0_ProjPres.cpp
// Projecting pressure
// James Liu, ColoState; 2014/10--2017/01

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


int Darcy3d_WG_TetraP0P0RT0_ProjPres(Vector &ProjPresEm, Vector &ProjPresFc,
                                     double (*fxnp)(PtVec3d),
                                     double (*fxnpD)(PtVec3d),
                                     const TetraMesh &mesh,
                                     const GaussQuad &GQTe,
                                     const GaussQuad &GQT)
{
  int j, k, labelc, labele;
  double pavg, pDavg, pval;
  double x,y,z;
  Tetra Te;
  Tri3d T;
  PtVec3d cntr, qp;

  // Setup
  int Num0E = mesh.numberElements();
  int Num1C = mesh.numberFaces();
  ProjPresEm.resize(Num0E);
  ProjPresFc.resize(Num1C);

  // Computing elementwise pressure averages
  for (labele=1; labele<=mesh.numberElements(); ++labele) {
    Te = mesh.element(labele);
    pavg = 0;
    for (k=0; k<GQTe.numberQuadraturePoints(); ++k) {
      qp = PtVec3d(0,0,0);
      for (j=0; j<GQTe.numberVertices(); ++j) {
        qp = qp + GQTe.baryCoordinate(k,j) * Te.vertex(j);
      }
      pval = fxnp(qp);
      pavg = pavg + GQTe.weight(k) * pval;
    }
    ProjPresEm(labele) = pavg;
  }

  // Computing facewise pressure averages
  for (labelc=1; labelc<=mesh.numberFaces(); ++labelc) {
    T = 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) * T.vertex(j);
      }
      pval = fxnpD(qp);
      pDavg = pDavg + GQT.weight(k) * pval;
    }
    ProjPresFc(labelc) = pDavg;
  }

  return(0);  // if successful
}