// EqnBC_LinElas3d_Mujica.cpp
// A simple example on (0,1)^3, derived from Mujica_ComputMech_2004
// Graham Harper, James Liu, ColoState; 2014/07--2017/11

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

#ifndef PI
#define PI 3.141592653589793
#endif

#ifndef lambda
#define lambda  4.008376878002496e+05
#endif

#ifndef mu
#define mu  80.194000225029669
#endif

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

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

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


// Dirichlet boundary conditions

PtVec3d fxnuD(PtVec3d pt)
{
  //double x = pt.xCrd();
  //double y = pt.yCrd();
  //double z = pt.zCrd();
  
  PtVec3d uD(0,0,0);
  return uD;
}

PtVec3d fxntN(PtVec3d pt)
{
  double x = pt.xCrd();
  double y = pt.yCrd();
  double z = pt.zCrd();

  PtVec3d tN(0,0,0);
  if (x<=0.5 && y<=0.5 && z==1)  tN = PtVec3d(0,0,-1);

  return -tN;  // Note the negative sign
}


// Function for the known displacement

PtVec3d fxnu(PtVec3d pt)
{
  //double x = pt.xCrd();
  //double y = pt.yCrd();
  //double z = pt.zCrd();
  
  return PtVec3d(0,0,0);
}


// Function for the stress tensor sigma

Mat3 fxnsigma(PtVec3d pt)
{
  //double x = pt.xCrd();
  //double y = pt.yCrd();
  //double z = pt.zCrd();
  
  Mat3 sigma;
  return sigma;
}


// Function for the known displacement divergence

double fxndivu(PtVec3d pt)

{
  //double x = pt.xCrd();
  //double y = pt.yCrd();
  //double z = pt.zCrd();
  
  double divu = 0;
  return divu;
}

// EqnBC_LinElas3d_Mujica.cpp
