function EqnBIC = EqnBIC_LinPoroElas_LiuWangICCS2021_Ex1
%% EqnBIC: Linear poroelasticity: A simple example on [0,1]^2 
% James Liu, ColoState; Zhuoran Wang, SYSU; 2020/09--2021/04 

% These parameters may change 
lambda = 1e6;  mu = 1;  % Lame constants 
kappa = 1;  c0 = 0;  
alpha = 1;  
T = 1;  % Final time for simulation 

EqnBIC = struct('T',T, ... 
  'lambda',lambda, 'mu',mu, 'kappa',kappa, 'c0',c0, 'alpha',alpha, ... 
  'fxnf',@fxnf, 'fxnuD',@fxnuD, 'fxntN',@fxntN, ...
  'fxnK',@fxnK, 'fxns',@fxns, 'fxnpD',@fxnpD, 'fxnuN',@fxnuN, ... 
  'fxnu',@fxnu, 'fxnp',@fxnp, 'fxndivu',@fxndivu, 'fxngradu',@fxngradu, ...
  'fxnv', @fxnv,'fxnsigma',@fxnsigma); 

% For solid: Body force: The right-hand side of the elasticity eqn. 
function f = fxnf(pt,t)
  x = pt(:,1);  y = pt(:,2);
  f1 =  (1 - 4*(sin(pi*x)).^2) .* sin(2*pi*y);
  f2 = -(1 - 4*(sin(pi*y)).^2) .* sin(2*pi*x);
  f1 = 2*mu*(pi^2) * f1;
  f2 = 2*mu*(pi^2) * f2;
  f1 = f1 - 2*mu/(1+lambda) * pi^2 * sin(pi*x).*sin(pi*y);
  f2 = f2 - 2*mu/(1+lambda) * pi^2 * sin(pi*x).*sin(pi*y);
  f1 = f1 + (lambda+mu-alpha)/(1+lambda) * (pi^2) * cos(pi*(x+y));
  f2 = f2 + (lambda+mu-alpha)/(1+lambda) * (pi^2) * cos(pi*(x+y));
  f = -sin(pi/2*t) * [f1,f2];
end

% For solid (displacement): Dirichlet boundary conditons 
function uD = fxnuD(pt,t)
  uD = fxnu(pt,t);
end

% For solid (displacement): Neumann boundary conditons (traction) 
% function tN = fxntN(pt,t) 
% end 

% For fluid: Permeability tensor 
function K = fxnK(pt)
  NumPts = size(pt,1);
  K = zeros(NumPts,2,2);
  K(:,1,1) = kappa;  
  K(:,2,2) = kappa;  
end

% For fluid source: The source for the Darcy equation 
function s = fxns(pt,t) 
  x = pt(:,1);  y = pt(:,2);
  s = (alpha+c0) * cos(pi/2*t)*(pi/2) * pi/(1+lambda) * sin(pi*(x+y));
  s = s + sin(pi/2*t) * kappa*(2*pi^2)* pi/(1+lambda) * sin(pi*(x+y));
end

% For fluid (pressure): Dirichlet boundary conditons 
function pD = fxnpD(pt,t) 
  pD = fxnp(pt,t);
  % x = pt(:,1);  y = pt(:,2);
  % pD = sin(pi/2*t) * pi/(1+lambda) * sin(pi*(x+y));
end 

% For fluid (pressure): Neumann boundary conditons (flux) 
% function uN = fxnuN(pt,t) 
% end 

% For solid displacement: Known analytical solution 
function u = fxnu(pt,t)
  x = pt(:,1);  y = pt(:,2);
  u1 =  (sin(pi*x)).^2 .* sin(2*pi*y);
  u2 = -(sin(pi*y)).^2 .* sin(2*pi*x);
  u1 = u1 + 1/(1+lambda) * sin(pi*x).*sin(pi*y);
  u2 = u2 + 1/(1+lambda) * sin(pi*x).*sin(pi*y);
  u = sin(pi/2*t) * [u1,u2];
end

% For fluid pressure: Known analytical solution 
function p = fxnp(pt,t)
  x = pt(:,1);  y = pt(:,2);
  p = sin(pi/2*t) * pi/(1+lambda) * sin(pi*(x+y));
end

% For solid dilation: Known analytical solution 
function divu = fxndivu(pt,t)
  x = pt(:,1);  y = pt(:,2);
  divu = sin(pi/2*t) * pi/(1+lambda) * sin(pi*(x+y));
end

% For solid displacement gradient: Known analytical solution 
function gradu = fxngradu(pt,t)
  x = pt(:,1);  y = pt(:,2);
  NumPts = size(pt,1);
  gradu = zeros(NumPts,2,2);
  gradu(:,1,1) = sin(pi/2*t) * ( pi*sin(2*pi*x).*sin(2*pi*y) + pi/(1+lambda)*cos(pi*x).*sin(pi*y));
  gradu(:,1,2) = sin(pi/2*t) * ( (sin(pi*x)).^2.*cos(2*pi*y)*(2*pi) +  pi/(1+lambda)*sin(pi*x).*cos(pi*y));
  gradu(:,2,1) = sin(pi/2*t) * (-(sin(pi*y)).^2.*cos(2*pi*x)*(2*pi) +  pi/(1+lambda)*cos(pi*x).*sin(pi*y));
  gradu(:,2,2) = sin(pi/2*t) * (-pi*sin(2*pi*y).*sin(2*pi*x) + pi/(1+lambda)*sin(pi*x).*cos(pi*y));
end

% For fluid velocity: Known analytical solution   
function v = fxnv(pt,t)
  x = pt(:,1); y = pt(:,2);
  v1 = sin(pi/2*t) * (-kappa) * pi/(1+lambda) * pi*cos(pi*(x+y));
  v2 = sin(pi/2*t) * (-kappa) * pi/(1+lambda) * pi*cos(pi*(x+y));
  v = [v1,v2];
end

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

end