function EqnBC = EqnBC_LinElas_ArnoldAwanouQiuAdvComputMath2015_Ex
%% A simple example on [0,1]^2 for linear elasticity 
% See Arnold,Awanou,Qiu, Adv.Comput.Math, 2015, p.568 
% James Liu, Zhuoran Wang, ColoState; 2012/07--2017/02 

% Parameters 
lambda = 1.23e2;  % Original value: lambda = 123;
mu = 79.3;

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

% The right-hand side vector function in the linear elasticity equation
function f = fxnf(pt)
  x = pt(:,1);  y = pt(:,2);
  u1 = cos(pi*x) .* sin(2*pi*y);
  u2 = sin(pi*x) .* cos(pi*y);
  u1dxdx =   -(pi*pi)*u1;
  u1dydy = -(4*pi*pi)*u1;
  u1dxdy = -(2*pi*pi)*sin(pi*x).*cos(2*pi*y);
  u2dxdx =   -(pi*pi)*u2;
  u2dydy =   -(pi*pi)*u2;
  u2dydx =   -(pi*pi)*cos(pi*x).*sin(pi*y);
  f1 = -(lambda+mu)*(u1dxdx + u2dydx) - mu*(u1dxdx+u1dydy);
  f2 = -(lambda+mu)*(u1dxdy + u2dydy) - mu*(u2dxdx+u2dydy);
  f = [f1,f2];
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 
function u = fxnu(pt)
  x = pt(:,1);  y = pt(:,2);
  u1 = cos(pi*x) .* sin(2*pi*y);
  u2 = sin(pi*x) .* cos(pi*y);
  u = [u1,u2];
end

% Known analytical solution for divergence of displacement
function divu = fxndivu(pt)
  x = pt(:,1);  y = pt(:,2);
  u1dx = -pi*sin(pi*x).*sin(2*pi*y);
  u2dy = -pi*sin(pi*x).*sin(pi*y);
  divu = u1dx + u2dy;
end

% JL20160915: NEED DOUBEL-CHECK (Finished by Sophy) 
% 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) = - lambda*pi*sin(pi*x).*sin(pi*y) ...
                 - lambda*pi*sin(pi*x).*sin(2*pi*y) ...
                 - mu*2*pi*sin(pi*x).*sin(2*pi*y);
  sigma(:,1,2) = mu*pi*cos(pi*x).*cos(pi*y) + mu*2*pi*cos(pi*x).*cos(2*pi*y);
  sigma(:,2,1) = mu*pi*cos(pi*x).*cos(pi*y) + mu*2*pi*cos(pi*x).*cos(2*pi*y);
  sigma(:,2,2) = - lambda*pi*sin(pi*x).*sin(pi*y) ...
                 - lambda*pi*sin(pi*x).*sin(2*pi*y) ...
                 - mu*2*pi*sin(pi*x).*sin(pi*y);
end

end