// EqnBC_LinElas3d_Ex11sincoscos.h
// A self-made simple example on (0,1)^3
// Dilation fully visible, displacement components partially visible on bndry.
// lambda=1, mu=1
// James Liu, Graham Harper, ColoState; 2014/07--2018/02

#include <cmath>
#include "mat3.h"
#include "PtVec3d.h"

#ifndef PI
#define PI 3.141592653589793
#endif

#ifndef lambda
#define lambda  1
#endif

#ifndef mu
#define mu  1
#endif


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

// The right-hand side function in the linear elasticity equation

PtVec3d fxnf(PtVec3d pt)
{
  double x = pt.xCrd();
  double y = pt.yCrd();
  double z = pt.zCrd();
  PtVec3d f(sin(PI*x) * cos(PI*y) * cos(PI*z),
            cos(PI*x) * sin(PI*y) * cos(PI*z),
            cos(PI*x) * cos(PI*y) * sin(PI*z));
  return (PI*(lambda+2*mu))*f;
}


// Dirichlet boundary conditions

PtVec3d fxnuD(PtVec3d pt)
{
  double x = pt.xCrd();
  double y = pt.yCrd();
  double z = pt.zCrd();
  PtVec3d uD(sin(PI*x) * cos(PI*y) * cos(PI*z),
             cos(PI*x) * sin(PI*y) * cos(PI*z),
             cos(PI*x) * cos(PI*y) * sin(PI*z));
  return 1.0/(3*PI) * uD;
}


// JL20170427: None for now 
// Neumann boundary conditions

PtVec3d fxntN(PtVec3d pt)
{
  return PtVec3d(0,0,0);
}


// Function for the known displacement

PtVec3d fxnu(PtVec3d pt)
{
  double x = pt.xCrd();
  double y = pt.yCrd();
  double z = pt.zCrd();
  PtVec3d u(sin(PI*x) * cos(PI*y) * cos(PI*z),
            cos(PI*x) * sin(PI*y) * cos(PI*z),
            cos(PI*x) * cos(PI*y) * sin(PI*z));
  return 1.0/(3*PI) * u;
}


// JL20170519: TO BE COMPLETED BY Graham
// Function for the stress tensor sigma

Mat3 fxnsigma(PtVec3d pt)
{
  double x = pt.xCrd();
  double y = pt.yCrd();
  double z = pt.zCrd();
  //
  Mat3 sigma;
  //
  sigma(1,1) = (3*lambda+2*mu) * cos(PI*x) * cos(PI*y) * cos(PI*z);
  sigma(1,2) =         (-2*mu) * sin(PI*x) * sin(PI*y) * cos(PI*z);
  sigma(1,3) =         (-2*mu) * sin(PI*x) * cos(PI*y) * sin(PI*z);
  //
  // Since sigma is symmetric, do not recompute these values
  sigma(2,1) = sigma(1,2);
  sigma(2,2) = (3*lambda+2*mu) * cos(PI*x) * cos(PI*y) * cos(PI*z);
  sigma(2,3) =         (-2*mu) * cos(PI*x) * sin(PI*y) * sin(PI*z);
  //
  sigma(3,1) = sigma(1,3);
  sigma(3,2) = sigma(2,3);
  sigma(3,3) = (3*lambda+2*mu) * cos(PI*x) * cos(PI*y) * cos(PI*z);
  //

  return (1.0/3)*sigma;
}


// Function for the known dilation or displacement divergence

double fxndivu(PtVec3d pt)

{
  double x = pt.xCrd();
  double y = pt.yCrd();
  double z = pt.zCrd();
  double divu = cos(PI*x) * cos(PI*y) * cos(PI*z);
  return divu;
}

// EqnBC_LinElas3d_Ex11sincoscos.cpp
