function EqnBC = EqnBC_LinElas_CarsScheIMAJNA2015_Ex1
%% EqnBC: Linear Elasticity: CarsSche_IMAJNA_2015, Ex.1, p.1616 
% Example on the unit square 
% James Liu, ColoState; 2012/07--2017/02

% Parameters (Changing nu to test locking)
E = 1e5;  % Elastic modulus 
nu = 0.499999;  % Poisson ratio: Choices: 0.49, 0.4999, 0.499999 
lambda = E*nu/((1-2*nu)*(1+nu));
mu = E/(2*(1+nu));
omega = (3/4)*pi; 
a = 0.544483737;  % Solution of sin(2*omega)*a + sin(2*omega*a) = 0 
C1 = -cos((1+a)*omega)/cos((1-a)*omega); 
C2 = 2*(lambda+2*mu)/(lambda+mu);

EqnBC = struct('lambda',lambda, 'mu',mu, ...
  'fxnf',@fxnf, 'fxnuD',@fxnuD, ...
  'fxnu',@fxnu, 'fxndivu',@fxndivu);

% The right-hand side vector function in the linear elasticity equation
function f = fxnf(pt)
  x = pt(:,1);  y = pt(:,2);
  f = [-mu*(pi^3)*sin(2*pi*y).*(2*cos(2*pi*x)-1), ...
        mu*(pi^3)*sin(2*pi*x).*(2*cos(2*pi*y)-1)];
end

% Dirichlet boundary condition
function uD = fxnuD(pt) 
  % x = pt(:,1);  y = pt(:,2);
  uD = fxnu(pt);
end

% Neumann boundary condition
% None 

% Known analytical solution for displacement/deformation 
function u = fxnu(pt)
  x = pt(:,1);  y = pt(:,2);
  u1 =  0.5*pi * (sin(pi*x)).^2 .* sin(2*pi*y);
  u2 = -0.5*pi * sin(2*pi*x) .* (sin(pi*y)).^2;
  u = [u1,u2];
end

% Divergence of displacement/deformation unknown or difficult to compute 
function divu = fxndivu(pt)
  x = pt(:,1);  % y = pt(:,2);
  divu = zeros(size(x));
end

end