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

// Brick: Coordinates P0(x0,y0,z0), P1(x1,y1,z1)
// RT[0] normalized Basis:
//   [1;0;0], [0;1;0], [0;0;1], [X;0;0], [0;Y;0], [0;0;Z]
//   where X=x-xc, Y=y-yc, Z=z-zc, and (xc,yc,zc) is the element center

#include "matrix.h"

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


////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////


// WG3d: Brick(Q0,Q0,RT[0]): Element mass matrix (actually just a number)

void WG3d_BrickQ0Q0RT0_EltMassMat(
  double x0, double y0, double z0,
  double x1, double y1, double z1,
  double EMM[][1])
{
  // Actually = 1*volume
  double Deltax = x1 - x0;
  double Deltay = y1 - y0;
  double Deltaz = z1 - z0;
  EMM[0][0] = Deltax * Deltay * Deltaz;
}


// WG3d: Brick(Q0,Q0,RT0): Coeff. of disc.wk.grad. in RT[0] nmlz.bas.

void WG3d_BrickQ0Q0RT0_CofDiscWkGrad_NmlzBas(
  double x0, double y0, double z0,
  double x1, double y1, double z1,
  double CDWG[][6])
{
  double Deltax, Deltay, Deltaz;
  double Deltax2, Deltay2, Deltaz2;
  // double vol;

  Deltax = x1 - x0;
  Deltay = y1 - y0;
  Deltaz = z1 - z0;

  Deltax2 = Deltax * Deltax;
  Deltay2 = Deltay * Deltay;
  Deltaz2 = Deltaz * Deltaz;

  // vol = Deltax * Deltay * Deltaz;

  // double CDWG[7][6];
  // 1st index:
  //   0 for the element interior bas.fxn.
  //   1-6 for face bas.fxn.: 1,2 for x-faces; 3,4 for y-faces; 5,6 for z-faces
  // 2nd index: 0-5 for the six normalized basis function of RT[0]

  for (int i=0; i<7; ++i)
    for (int j=0; j<6; ++j)
      CDWG[i][j] = 0.0;

  CDWG[0][3] = -12/Deltax2;  CDWG[0][4] = -12/Deltay2;  CDWG[0][5] = -12/Deltaz2;

  CDWG[1][0] = -1/Deltax;  CDWG[1][3] =  6/Deltax2;
  CDWG[2][0] =  1/Deltax;  CDWG[2][3] =  6/Deltax2;

  CDWG[3][1] = -1/Deltay;  CDWG[3][4] =  6/Deltay2;
  CDWG[4][1] =  1/Deltay;  CDWG[4][4] =  6/Deltay2;

  CDWG[5][2] = -1/Deltaz;  CDWG[5][5] =  6/Deltaz2;
  CDWG[6][2] =  1/Deltaz;  CDWG[6][5] =  6/Deltaz2;

  return;
}


// WG3d: Brick(Q0,Q0,RT[0]): Element grad-grad matrix

void WG3d_BrickQ0Q0RT0_EltGradGradMat(double EGGM[][7],
  double x0, double y0, double z0,
  double x1, double y1, double z1)
{
  int i, j;
  double Deltax, Deltay, Deltaz;
  double Deltax2, Deltay2, Deltaz2;
  double vol;

  Deltax = x1 - x0;
  Deltay = y1 - y0;
  Deltaz = z1 - z0;

  Deltax2 = Deltax * Deltax;
  Deltay2 = Deltay * Deltay;
  Deltaz2 = Deltaz * Deltaz;

  vol = Deltax * Deltay * Deltaz;
  
  // std::cout << "Deltax,Deltay,Deltaz=" << Deltax << " " << Deltay << " " << Deltaz << "\n";

  // double EGGM[7][7];
  // 1st or 2nd index:
  //   0 for the element interior bas.fxn.
  //   1-6 for face bas.fxn.: 1,2 for x-faces; 3,4 for y-faces; 5,6 for z-faces
  for (i=0; i<7; ++i)
    for (j=0; j<7; ++j)
      EGGM[i][j] = 0.0;

  EGGM[0][0] =  vol*12.0*(1/Deltax2 + 1/Deltay2 + 1/Deltaz2);
  EGGM[0][1] = -vol*6/Deltax2;  EGGM[0][2] = EGGM[0][1];
  EGGM[0][3] = -vol*6/Deltay2;  EGGM[0][4] = EGGM[0][3];
  EGGM[0][5] = -vol*6/Deltaz2;  EGGM[0][6] = EGGM[0][5];
  for (i=1; i<=6; i++)  EGGM[i][0] = EGGM[0][i];

  EGGM[1][1] =  vol*4/Deltax2;  EGGM[1][2] =  vol*2/Deltax2;
  EGGM[2][1] =  EGGM[1][2];     EGGM[2][2] =  vol*4/Deltax2;

  EGGM[3][3] =  vol*4/Deltay2;  EGGM[3][4] =  vol*2/Deltay2;
  EGGM[4][3] =  EGGM[3][4];     EGGM[4][4] =  vol*4/Deltay2;

  EGGM[5][5] =  vol*4/Deltaz2;  EGGM[5][6] =  vol*2/Deltaz2;
  EGGM[6][5] =  EGGM[5][6];     EGGM[6][6] =  vol*4/Deltaz2;

  /*
  for (i=0; i<7; ++i) {
    for (j=0; j<7; ++j) {;
      std::cout << "EGGM[i][j]=" << EGGM[i][j] << "  ";
    }
    std::cout << "\n";
  }
  */

  return;
}


////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////


// WG3d: Brick(Q0,Q0,RT[0]): 6x7 matrix of coeff.disc.wk.grad. in RT[0] nmlz.bas.

FullMatrix WG3d_BrickQ0Q0RT0_CofDiscWkGrad_NmlzBas(const Brick &brick)
{
  double Deltax, Deltay, Deltaz;
  double Deltax2, Deltay2, Deltaz2;
  PtVec3d P0, P1;
  FullMatrix CDWG(6,7);
  // 1st index: 1-6 for the six normalized basis function of RT[0]
  // 2nd index:
  //   1 for the element interior bas.fxn.
  //   2-7 for face bas.fxn.: 2,3 for x-faces; 4,5 for y-faces; 6,7 for z-faces

  brick.getVertices(P0, P1);
  
  Deltax = P1.xCrd() - P0.xCrd();
  Deltay = P1.yCrd() - P0.yCrd();
  Deltaz = P1.zCrd() - P0.zCrd();
  
  Deltax2 = Deltax * Deltax;
  Deltay2 = Deltay * Deltay;
  Deltaz2 = Deltaz * Deltaz;
  
  
  CDWG(4,1) = -12/Deltax2;
  CDWG(5,1) = -12/Deltay2;
  CDWG(6,1) = -12/Deltaz2;
  
  CDWG(1,2) = -1/Deltax;
  CDWG(4,2) =  6/Deltax2;
  
  CDWG(1,3) =  1/Deltax;
  CDWG(4,3) =  6/Deltax2;
  
  CDWG(2,4) = -1/Deltay;
  CDWG(5,4) =  6/Deltay2;
  
  CDWG(2,5) =  1/Deltay;
  CDWG(5,5) =  6/Deltay2;
  
  CDWG(3,6) = -1/Deltaz;
  CDWG(6,6) =  6/Deltaz2;
  
  CDWG(3,7) =  1/Deltaz;
  CDWG(6,7) =  6/Deltaz2;
  
  return CDWG;
}


// WG3d: Brick(Q0,Q0,RT[0]): Element grad-grad matrix

FullMatrix WG3d_Brick0_EltGradGradMatK(const Brick &brick, Mat3 &MatK)
{
  int i, j, k, l;
  FullMatrix CDWG(6,7), EltMat(7,7), GMK(6,6);

  CDWG = WG3d_BrickQ0Q0RT0_CofDiscWkGrad_NmlzBas(brick);
  GMK = Hdiv_BrickRT0_GramMatK_NmlzBas(brick, MatK);

  // JL20150619: NEEDS DOUBLE-CHECK!!!
  for (i=1; i<=7; ++i) {
    for (j=1; j<=7; ++j) {
      EltMat(i,j) = 0;
      for (k=1; k<=6; ++k)
        for (l=1; l<=6; ++l)
          EltMat(i,j) += CDWG(k,i) * GMK(k,l) * CDWG(l,j);
    }
  }
  
  return EltMat;
}


// WG3d: Brick(Q0,Q0,RT[0]): Normal fluxes of 7 WG basis functions

FullMatrix WG3d_Brick0_NmlFluxK_NmlzBas_Proj(const Brick &brick, Mat3 &MatK)
{
  int i, j, k;
  FullMatrix CDWG(6,7), NmlFluxProjRT0(6,6), NmlFluxRT0(6,6), NmlFluxWG(6,7);

  CDWG = WG3d_BrickQ0Q0RT0_CofDiscWkGrad_NmlzBas(brick);
  NmlFluxRT0 = Hdiv_BrickRT0_NmlFlux_NmlzBas(brick);

  for (i=1; i<=6; ++i) {
    NmlFluxProjRT0(i,1)
    = MatK(1,1)*NmlFluxRT0(i,1) + MatK(2,1)*NmlFluxRT0(i,2) + MatK(3,1)*NmlFluxRT0(i,3);

    NmlFluxProjRT0(i,2)
    = MatK(1,2)*NmlFluxRT0(i,1) + MatK(2,2)*NmlFluxRT0(i,2) + MatK(3,2)*NmlFluxRT0(i,3);

    NmlFluxProjRT0(i,3)
    = MatK(1,3)*NmlFluxRT0(i,1) + MatK(2,3)*NmlFluxRT0(i,2) + MatK(3,3)*NmlFluxRT0(i,3);

    NmlFluxProjRT0(i,4) = MatK(1,1) * NmlFluxRT0(i,4);
    NmlFluxProjRT0(i,5) = MatK(2,2) * NmlFluxRT0(i,5);
    NmlFluxProjRT0(i,6) = MatK(3,3) * NmlFluxRT0(i,6);
  }

  // JL20150620: NOTICE THE NEGATIVE SIGN!!
  NmlFluxWG = (-1) * (NmlFluxProjRT0 * CDWG);
  /*
  for (i=1; i<=6; ++i) {
    for (j=1; j<=7; ++j) {
      NmlFluxWG(i,j) = 0;
      for (k=1; k<=6; ++k)
        NmlFluxWG(i,j) -= NmlFluxProjRT0(i,k) * CDWG(k,j);
    }
  }
  */

  return NmlFluxWG;
}

// WG3d_Brick.cpp