function EqnBC = EqnBC_Darcy_Ex11ThreeChannels
%% EqnBC: Darcy: Fast and slow channels combined 
% Example on [0,1]^2 with 
%   Chanel#1: [0.2,0.8]x[0.8,0.9] w/ permeability= 
%   Chanel#2: [0.2,0.8]x[0.4,0.6] w/ permeability= 
%   Chanel#3: [0.2,0.8]x[0.1,0.2] w/ permeability= 
%   The rest w/ permeability=1 
% James Liu, ColoState; 2012/07--2017/02 

EqnBC = struct('fxnKs',@fxnKs, 'fxnK',@fxnK, 'fxnf',@fxnf, ...
  'fxnpD',@fxnpD, 'fxnuN',@fxnuN, ...
  'fxnp',@fxnp, 'fxnpg',@fxnpg, 'fxnu',@fxnu);

% Permeability as a scalar 
function Ks = fxnKs(pt)
  x = pt(:,1);  y = pt(:,2);
  Ks = (1e-3)*ones(size(x));
  % For Channel#1: 
  k1 = find(0.2<=x & x<=0.8 & 0.8<=y & y<=0.9);
  Ks(k1) = 2;
  % For Channel#2: 
  k2 = find(0.2<=x & x<=0.8 & 0.4<=y & y<=0.4);
  Ks(k2) = 1;
  % For Channel#3: 
  k3 = find(0.2<=x & x<=0.8 & 0.1<=y & y<=0.2);
  Ks(k3) = 2;
end

% Permeability as a 2x2 full SPD matrix 
function K = fxnK(pt)
  Ks = fxnKs(pt);
  NumPts = size(pt,1);
  K = zeros(NumPts,2,2);
  K(:,1,1) = Ks(:);
  K(:,2,2) = Ks(:);
end

% Source (the right-hand function)
function f = fxnf(pt)
  x = pt(:,1);  % y = pt(:,2);
  f = zeros(size(x));
end

% Dirichlet boundary condition
function pD = fxnpD(pt) 
  x = pt(:,1);  % y = pt(:,2);
  bnd0 = (x==0);
  bnd1 = (x==1);
  pD = 1*bnd0 + 0*bnd1;
end

% Neumann boundary condition
function uN = fxnuN(pt)
  x = pt(:,1);  y = pt(:,2);
  bnd = (y==0 | y==1);
  uN = 0*bnd;
end

% Exact pressure UNKNOWN or just set to zero
function p = fxnp(pt)
  x = pt(:,1);  % y = pt(:,2);
  p = zeros(size(x));
end

% Exact pressure gradient UNKNOWN or just set to zero
function pg = fxnpg(pt) 
  % x = pt(:,1);  y = pt(:,2);
  pg = zeros(size(pt));
end

% Exact velocity UNKNOWN or just set to zero
function u = fxnu(pt) 
  u = zeros(size(pt));
end

end