function EqnBIC = EqnBIC_NlinPoroElas_ICCS2021Ex3
%% EqnBIC: Nonlinear poroelasticity: Example on (0,1)^2 for ICCS2021
% James Liu ColoState; Zhuoran Wang, SYSU; 2012/07--2021/01 

% These parameters are allowed to change 
lambda = 1;  mu = 1;  % Lame constants 
kappa = 1;  % Adjustable for testing  
c0 = 0;  
alpha = 1;  
T = 1;  % Final time for simulation 
a = 1;  b = 1;  % Size of domain
c = 0; % rigit-plate-constraint
parama = 0.1;  % parameter a for dilation-dependent permeability 

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

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

% For solid: Displacement: Dirichlet boundary conditons 
function uD = fxnuD(pt,t)
  x = pt(:,1);  y = pt(:,2); 
  u1 = sin(pi/2*t) * sin(pi*x).*sin(pi*y);
  u2 = sin(pi/2*t) * sin(pi*x).*sin(pi*y);
  uD = [u1,u2];
end

function uGrad = fxnuGrad(pt,t)
  x = pt(:,1);  y = pt(:,2);
  uGrad(:,1,1) = sin(pi/2*t) * pi*cos(pi*x).*sin(pi*y);
  uGrad(:,1,2) = sin(pi/2*t) * pi*sin(pi*x).*cos(pi*y);
  uGrad(:,2,1) = uGrad(:,1,1);
  uGrad(:,2,2) = uGrad(:,1,2);    
end

% For solid: Displacement: Neumann boundary conditons (traction) 

% For fluid: Permeability in Darcy eqn. 
function K = fxnK(pt,t)
  NumPts = size(pt,1);
  x = pt(:,1); y = pt(:,2);
  K = zeros(NumPts,2,2);
  dilation = sin(pi/2*t) * pi*sin(pi*(x+y)); 
  K(:,1,1) = (1+parama*dilation) * kappa;  
  K(:,2,2) = (1+parama*dilation) * kappa;  
end

% For fluid: Source: RHS of Darcy eqn. 
function s = fxns(pt,t) 
  x = pt(:,1); y = pt(:,2);
  nonlinear_term = 1+parama*sin(pi/2*t)...
                 * (pi*cos(pi*x).*sin(pi*y)...
     			  + pi*sin(pi*x).*cos(pi*y));
  s = pi/2*cos(pi/2*t)*(cos(pi*y)+1)*c0 ...
    + pi/2*cos(pi/2*t)*alpha ...
    *(pi*cos(pi*x).*sin(pi*y) ...
    + pi*sin(pi*x).*cos(pi*y)) ...
    + sin(pi/2*t)...
    *(nonlinear_term.*pi^2.*cos(pi*y)...
    + pi*sin(pi*y).*(sin(pi/2*t).*parama ...
    * (pi^2*cos(pi*x).*cos(pi*y) ...
    - pi^2*sin(pi*x).*sin(pi*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) * (cos(pi*y)+1);
end

% For solid displacement: Known analytical solution 
function u = fxnu(pt,t)
  x = pt(:,1);  y = pt(:,2);
  u1 = sin(pi/2*t) .* sin(pi*x).*sin(pi*y);
  u2 = sin(pi/2*t) .* sin(pi*x).*sin(pi*y);
  u = [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) * (cos(pi*y)+1);
end

% For fluid velocity: Known analytical solution   
function v = fxnv(pt,t)
  v = zeros(size(pt));
  x = pt(:,1); y = pt(:,2);
  nonlinear_term = 1+parama*sin(pi/2*t) ...
                  *(pi*cos(pi*x).*sin(pi*y) ...
         	      + pi*sin(pi*x).*cos(pi*y));
  v(:,2) = sin(pi/2*t)*nonlinear_term*pi.*sin(pi*y);
end

% ZW20210109: Need to edit stress 
% For solid: Known analytical solution 
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*cos(pi*x).*sin(pi*y);
  gradu(:,1,2) = sin(pi/2*t)*pi*sin(pi*x).*cos(pi*y);
  gradu(:,2,1) = gradu(:,1,1);
  gradu(:,2,2) = gradu(:,1,2);
  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*cos(pi*x).*sin(pi*y) + pi*sin(pi*x).*cos(pi*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