function EqnBC = EqnBC_LinElas_BrenSINUM1993_Ex
%% A simple example for linear elasticity on [0,1]^2 
% See Brenner, SIAM J.Numer.Anal., 1993, p.134 
% James Liu, ColoState; 2012/07--2021/02 

% Parameters 
lambda = 1e6;
mu = 1;

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

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

% Known analytical solution for divergence of displacement
function divu = fxndivu(pt)
  x = pt(:,1);  y = pt(:,2);
  divu = pi/(1+lambda)*sin(pi*(x+y));
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);
  % The strain part 
  sigma(:,1,1) = -2*pi*sin(2*pi*x).*sin(2*pi*y) + pi/(1+lambda)*cos(pi*x).*sin(pi*y);
  sigma(:,2,2) =  2*pi*sin(2*pi*x).*sin(2*pi*y) + pi/(1+lambda)*sin(pi*x).*cos(pi*y);
  sigma(:,1,2) =  2*pi * (-sin(pi*x).*sin(pi*x).*cos(2*pi*y) ...
                        +sin(pi*y).*sin(pi*y).*cos(2*pi*x))... 
             + pi/(1+lambda)/2*sin(pi*(x+y));
  % Now the stress 
  sigma(:,1,1) = (2*mu)*sigma(:,1,1) + lambda*pi/(1+lambda)*sin(pi*(x+y));
  sigma(:,2,2) = (2*mu)*sigma(:,2,2) + lambda*pi/(1+lambda)*sin(pi*(x+y));
  sigma(:,1,2) = (2*mu)*sigma(:,1,2);
  % Finally symmetry 
  sigma(:,2,1) = sigma(:,1,2);
end

% JL20210206: TO BE FINISHED 
% 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);
end

end