// LinElas3d_CG_TetraP13_ProjDsplDivStrs.cpp
// Computing projected displacement, divergence, stress
// James Liu, ColoState; 2014/07--2016/10

// #include <cmath>

#include "vector.h"

#include "cell3d.h"
#include "GaussQuad.h"
#include "mat3.h"
#include "PtVec3d.h"
#include "TetraMesh.h"
#include "TetraMesh_ConvCP1toP0.h"


// JL20160901: TO BE REVISED FOR EFFECIENCY

int LinElas3d_CG_TetraP13_ProjDsplDivStrs(PtVec3d *DsplEmAvg, Vector &DivEmAvg,
                                          double lambda, double mu,
                                          PtVec3d (*fxnu)(PtVec3d),
                                          double (*fxndivu)(PtVec3d),
                                          Mat3 (*fxnsigma)(PtVec3d),
                                          const TetraMesh &mesh,
                                          const GaussQuad &GQTe)
{
  int labelNode[4];
  
  // Setup
  // DsplEmAvg: Size: mesh.numberElements()
  DivEmAvg.resize(mesh.numberElements());

  // Interpolating exact displacement at nodes
  PtVec3d *DsplNd = new PtVec3d[mesh.numberNodes()];
  for (int labeld=1; labeld<=mesh.numberNodes(); ++labeld) {
    DsplNd[labeld-1] = fxnu(mesh.node(labeld));
  }
  // Converting
  for (int labele=1; labele<=mesh.numberElements(); ++labele) {
    int le = labele - 1;
    mesh.getElementNode(labele, labelNode);
    DsplEmAvg[le] = PtVec3d(0,0,0);
    for (int j=0; j<4; ++j){
      DsplEmAvg[le] = DsplEmAvg[le] + 0.25 * DsplNd[labelNode[j]-1];
    }
  }
  // Freeing...
  delete[] DsplNd;
  
  // Interpolating exact divergence at nodes
  Vector DivNd(mesh.numberNodes());
  for (int labeld=1; labeld<=mesh.numberNodes(); ++labeld) {
    DivNd(labeld) = fxndivu(mesh.node(labeld));
  }
  // Converting
  TetraMesh_ConvCP1toP0(DivEmAvg, DivNd, mesh);

  return(0);  // If successful
}