function EqnBC = EqnBC_Darcy_StanfordPerm
%% EqnBC: Darcy with Stanford Permeability on [0,5]x[0,1]
% James Liu, ColoState; Victor Ginting, UWYO; 2012/07--2018/12 

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

% Permeability as an elementwise 2x2 SPD matrix 
function K = fxnK(pt)
  tmp = load('Coeff100x100log.txt');
  PERMK = tmp(2:101, 2:101);
  PERMK = exp(2.0*PERMK);
  % Temporarily 
  x = pt(:,1);  y = pt(:,2);
  I = 100 - floor(100*y);
  J = floor(100*x/5) + 1;
  NumPts = size(pt,1);
  Kscalar = zeros(NumPts,1);
  for k=1:NumPts
    Kscalar(k) = PERMK(I(k), J(k));
  end
  K = zeros(NumPts,2,2);
  K(:,1,1) = Kscalar(:);
  K(:,2,2) = Kscalar(:);
end

% Permeability as an elementwise scalar 
function Ks = fxnKs(pt)
  tmp = load('Coeff100x100log.txt');
  PERMK = temp(2:121, 2:121);
  PERMK = exp(PERMK);
  % Temporarily 
  x = pt(:,1);  y = pt(:,2);
  I = 120 - floor(120*y);
  J = floor(120*x/5) + 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 in the Darcy equation)
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==5);
  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) 
  % x = pt(:,1);  y = pt(:,2);
  u = zeros(size(pt));
end

end