// Darcy3d_WG_TetraP0P0RT0_Err.cpp
// JL20150629: Using quadratures independent of other subroutines
// James Liu, ColoState; 2014/07--2017/02

#include <cmath>

#include "vector.h"

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


int Darcy3d_WG_TetraP0P0RT0_Err(double &L2ErrPres, double &L2ErrVel,
                                const TetraMesh &mesh,
                                const Vector &NumerPresEm,
                                const FullMatrix &NumerVelCofRT0,
                                double (*fxnp)(PtVec3d),
                                PtVec3d (*fxnu)(PtVec3d))
//                              const GaussQuad &GQH, const GaussQuad &GQQ)
{
  int i, j, k, labele, le;
  int labelFace[4], sign[4];
  double ferr, flux, jac, perr, pval, sa, sum, sump, sumv, vol, vtmp;
  double x, y, z, xc, yc, zc, X, Y, Z;
  PtVec3d cntr, nml, NumerVel, qp, vel, verr;
  PtVec3d w[5];  // w[0] not used
  Tri3d tri;
  Tetra tetra;

  // Using its own quadratures
  GaussQuad GQTe, GQT;
  GQTe.setForTetrahedron(11);
  GQT.setForTriangle(7);
  
  // Setup
  int Num0E = mesh.numberElements();
  int Num1C = mesh.numberFaces();
  Vector ErrPresEm(Num0E);
  Vector ErrVelEm(Num0E);
  w[1] = PtVec3d(1,0,0);
  w[2] = PtVec3d(0,1,0);
  w[3] = PtVec3d(0,0,1);

  // Computing L2-norm of errors in pressure, velocity 
  for (le=0; le<Num0E; ++le) {
    labele = le + mesh.beginLabelElement();
    tetra = mesh.element(labele);
 
    cntr = tetra.center();
    xc = cntr.xCrd();
    yc = cntr.yCrd();
    zc = cntr.zCrd();
    vol = tetra.volume();

    sump = 0;
    sumv = 0;
    for (i=0; i<GQTe.numberQuadraturePoints(); ++i) {
      qp = PtVec3d(0,0,0);
      for (j=0; j<GQTe.numberVertices(); ++j) {
       qp = qp + GQTe.baryCoordinate(i,j) * tetra.vertex(j);
      }
      pval = fxnp(qp);
      perr = pval - NumerPresEm(labele);
      sump += GQTe.weight(i) * (perr*perr);
      //
      x = qp.xCrd();  y = qp.yCrd();  z = qp.zCrd();
      X = x - xc;  Y = y - yc;  Z = z - zc;
      w[4] = PtVec3d(X,Y,Z);
      PtVec3d NV(0,0,0);
      for (j=1; j<=4; ++j) {
        NV = NV + NumerVelCofRT0(labele,j) * w[j];
	  }
      vel = fxnu(qp);
      verr = vel - NV;
      vtmp = verr.l2norm();
      sumv += GQTe.weight(i) * (vtmp*vtmp);
    }
    ErrPresEm(labele) = sqrt(sump*vol);
    ErrVelEm(labele) = sqrt(sumv*vol);
  }
  L2ErrPres = ErrPresEm.l2norm();
  L2ErrVel = ErrVelEm.l2norm();

  // Computing the discrete L2-norm of the pressure error
  /*
  Vector discL2ErrPresElt(Num0E);
  for (le=0; le<Num0E; ++le) {
    labele = le + mesh.beginLabelElement();
    tetra = mesh.element(labele);
    cntr = tetra.center();
    vol = tetra.volume();
    pval = fxnp(cntr);
    perr = pval - NumerPresEm(labele);
    discL2ErrPresElt(labele) = sqrt(vol) * fabs(perr);
  }
  discL2ErrPres = discL2ErrPresElt.l2norm();
  */
  
  // Computing the discrete L2-norm of the normal flux error
  // Just using the analytical and numerical velocity normal components
  /*
  sum = 0;
  for (le=0; le<Num0E; ++le) {
    labele = le + mesh.beginLabelElement();
    mesh.getElementFace(labele, labelFace);
    
    tetra = mesh.element(labele);
    vol = tetra.volume();
    cntr = tetra.center();
    xc = cntr.xCrd();  yc = cntr.yCrd();  zc = cntr.zCrd();

    for (i=0; i<4; ++i) {
      tri = mesh.face(labelFace[i]);
      ferr = 0;
      for (k=0; k<GQTe.numberQuadraturePoints(); ++k) {
        nml = tri.normal();
        qp = PtVec3d(0,0,0);
        for (j=0; j<GQTe.numberVertices(); ++j) {
          qp = qp + GQTe.baryCoordinate(i,j) * tetra.vertex(j);
        }
        x = qp.xCrd();  y = qp.yCrd();  z = qp.zCrd();
        X = x - xc;  Y = y - yc;  Z = z - zc;
        w[4] = PtVec3d(X,Y,Z);
        NumerVel = PtVec3d(0,0,0);
        for (j=1; j<=4; ++j)
          NumerVel = NumerVel + NumerVelCofRT0(labele,j) * w[j];
        vel = fxnu(qp);
        verr = vel - NumerVel;
        ferr += dotProduct(verr,nml);
      }
      sum += (ferr*ferr) * vol;
    }
  }
  discL2ErrFlux = sqrt(sum);
  */
  
  return(0);  // if successful
}