// JL20161015: TO BE FINISHED
// LinElas3d_CG_TetraP13_DsplDivStrs.cpp
// Postprocessing to get elementwise displacement, divergence, and stress
// James Liu, ColoState; 2014/07--2016/10

#include "vector.h"

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


// JL20160901: TO BE REVISED FOR EFFECIENCY

int LinElas3d_CG_TetraP13_DsplDivStrs(PtVec3d *DsplEm, double *DivEm,
                                      Mat3 *StrsEm, double lambda, double mu,
                                      const TetraMesh &mesh, Vector &sln)
{
  int labelVertex[4];
  double cof[12], BasFxnDiv[12], tmpx, tmpy, tmpz;
  PtVec3d BasFxnVal[12];
  Mat3 BasFxnStrs[12];

  // Computing displacement elementwise average
  // DsplEm: Type = PtVec3d;  Size = mesh.numberElements()
  for (int labele=1; labele<=mesh.numberElements(); ++labele) {
    int le = labele - 1;
    mesh.getElementNode(labele, labelVertex);
    
    // The 12 coefficients for CG P1^3 basis functions
    for (int i=0; i<4; ++i)
      for (int j=0; j<3; ++j)
        cof[3*i+j] = sln[3*(labelVertex[i]-1)+j];

    //
    tmpx = 0;  tmpy = 0;  tmpz = 0;
    for (int i=0; i<4; ++i) {
      tmpx = tmpx + 0.25 * cof[3*i+0];
      tmpy = tmpy + 0.25 * cof[3*i+1];
      tmpz = tmpz + 0.25 * cof[3*i+2];
    }
    DsplEm[le] = PtVec3d(tmpx,tmpy,tmpz);
  }

  // JL20161015: TO BE FINISHED
  // Computing divergence elementwise (constant)
  // DivEm: Type = double;  Size = mesh.numberElements()
  for (int labele=1; labele<=mesh.numberElements(); ++labele) {
    int le = labele - 1;
    DivEm[le] = 0;
  }

  // JL20161015: TO BE FINISHED
  // Computing stress elementwise (constant)
  // StrsEm: Type = Mat3;  Size = mesh.numberElements()
  for (int labele=1; labele<=mesh.numberElements(); ++labele) {
    int le = labele - 1;
    StrsEm[le].reset();
  }

  return(0);  // If successful
}