function EqnBC = EqnBC_LinElas_Wang3ZhangJCAM2016_Test3Locking
%% EqnBC: Lin. Elasticity: WangWangWangZhang_JCAM_2016, Test3: Locking-free
% On [0,1]^2, Dirichlet (and Neumann) bndry.conds. 
% James Liu, Zhuoran Wang, ColoState; 2012/07--2021/02 

% These two parameters (lambda,mu) can be adjusted  
lambda = 1e0;  % Choices: 1, 1e3, 1e6, ... test for lcoking
mu = 0.5;

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

% Body force or the RHS vector function in the linear elasticity equation
function f = fxnf(pt)
  x = pt(:,1);  y = pt(:,2);
  f = (2*mu) * [sin(x).*sin(y), cos(x).*cos(y)];
end

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

% Neumann boundary condition: The right side (x=1)
% NOTE: negative sign 
function tN = fxntN(pt)
  x = pt(:,1);  y = pt(:,2);
  tN = zeros(size(pt));
  if (x==1)
    tN(:,1) = -( 2*mu*(cos(x).*sin(y)) + 2*(mu/lambda+1) );
  end
end

% Known analytical solution for displacement
function u = fxnu(pt)
  x = pt(:,1);  y = pt(:,2);
  u = [sin(x).*sin(y), cos(x).*cos(y)] + (1/lambda)*[x,y];
end

% Known analytical solution for divergence of displacement (dilation) 
function divu = fxndivu(pt)
  % x = pt(:,1);  y = pt(:,2);
  divu = 2/lambda;
end

% Known analytical solution for stress 
function sigma = fxnsigma(pt) 
  x = pt(:,1);  y = pt(:,2);
  NumPts = size(pt,1);
  sigma = zeros(NumPts,2,2);
  sigma(:,1,1) =  2*mu*cos(x).*sin(y) + 2*(mu/lambda+1);
  sigma(:,2,2) = -2*mu*cos(x).*sin(y) + 2*(mu/lambda+1);
  % sigma(:,1,2) = sigma(:,2,1) = 0 
end

% Known analytical solution for displacement gradient  
function ug = fxnug(pt) 
  x = pt(:,1);  y = pt(:,2);
  NumPts = size(pt,1);
  ug = zeros(NumPts,2,2);
  ug(:,1,1) =  cos(x).*sin(y) + 1/lambda;
  ug(:,1,2) =  sin(x).*cos(y);
  ug(:,2,1) = -sin(x).*cos(y);
  ug(:,2,2) = -cos(x).*sin(y) + 1/lambda;
end

end