function EqnBC = EqnBC_StokesDarcy_ChenWangWangIMAJNA2016_Ex 
%% Stokes dom.(0,\pi)x(0,1); Darcy dom.(0,\pi)x(-1,0); Interface y=0 
% Adopted from Chen,Wang,Wang, IMAJNA, 2016, Ex. 
% By Graham Harper, James Liu, ColoState; 2014/07--2020/07 

%% 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: A bunch of auxiliary functions 
function v = v_Stokes(y) 
  v = (1/pi^2)*sin(pi*y).^2 - 2;
end
function v = v_p(y) 
  v = (1/pi)*sin(2*pi*y);
end
function v = v_pp(y) 
  v = 2*cos(2*pi*y);
end
function v = v_ppp(y) 
  v = -4*pi*sin(2*pi*y);
end

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

%% Stokes: Dirichlet boundary condition 
function uD = fxnuD_Stokes(pt)
  x = pt(:,1);  y = pt(:,2);  
  u1 = cos(x) .* ((1/pi)  *sin(2*pi*y));
  u2 = sin(x) .* ((1/pi^2)*sin(pi*y).^2 - 2);
  uD = [u1, u2];
end

%% Stokes: Neumann or traction boundary condition 
function tN = fxntN_Stokes(pt)
  x = pt(:,1);  % y = pt(:,2);
  tN1 = zeros(size(x));
  tN2 = zeros(size(x));
  tN = [tN1, tN2];  % GH Remarks: TO BE COMPUTED
end

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

%% Stokes: Known exact solution for pressure
function p = fxnp_Stokes(pt)
  x = pt(:,1);  y = pt(:,2);
  p = sin(x) .* sin(y);
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 = zeros(size(x));
end

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

%% Darcy: Dirichlet boundary condition
function pD = fxnpD_Darcy(pt)
  x = pt(:,1);  y = pt(:,2);  
  pD = sin(x) .* (exp(y) - exp(-y));
end

%% Darcy: Neumann boundary condition
function uN = fxnuN_Darcy(pt)
  x = pt(:,1);  % y = pt(:,2);
  uN = zeros(size(x));  % GH Remarks: TO BE COMPUTED
end

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

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

end