// LinElas3d_WG_HexaQ03Q03RT03Q0_AsmSource1.cpp
// Solving linear elasticity in 3d by WG(Q0^3,Q0^3;RT[0]^3,Q0) on a hexa.mesh
// The single-matrix approach
// James Liu, Graham Harper, ColoState; 2014/07--2017/05

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

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


int LinElas3d_WG_HexaQ03Q03RT03Q0_AsmSource1(Vector &GlbVecSource,
                                             PtVec3d (*fxnf)(PtVec3d),
                                             const HexaMesh &mesh,
                                             const GaussQuad &GQH)
{
  // Setup
  int Num0E = mesh.numberElements();
  int Num1C = mesh.numberFaces();
  int DOFs = 3*Num0E + 3*Num1C;
  GlbVecSource.resize(DOFs);

  // Assembling the source term 
  for (int labele=1; labele<=Num0E; ++labele) {
	  int le = labele - 1;
    Hexa hexa = mesh.element(labele);
    hexa.enrich();
    PtVec3d fintgrl(0,0,0);
    for (int k=0; k<GQH.numberQuadraturePoints(); ++k) {
      double xhat = GQH.CartesianCoordinate(k,0);
      double yhat = GQH.CartesianCoordinate(k,1);
      double zhat = GQH.CartesianCoordinate(k,2);
      double jac = hexa.JacobianDeterminant(xhat, yhat, zhat);
      jac = fabs(jac);  // JL20150422: DON'T USE fabs?
      /*
      PtVec3d qp(0,0,0);
      for (int j=0; j<GQH.numberVertices(); ++j) {
        qp = qp + GQH.baryCoordinate(k,j)*hexa.vertex(j);
      }
      */
      PtVec3d qp = hexa.trilinearmapping(xhat,yhat,zhat);
      PtVec3d fval = fxnf(qp);
      fintgrl = fintgrl + (GQH.weight(k)*jac) * fval;
    }
    GlbVecSource(3*le+1) = fintgrl.xCrd();
    GlbVecSource(3*le+2) = fintgrl.yCrd();
    GlbVecSource(3*le+3) = fintgrl.zCrd();
  }

  return(0);  // If successful
}

// LinElas3d_WG_HexaQ03Q03RT03Q0_AsmSource1.cpp
