function EqnBC = EqnBC_StokesDarcy_Harper3CMAME2021_Ex3 
%% Stokes dom.(0,2)\times(0,1); Darcy dom.(0,2)\times(-1,0); interface y=0 
% Lid-driven cavity + Heterogeneous hydraulic conductivity  
% BndryCondType = [0; 2; 2; 1; 1; 1; 2];  % [C; D; D; S; S; S; D]
% James Liu, ColoState; 2019/07--2021/04 

%% 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: Zero 
function f = fxnf_Stokes(pt)
  f = zeros(size(pt));
end

%% Stokes: Dirichlet boundary conditions 
function uD = fxnuD_Stokes(pt)
  x = pt(:,1);  y = pt(:,2);
  uD = zeros(size(pt));
  idx3 = find(abs(y-1)<1e-12);  % Top boundary y=1 (top) 
  uD(idx3,1) = 1;
end

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

%% Stokes: No known exact solution for velocity, so set to 0 
function u = fxnu_Stokes(pt)
  u = zeros(size(pt));
end

%% Stokes: No known exact solution for pressure yet , so set to 0 
function p = fxnp_Stokes(pt)
  x = pt(:,1);  % y = pt(:,2);
  p = zeros(size(x));
end

%% Darcy: Permeability as an elementwise constant scalar matrix 
function K = fxnK_Darcy(pt)
  kappa = 10^(-6);
  PERMK = ones(10,10);
  PERMK(7,[2,4,7,9]) = kappa;
  PERMK(8,[2,5]) = kappa;
  x = pt(:,1);  y = pt(:,2);
  I = 10 - floor(10*(y+1)/2);  % Note: y \in (-1,0) 
  J = floor(5*x) + 1;          % Note: x \in (0,2) 
  NumPts = size(pt,1);
  Ks = ones(NumPts,1);
  for k=1:NumPts
    Ks(k) = PERMK(I(k),J(k));
  end
  K = zeros(NumPts,2,2);
  K(:,1,1) = Ks;
  K(:,2,2) = Ks;
end

%% Darcy: Source term: Zero 
function s = fxns_Darcy(pt)
  x = pt(:,1);  % y = pt(:,2);
  s = zeros(size(x));
end

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

%% Darcy: Dirichlet boundary conditions for bottom (y=-1)
function pD = fxnpD_Darcy(pt)
  x = pt(:,1);  y = pt(:,2);
  pD = zeros(size(x));
  idx1 = find(abs(y+1)<1e-12);  % Bootom boundary y=-1 
  pD(idx1) = ones(size(idx1));
end

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

%% Darcy: No known exact solution for pressure, so set to 0 
function p = fxnp_Darcy(pt)
  x = pt(:,1);  % y = pt(:,2);
  p = zeros(size(x));
end

%% Darcy: No known exact solution for velocity, so set to 0 
function u = fxnu_Darcy(pt)
  % x = pt(:,1);  y = pt(:,2);  
  u = zeros(size(pt));
end

end