function EqnBIC = EqnBIC_LinPoroElas_Mandel
%% EqnBIC: Linear Poroelasticity: 
% Most parameters are allowed to change 
% Zhuoran Wang, SYSU; James Liu, ColoState; 2017/07--2021/04 

%% For elasticity 
% nu = (3*K-2*mu)/(2*(3*K+mu));  % Poisson's coefficient
% E = 9*K/(3*K+mu)*mu;  % Young's modulus 
E = 1e4;  nu = 0.4;  lambda = E*nu/(1-2*nu)/(1+nu);  mu = E/2/(1+nu); 

%% For Darcy 
kappa = 1e-6;  % Choices: 1e-6, 1e-3, 1, 1e3 
c0 = 0;  

%% For ... 
% K = lambda + 2/3*mu;  % Skeleton bulk modulus
% K = 1e-06;
% Ku = K+alpha^2/c0; % undrained version
% cf = kappa/c0*(K+4/3*mu)/(Ku+4/3*mu); % Fluid diffusivity coeff.
cf = kappa*(lambda + 2*mu);  
% B = alpha/(c0*Ku); % Skempton's coeff.
% B = 3/(alpha*(1-2*nu))*(nuu-nu)/(1+2*nuu);
% nuu = ((B-2*B*nu)/3+nu)/(1-2*B*(1-2*nu)/3);
B = 1;  nuu = (3*nu+B*(1-2*nu))/(3-B*(1-2*nu));

%% For coupling 
alpha = 1; 

%% More parameters 
a = 1; b = 1;  % Domain size 
F = 1;  % External force 
T = 50;  % Final time 

EqnBIC = struct('T',T, ...
  'lambda',lambda, 'mu',mu, 'alpha',alpha, 'c0',c0,'nuu',nuu,'kappa',kappa, ...
  'nu',nu,'E',E,'cf',cf,'B',B,'F',F,'a',a,'b',b, ...
  'fxnf',@fxnf, 'fxnuDpartial',@fxnuDpartial, 'fxntN',@fxntN, ...
  'fxnK',@fxnK, 'fxns',@fxns, 'fxnpD',@fxnpD,'fxnuN',@fxnuN, ...
  'fxnu',@fxnu, 'fxnuGrad',@fxnuGrad,'fxnp',@fxnp, 'fxnsigma',@fxnsigma);

%% Import values of alphan
terms = 1000;
[alphan,An,Bn] = Mandel_Newton;

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

%% For solid: Displacement: partial Dirichlet boundary conditons 
function uDp = fxnuDpartial(pt,t,type)
  x = pt(:,1);  y = pt(:,2);
  if type == 3 & x == 0
      uDp(:,1) = 0;
  end
  if type == 4 & y == 0
      uDp(:,2) = 0; 
  end
  if type == 5 
     x = pt(:,1); y = pt(:,2);
     NumPts = size(pt,1);
     for i=1:NumPts  % modify
       xhat = x(i)/a; that = cf*t/(a^2);
       for j = 1:terms
             u2(i,j) = An(j)*exp(-alphan(j)^2*that)*cos(alphan(j));
       end
     end
  uDp(:,2) =  (-F*(1-nu)/(2*mu*a)+F*(1-nuu)/(mu*a)*sum(u2,2)).*y;
  end
end

%% For solid: Displacement: Neumann boundary conditons (traction)
% this is for test mandel's problem
function tN = fxntN(pt,t)
  x = pt(:,1);  y = pt(:,2);
  bndry1 = (y==0);  bndry2 = (x==1);  bndry3 = (y==1);  bndry4 = (x==0);
  tN1 = 0*bndry2;
  tN2 = 0*bndry2;
  tN = [tN1, tN2];  
end

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

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

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

%% For solid displacement: 
function uhat = fxnu(pt,t)
  x = pt(:,1);  y = pt(:,2);
  NumPts = size(pt,1);
  for i=1:NumPts  % modify
    xhat = x(i)/a; that = cf*t/(a^2);
    for j = 1:terms
        u1_1(i,j) = An(j)*exp(-alphan(j)^2*that)*cos(alphan(j));
        u1_2(i,j) = Bn(j)*exp(-alphan(j)^2*that)*sin(alphan(j)*xhat);
        u2(i,j) = An(j)*exp(-alphan(j)^2*that)*cos(alphan(j));
    end
  end
  u1hat = nu/(2*mu).*x - nuu/mu*sum(u1_1,2).*x + 1/mu*sum(u1_2,2);
  u2hat = -(1-nu)/(2*mu).*y + (1-nuu)/mu*sum(u2,2).*y;
  uhat = [u1hat*F/a,u2hat*F/a];  % original
end

%% For solid displacement gradient: 
function uGrad = fxnuGrad(pt,t)
  x = pt(:,1);  y = pt(:,2);
  NumPts = size(pt,1);
  for i=1:NumPts  % modify
    xhat = x(i)/a; that = cf*t/(a^2);
    for j = 1:terms
        u1_1(i,j) = An(j)*exp(-alphan(j)^2*that)*cos(alphan(j));
        u1_2(i,j) = Bn(j)*exp(-alphan(j)^2*that)*alphan(j)/a*cos(alphan(j)*xhat);
        u2(i,j) = An(j)*exp(-alphan(j)^2*that)*cos(alphan(j));
    end
  end
  u1hat_gradx = nu/(2*mu) - nuu/mu*sum(u1_1,2) + 1/mu*sum(u1_2,2);
  u2hat_grady = -(1-nu)/(2*mu) + (1-nuu)/mu*sum(u2,2);
  u1hat_grady = 0;
  u2hat_gradx = 0;
  uGrad(:,1,1) = u1hat_gradx * F/a;
  uGrad(:,2,2) = u2hat_grady * F/a;
  uGrad(:,1,2) = 0;
  uGrad(:,2,1) = 0;
end

%% For fluid pressure: Series for the exact solution 
function phat = fxnp(pt,t)
  x = pt(:,1);  y = pt(:,2);
  NumPts = size(pt,1);
  for i=1:NumPts  % modify
    xhat = x(i)/a; that = cf*t/(a^2);
    for j = 1:terms
        px(i,j) = An(j)*exp(-alphan(j)^2*that)* ...
                 (cos(alphan(j)*xhat) - cos(alphan(j)));
    end
  end
  phat = 2*sum(px,2)*F/a*B*(1+nuu)/3;  % Original
  % phat = 2*sum(px,2);  % Normalized
end

% % This is the normalized version 
% function sigmahat = fxnsigma(pt,t) 
%   x = pt(:,1);  y = pt(:,2);
%   NumPts = size(pt,1);
%   for i=1:NumPts  % modify
%     xhat = x(i)/a; that = cf*t/(a^2);
%     for j = 1:terms
%         sigma(i,j) = An(j)*exp(-alphan(j)^2*that)*(cos(alphan(j)) ...
%                    - B*(nuu-nu)/(1-nu)*cos(alphan(j)*xhat));
%     end
%   end
%   sigmahat = zeros(NumPts,2,2);
%   sigmahat(:,1,1) = 0;
%   sigmahat(:,1,2) = 0;
%   sigmahat(:,2,1) = 0;
%   sigmahat(:,2,2) = -1+2*sum(sigma,2);
% end

end