// Darcy3d_WG_TetraP0P0RT0_AsmSource1.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_AsmSource1(Vector &GlbVecSource,
                                       double (*fxnf)(PtVec3d),
                                       const TetraMesh &mesh,
                                       const GaussQuad &GQTe)
{
  int i, j, labele, le;
  double ftmp, fval, vol;
  PtVec3d qp;
  Tetra tetra;

  // Setup
  int Num0E = mesh.numberElements();
  int Num1C = mesh.numberFaces();
  int DOFs = Num0E + Num1C;
  GlbVecSource.resize(DOFs);
  
  // Assembling GlbVecSource: For elements
  for (le=0; le<Num0E; ++le) {
    labele = le + mesh.beginLabelElement();
    tetra = mesh.element(labele);
    vol = tetra.volume();
    ftmp = 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);
      }
      fval = fxnf(qp);
      ftmp += GQTe.weight(i) * fval;
    }
    GlbVecSource(labele) = ftmp * vol;
  }

  return(0);  // if successful
}