function EqnBIC = EqnBIC_LinPoroElas_YiSINUM2017_Ex1
%% EqnBIC: Linear poroelasticity: An example on [0,1]^2 
% Adopted from Yi,SINUM (2017) p.1932 Ex.1 
% Zhuoran Wang, James Liu, ColoState; 2020/03--2020/12 

% These parameters may change 
lambda = 1e4;  % Choices: 1, 1e4, 1e8, etc. 
mu = 1;  % Lame constants 
kappa = 1;  c0 = 0;  
alpha = 1;  
T = 1;  % Final time for simulation 
a = 1;  b = 1;  c = 1;

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

% For solid: Body force: The right-hand side of the elasticity eqn. 
function f = fxnf(pt,t)
  x = pt(:,1);  y = pt(:,2);
end

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

% For solid: 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 right-hand side of the Darcy equation 
function s = fxns(pt,t) 
  x = pt(:,1); y = pt(:,2);
end

% For fluid: Dirichlet boundary conditons 
function pD = fxnpD(pt,t) 
  pD = fxnp(pt,t);
  % x = pt(:,1);  y = pt(:,2);
  % pD = exp(-t) * sin(pi*x).*sin(pi*y); 
end 

% For fluid: 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 = -2*(sin(pi*x)).^2 .* sin(2*pi*y) + sin(pi*x).*sin(pi*y)/(lambda+mu); 
  u2 =  2*(sin(pi*y)).^2 .* sin(2*pi*x) + sin(pi*x).*sin(pi*y)/(lambda+mu);
  u1 = exp(-t)*u1;  
  u2 = exp(-t)*u2; 
  u = [u1, u2];
end

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

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

end