function EqnBC = EqnBC_StokesDarcy_DiscaQuarteRevMatComplut2009_ExP383Var 
%% Stokes dom.(0,1)\times(1,2); Darcy dom.(0,1)\times(0,1); Interface y=1 
% DiscaQuarte_RevMatComplut_2009, Ex. on P.383: Known analytical solutions 
% Slight modification in Stokes velocity for nonzero BJS condition 
% James Liu, Graham Harper, ColoState; 2019/07--2020/08 

%% Parameters 
nu = 1;  % For Stokes 
alpha = 1;  % For BJS interface condition 

EqnBC = struct('nu',nu, 'alpha',alpha, ... 
  'fxnf_Stokes',@fxnf_Stokes, ...
  'fxnuD_Stokes',@fxnuD_Stokes, 'fxntN_Stokes',@fxntN_Stokes, ...
  'fxnu_Stokes',@fxnu_Stokes, 'fxnp_Stokes',@fxnp_Stokes, ... 
  'fxnK_Darcy',@fxnK_Darcy, 'fxns_Darcy',@fxns_Darcy, ... 
  'fxng_Darcy',@fxng_Darcy, ...
  'fxnpD_Darcy',@fxnpD_Darcy, 'fxnuN_Darcy',@fxnuN_Darcy, ...
  'fxnp_Darcy',@fxnp_Darcy, 'fxnu_Darcy',@fxnu_Darcy);

%% Stokes: Body force: 
function f = fxnf_Stokes(pt)
  x = pt(:,1);  y = pt(:,2);
  f1 =  pi^2/2 * sin(pi/2*x) .* cos(pi/2*y) + 1;
  f2 = -pi^2/2 * cos(pi/2*x) .* sin(pi/2*y);
  f = -[f1,f2];
end

%% Stokes: Dirichlet boundary conditions 
function uD = fxnuD_Stokes(pt)
  % uD = fxnu_Stokes(pt);
  x = pt(:,1);  y = pt(:,2);  
  u1 = -sin(pi/2*x) .* cos(pi/2*y) + 1;  % JL20200728 modification: "+1" 
  u2 =  cos(pi/2*x) .* sin(pi/2*y) - (1-x);
  uD = [u1,u2];
end

%% JL20200810: To be revised for a general treatment 
%% Stokes: Neumann or traction boundary conditions 
function tN = fxntN_Stokes(pt)
  x = pt(:,1);  y = pt(:,2);
  NumPts = size(pt,1);
  sigma = zeros(NumPts,2,2);
  sigma(:,1,1) = 2*nu*(-pi/2) * cos(pi/2*x) .* cos(pi/2*y) - (1-x);
  sigma(:,2,2) = 2*nu*( pi/2) * cos(pi/2*x) .* cos(pi/2*y) - (1-x);
  sigma(:,1,2) =   nu; 
  sigma(:,2,1) =   nu; 
  % Assuming (y==2) the top-side is Neumann 
  tN1 = sigma(:,1,2);  tN2 = sigma(:,2,2);
  % % Assuming (x==0) the left-side is Neumann 
  % tN1 = -sigma(:,1,1);  tN2 = -sigma(:,2,1);
  % % Assuming (x==1) the right-side is Neumann 
  % tN1 = sigma(:,1,1);  tN2 = sigma(:,2,1);
  tN = [tN1, tN2];
end

%% Stokes: Known exact solution for velocity 
function u = fxnu_Stokes(pt)
  x = pt(:,1);  y = pt(:,2);  
  u1 = -sin(pi/2*x) .* cos(pi/2*y) + 1;  % JL20200728 modification: "+1" 
  u2 =  cos(pi/2*x) .* sin(pi/2*y) - (1-x);
  u = [u1,u2];
end

%% Stokes: Known exact solution for pressure 
function p = fxnp_Stokes(pt)
  x = pt(:,1);  % y = pt(:,2);
  p = 1 - x;
end

%% Darcy: Permeability as an elementwise constant 2x2 SPD matrix: [1,0;0,1] 
function K = fxnK_Darcy(pt)
  x = pt(:,1);  % y = pt(:,2);
  K(:,1,1) = ones(size(x));
  K(:,2,2) = ones(size(x));
  K(:,1,2) = zeros(size(x));
  K(:,2,1) = zeros(size(x));
end

%% Darcy: Source term 
function s = fxns_Darcy(pt)
  x = pt(:,1);  y = pt(:,2);
  s = pi * cos(pi/2*x) .* cos(pi/2*y);
end

%% Darcy: Additional term 
function g = fxng_Darcy(pt)
  g = zeros(size(pt));
end

%% Darcy: Dirichlet boundary condition for y=0 
function pD = fxnpD_Darcy(pt)
  % pD = fxnp_Darcy(pt);
  x = pt(:,1);  y = pt(:,2);
  pD = (2/pi)*cos(pi/2*x).*cos(pi/2*y) + (1-x).*y;
end

%% Darcy: Neumann boundary conditions: 
function uN = fxnuN_Darcy(pt)
  x = pt(:,1);  y = pt(:,2);
  u = fxnu_Darcy(pt);
  uN = (-1).*(x==0).*u(:,1) ...
     +       (x==1).*u(:,1) ... 
     + (-1).*(y==0).*u(:,2) ...
     +       (y==1).*u(:,2) ;  
end

%% Darcy: Known exact solution for pressure 
function p = fxnp_Darcy(pt)
  x = pt(:,1);  y = pt(:,2);
  p = (2/pi)*cos(pi/2*x).*cos(pi/2*y) + (1-x).*y;
end

%% Darcy: Known exact solution for velocity 
function u = fxnu_Darcy(pt)
  x = pt(:,1);  y = pt(:,2);  
  u1 = sin(pi/2*x).*cos(pi/2*y) + y;
  u2 = cos(pi/2*x).*sin(pi/2*y) - (1-x);  
  u = [u1, u2];
end

end