// Darcy3d_WG_TetraP0P0RT0_PresVelFlux.cpp
// Computing numerical pressure, velocity, flux
// James Liu, ColoState; 2014/07--2017/02

#include <cmath>

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

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


int Darcy3d_WG_TetraP0P0RT0_PresVelFlux(Vector &NumerPresEm,
                                        FullMatrix &NumerVelCofRT0,
                                        FullMatrix &NumerFlux,
                                        const TetraMesh &mesh,
                                        Mat3 *PermK, Vector &sln,
                                        const GaussQuad &GQTe,
                                        const GaussQuad &GQT)
{
  int i, j, k, labele, le;
  int labelFace[4], sign[4];
  PtVec3d EmCntr, FcCntr, nml;
  Mat3 MatK;
  Tri3d tri, STF[4];
  Tetra tetra;
  Vector cof(5), flxtmp(4), vctmp(4);
  FullMatrix CDWG(5,4), TransCDWG(4,5);
  FullMatrix KtransMat(12,4), NumerFluxRT0(4,4), ProjMat(4,12);

  // Setup
  int Num0E = mesh.numberElements();

  // Computing/extracting numerical pressure: Elementwise constant
  NumerPresEm = sln.getSubVector(1,mesh.numberElements());

  // Computing
  // -- numerical velocity cofficients in RT0 normalized basis
  // -- normal fluxes (all elements, 4 faces per element)
  NumerVelCofRT0.resize(Num0E,4);
  NumerFlux.resize(Num0E,4);
  for (le=0; le<Num0E; ++le) {
    labele = le + mesh.beginLabelElement();
    tetra = mesh.element(labele);

	mesh.getElementFace(labele, labelFace);
    for (j=0; j<4; ++j)  STF[j] = mesh.face(labelFace[j]);

    // Determining signs by using centers and normals
    EmCntr = tetra.center();
    for (j=0; j<4; ++j) {
      sign[j] = 1;
      tri = STF[j];
      FcCntr = tri.center();
      nml = tri.normal();
      if (dotProduct(nml,FcCntr-EmCntr)<0)  sign[j] = -1;
    }

    MatK = PermK[le];
    
    cof(1) = sln(labele);
    for (k=2; k<=5; ++k)  cof(k) = sln(Num0E+labelFace[k-2]);

    WG3d_TetraP0P0RT0_CofDiscWkGrad_NmlzBas(CDWG, tetra, STF, sign, GQTe);
    TransCDWG = transpose(CDWG);

    Hdiv_Tetra01_KtransMat(KtransMat, MatK);
    Hdiv_Tetra01_ProjP1v3RT0(ProjMat, tetra, GQTe);

    vctmp = ((ProjMat*KtransMat)*TransCDWG)*cof;
    vctmp = -vctmp;
    for (j=1; j<=4; ++j)  NumerVelCofRT0(labele,j) = vctmp(j);

    Hdiv_TetraRT0_NmlFlux_NmlzBas(NumerFluxRT0, tetra);
    
    flxtmp = NumerFluxRT0 * vctmp;
    for (i=1; i<=4; ++i)  NumerFlux(labele,i) = flxtmp(i);
  }

  return(0);  // if successful
}