// LinElas3d_CG_TetraP13_AsmSource.cpp
// James Liu, ColoState; 2014/07--2018/02

#include "vector.h"

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


// JL20160901: TO BE REVISED FOR EFFECIENCY

int LinElas3d_CG_TetraP13_AsmSource(Vector &GlbVecSource,
                                    PtVec3d (*fxnf)(PtVec3d),
                                    const TetraMesh &mesh,
                                    const GaussQuad &GQTe)
{
  int labelVertex[4];
  PtVec3d qp;
  Vector EltRHS(12);

  // Setup
  int DOFs = 3*mesh.numberNodes();
  GlbVecSource.resize(DOFs);

  // Preparation
  PtVec3d *data = new PtVec3d[GQTe.numberQuadraturePoints()];
  
  // Computing and assembling
  for (int labele=1; labele<=mesh.numberElements(); ++labele) {
    mesh.getElementNode(labele, labelVertex);
    Tetra tetra = mesh.element(labele);

    for (int k=0; k<GQTe.numberQuadraturePoints(); ++k) {
      qp = PtVec3d(0,0,0);
      for (int j=0; j<GQTe.numberVertices(); ++j) {
        qp = qp + GQTe.baryCoordinate(k,j) * tetra.vertex(j);
      }
      data[k] = fxnf(qp);
    }

    // Element RHS for source or Element-level linear form
    EltRHS = FE_Tetra_CP13_IntgrlBasFxnData(tetra, data, GQTe);

    // Now assembling
    for (int j=0; j<4; ++j) {
      int k = labelVertex[j];
      GlbVecSource(3*(k-1)+1) += EltRHS(3*j+1);
      GlbVecSource(3*(k-1)+2) += EltRHS(3*j+2);
      GlbVecSource(3*(k-1)+3) += EltRHS(3*j+3);
    }
  }

  // Releasing
  delete[] data;

  return(0);  // If successful
}

// LinElas3d_CG_TetraP13_AsmSource.cpp
