function EqnBIC = EqnBIC_LinPoroElas_CantileverBracket
%% EqnBIC: Linear Poroelasticity: 
% Example on (0,1)^2 
% James Liu, ColoState; 2012/07--2018/01

% These parameters shoul be allowed to change 
E = 1e5;  nu = 0.5 - 1e-1;
lambda = E*nu/((1-2*nu)*(1+nu));  mu = E/(2*(1+nu));  % Lame constants 
alpha = 0.93;  c0 = 0.0;  
kappa = 1e-7;  % kappa = 1e-7;
T = 5e-3;  % ? Final time for simulation, Deltat = 1e-3 

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

% For solid: Body force: The right-hand side in the 1st PDE 
function f = fxnf(pt,t)
  % x = pt(:,1);  y = pt(:,2);
  f = zeros(size(pt));
end

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

% For solid: Displacement: Neumann boundary conditons (traction)
function tN = fxntN(pt,t)
  x = pt(:,1);  y = pt(:,2);
  bndry1 = (y==0);  bndry2 = (x==1);  bndry3 = (y==1);  bndry4 = (x==0);
%   t1N =    0*bndry3;
%   t2N = (-1)*bndry3;
  tN1 = 0*bndry1 + 0*bndry2 +    0*bndry3 + 0*bndry4;
  tN2 = 0*bndry1 + 0*bndry2 + (-1)*bndry3 + 0*bndry4;
  tN = -[tN1, tN2];  % Note: The "-" negative sign 
end

% For fluid: Permeability: A function in the 2nd PDE 
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 right-hand side in the 2nd PDE  
function s = fxns(pt,t) 
  x = pt(:,1);  % y = pt(:,2);
  s = zeros(size(x));
end

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

% For fluid: Pressure: Neumann boundary conditons 
function uN = fxnuN(pt,t)
  x = pt(:,1);  y = pt(:,2);
  bndry1 = (y==0);  bndry2 = (x==1);  bndry3 = (y==1);  bndry4 = (x==0);
  uN = 0*bndry1 + 0*bndry2 + 0*bndry3 + 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));
end

% For fluid: Pressure: No known analytical solution   
function p = fxnp(pt,t)
  x = pt(:,1);  % y = pt(:,2);
  p = zeros(size(x));
end

end