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

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

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


int LinElas3d_WG_HexaQ03Q03RT03Q0_AsmSource(Vector &GlbVecE,
                                            PtVec3d (*fxnf)(PtVec3d),
                                            const HexaMesh &mesh,
                                            const GaussQuad &GQH)
{
  // Setup
  int Num0E = mesh.numberElements();
  int DOFs = 3*Num0E;
  GlbVecE.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);  // Needed?
      PtVec3d qp = hexa.trilinearmapping(xhat,yhat,zhat);
      PtVec3d fval = fxnf(qp);
      fintgrl = fintgrl + (GQH.weight(k)*jac) * fval;
    }
    GlbVecE(3*le+1) = fintgrl.xCrd();
    GlbVecE(3*le+2) = fintgrl.yCrd();
    GlbVecE(3*le+3) = fintgrl.zCrd();
  }
  
  return(0);  // If successful
}

// LinElas3d_WG_HexaQ03Q03RT03Q0_AsmSource.cpp
