// FE_Tetra_CP1.cpp
// James Liu, ColoState; 2014/07--2016/05

#include "matrix.h"

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


// JL20150513: TO BE REVISED FOR EFFECIENCY
// The scalar values of the 4 basis functions at a point

int FE_Tetra_CP1_BasFxnVal(double val[4], const Tetra &Te, PtVec3d pt)
{
  int i;
  double cof[4][4], dd, x[4], y[4], z[4];
  double xx, yy, zz;
  PtVec3d P[4], U[4], V[4], W[4];

  // The determinant for the tetra
  dd = 6 * Te.volume();
  
  // Primary points/vectors
  for (i=0; i<4; ++i) {
    P[i] = Te.vertex(i);
    x[i] = P[i].xCrd();
    y[i] = P[i].yCrd();
    z[i] = P[i].zCrd();
  }
  
  // Auxiliary points/vectors
  for (i=0; i<4; ++i) {
    U[i] = PtVec3d(1,y[i],z[i]);
    V[i] = PtVec3d(1,x[i],z[i]);
    W[i] = PtVec3d(1,x[i],y[i]);
  }

  // For phi_0
  cof[0][0] =  det3P(P[1],P[2],P[3])/dd;
  cof[0][1] = -det3P(U[1],U[2],U[3])/dd;
  cof[0][2] =  det3P(V[1],V[2],V[3])/dd;
  cof[0][3] = -det3P(W[1],W[2],W[3])/dd;

  // For phi_1
  cof[1][0] = -det3P(P[0],P[2],P[3])/dd;
  cof[1][1] =  det3P(U[0],U[2],U[3])/dd;
  cof[1][2] = -det3P(V[0],V[2],V[3])/dd;
  cof[1][3] =  det3P(W[0],W[2],W[3])/dd;

  // For phi_2
  cof[2][0] =  det3P(P[0],P[1],P[3])/dd;
  cof[2][1] = -det3P(U[0],U[1],U[3])/dd;
  cof[2][2] =  det3P(V[0],V[1],V[3])/dd;
  cof[2][3] = -det3P(W[0],W[1],W[3])/dd;

  // For phi_3
  cof[3][0] = -det3P(P[0],P[1],P[2])/dd;
  cof[3][1] =  det3P(U[0],U[1],U[2])/dd;
  cof[3][2] = -det3P(V[0],V[1],V[2])/dd;
  cof[3][3] =  det3P(W[0],W[1],W[2])/dd;

  //
  xx = pt.xCrd();  yy = pt.yCrd();  zz = pt.zCrd();
  for (i=0; i<4; ++i)
    val[i] = cof[i][0] + cof[i][1]*xx + cof[i][2]*yy + cof[i][3]*zz;
  
  return(0);  // if successful
}


// The constant gradient (3-dim) vectors of the 4 basis functions
// regardless of the point(s) on the tetrahedral element

int FE_Tetra_CP1_BasFxnGrad(PtVec3d BasFxnGrad[4], const Tetra &Te)
{
  int i;
  double cof[4][4], dd, x[4], y[4], z[4];
  PtVec3d P[4], U[4], V[4], W[4];

  // The determinant for the tetra
  dd = 6 * Te.volume();
  
  // Primary points/vectors
  for (i=0; i<4; ++i) {
    P[i] = Te.vertex(i);
    x[i] = P[i].xCrd();
    y[i] = P[i].yCrd();
    z[i] = P[i].zCrd();
  }
  
  // Auxiliary points/vectors
  for (i=0; i<4; ++i) {
    U[i] = PtVec3d(1,y[i],z[i]);
    V[i] = PtVec3d(1,x[i],z[i]);
    W[i] = PtVec3d(1,x[i],y[i]);
  }
  
  // For phi_0
  cof[0][1] = -det3P(U[1],U[2],U[3])/dd;
  cof[0][2] =  det3P(V[1],V[2],V[3])/dd;
  cof[0][3] = -det3P(W[1],W[2],W[3])/dd;
  
  // For phi_1
  cof[1][1] =  det3P(U[0],U[2],U[3])/dd;
  cof[1][2] = -det3P(V[0],V[2],V[3])/dd;
  cof[1][3] =  det3P(W[0],W[2],W[3])/dd;
  
  // For phi_2
  cof[2][1] = -det3P(U[0],U[1],U[3])/dd;
  cof[2][2] =  det3P(V[0],V[1],V[3])/dd;
  cof[2][3] = -det3P(W[0],W[1],W[3])/dd;
  
  // For phi_3
  cof[3][1] =  det3P(U[0],U[1],U[2])/dd;
  cof[3][2] = -det3P(V[0],V[1],V[2])/dd;
  cof[3][3] =  det3P(W[0],W[1],W[2])/dd;

  // For any point
  for (i=0; i<4; ++i)
    BasFxnGrad[i] = PtVec3d(cof[i][1], cof[i][2], cof[i][3]);

  return(0);  // if successful
}


// The 4x4 element grad-grad matrix with full diffusion/permeability matrix

FullMatrix FE_Tetra_CP1_EltGradGradMatK(const Tetra &Te, const Mat3 &MatK)
{
  int i, j;
  FullMatrix EGGMK(4,4);
  PtVec3d BFG[4];

  FE_Tetra_CP1_BasFxnGrad(BFG, Te);

  for (i=1; i<=4; ++i)
    for (j=1; j<=4; ++j)
      EGGMK(i,j) = dotProduct(BFG[i-1], MatK*BFG[j-1]) * Te.volume();

  return EGGMK;
}

// FE_Tetra_CP1.cpp