function EqnBC = EqnBC_Darcy_Ex13HetePerm
%% EqnBC: Darcy: Heterogeneous permeability 
% Example on the unit square with a heterogeneous permeability 
% James Liu, Farrah Sadre-Marandi, ColoState; 2012/07--2018/12 

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

% Permeability as a 2x2 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

% Permeability as a scalar 
function Ks = fxnKs(pt)
  % Loading 20x20 matrix PERMK 
  % load Darcy_Ex13HetePerm.mat 
  % Temporarily 
  PERMK = ones(20,20);
  PERMK(1,[1,4,5,10,11,16,18]) = 1e-6;
  PERMK(2,[5,9,10,11]) = 1e-6;
  PERMK(3,[6]) = 1e-6;
  PERMK(4,[10,16,20]) = 1e-6;
  PERMK(5,[7,12,19]) = 1e-6;
  PERMK(6,[4,8,14]) = 1e-6;
  PERMK(7,[9,14,15,18]) = 1e-6;
  PERMK(8,[18,19]) = 1e-6;
  PERMK(9,[1,3,8,9,10,13]) = 1e-6;
  PERMK(10,[10,12,18]) = 1e-6;
  PERMK(11,[1,11,16]) = 1e-6;
  PERMK(12,[5,7,8,10,20]) = 1e-6;
  PERMK(13,[3,6,7,8,11,15,19]) = 1e-6;
  PERMK(14,[1,3,6,8,9]) = 1e-6;
  PERMK(15,[10,16,17]) = 1e-6;
  PERMK(16,[2,3,6,7,12,19]) = 1e-6;
  PERMK(17,[2,3,8,15,19]) = 1e-6;
  PERMK(18,[5,19]) = 1e-6;
  PERMK(19,[10,11,14,18,20]) = 1e-6;
  PERMK(20,[12,14,16]) = 1e-6;
  % 
  x = pt(:,1);  y = pt(:,2);
  I = 20 - floor(20*y);
  J = floor(20*x) + 1;
  NumPts = size(pt,1);
  Ks = zeros(NumPts,1);
  for k=1:NumPts
    Ks(k) = PERMK(I(k),J(k));
  end
end

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

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

% Neumann boundary conditions
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