function EqnBIC = EqnBIC_LinPoroElas_KansRivieJSC2018_ExVar
%% EqnBIC: Linear Poroelasticity: 
% KansRivie,JSC,2018, Ex. variant: An example on [0,1]^2  
% James Liu, ColoState; Zhuoran Wang, SYSU; 2017/07--2021/02 

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

%% Auxiliary function: phi(x,y) and its derivatives 
function z = phi(x,y) 
  z = sin(2*pi*x) .* sin(2*pi*y);
end
% function dx = dxphi(x,y) 
%   dx = (2*pi) * cos(2*pi*x) .* sin(2*pi*y);
% end  
% function dy = dyphi(x,y) 
%   dy = sin(2*pi*x) * (2*pi) .* cos(2*pi*y);
% end 
function z = psi(t) 
  z = (8*pi^2)*sin(2*pi*t) - (2*pi)*cos(2*pi*t) + (2*pi)*exp(-8*pi^2*t);
  z = z/(64*pi^4+4*pi^2);
end
function z = psi1(t)
  z = (16*pi^3)*cos(2*pi*t) + (4*pi^2)*sin(2*pi*t) ... 
    - (16*pi^3)*exp(-8*pi^2*t);
  z = z/(64*pi^4+4*pi^2);
end

%% For linear poroelasticity 
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); 

% For solid: Body force: The right-hand side in the elasticity PDE 
function f = fxnf(pt,t)
  x = pt(:,1);  y = pt(:,2);
  f = zeros(size(pt));
  f(:,1) = psi(t)*(2*pi) * cos(2*pi*x).*sin(2*pi*y) * (2*mu+lambda);
  f(:,2) = psi(t)*(2*pi) * cos(2*pi*y).*sin(2*pi*x) * (2*mu+lambda);
  f(:,1) = f(:,1) + alpha*psi(t) * (2*pi) * cos(2*pi*x).*sin(2*pi*y);
  f(:,2) = f(:,2) + alpha*psi(t) * (2*pi) * cos(2*pi*y).*sin(2*pi*x);
end

% For solid: Displacement: Dirichlet boundary conditons 
function uD = fxnuD(pt,t)
  x = pt(:,1);  y = pt(:,2);
  bndry1 = (y==0);  bndry2 = (x==1);  bndry4 = (x==0);
  u1D = 0*bndry1 + 0*bndry2 + 0*bndry4;
  u2D = 0*bndry1 + 0*bndry2 + 0*bndry4;
  uD = [u1D, u2D];
end

% For solid: Displacement: Neumann boundary conditons (traction)
function tN = fxntN(pt,t)
  x = pt(:,1);  y = pt(:,2);
  bndry3 = (y==1);
  tN = -[0*bndry3, -1*bndry3];  % Note: The first negative sign 
end

% For fluid: Permeability tensor in the Darcy eqn. 
function K = fxnK(pt)
  % x = pt(:,1);  y = pt(:,2);
  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 = (-psi1(t)+(8*pi^2)*psi(t)) * sin(2*pi*x).*sin(2*pi*y);
end

% For fluid: Pressure: Dirichlet boundary conditons 
function pD = fxnpD(pt,t)
  x = pt(:,1);  y = pt(:,2);
  bndry3 = (y==1);
  pD = 0*bndry3;
end

% For fluid: Pressure: Neumann boundary conditons 
function uN = fxnuN(pt,t)
  x = pt(:,1);  y = pt(:,2);
  bndry1 = (y==0);  bndry2 = (x==1);  bndry4 = (x==0);
  uN = 0*bndry1 + 0*bndry2 + 0*bndry4;
end

% For solid: Displacement: No known analytical solution 
function u = fxnu(pt,t)
  x = pt(:,1);  y = pt(:,2);
  u = zeros(size(pt));
  u(:,1) = psi(t)/(4*pi) * cos(2*pi*x).*sin(2*pi*y);
  u(:,2) = psi(t)/(4*pi) * cos(2*pi*y).*sin(2*pi*x);
end

% For fluid: Pressure: Known analytical solution   
function p = fxnp(pt,t)
  x = pt(:,1);  y = pt(:,2);
  p = psi(t) * phi(x,y);
end

% For solid displacement gradient: Known analytical solution 
function gradu = fxngradu(pt,t)
  NumPts = size(pt,1);
  gradu = zeros(NumPts,2,2);
  % x = pt(:,1);  y = pt(:,2);
end

end