function EqnBC = EqnBC_Poisson_MuWangYeIJNAM2015_Case2
%% EqnBC: Poisson: Exact sln. p(x,y)=x(1-x)y(1-y) w/ degenerate diffusion  
% Example on (0,1)^2 with homogeneous Dirichlet and/or Neumann bndry.cond. 
% James Liu, ColoState; 2012/07--2018/07 

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

% Diffusion coeff. or permeability as a scalar 
function Ks = fxnKs(pt)
  x = pt(:,1);  y = pt(:,2);
  Ks = x.*y;
end

% Diffusion coeff. or permeability as 2x2 SPD matrix 
function K = fxnK(pt)
  NumPts = size(pt,1);
  K = zeros(NumPts,2,2);
  x = pt(:,1);  y = pt(:,2);
  K(:,1,1) = x.*y;
  K(:,2,2) = x.*y;
end

% The right-hand sied function in the Poisson/Darcy equation
function f = fxnf(pt)
  x = pt(:,1);  y = pt(:,2);
  f = (1-4*x).*(y.^2).*(1-y) + (x.^2).*(1-x).*(1-4*y);
  f = -f;
end

% Dirichlet boundary condition
function  pD = fxnpD(pt) 
  x = pt(:,1);  % y = pt(:,2);
  pD = zeros(size(x));
end

% Neumann boundary condition
EPSILON = 1e-12;
function uN = fxnuN(pt)
  x = pt(:,1);  y = pt(:,2);
  uN = (abs(x-0)<EPSILON).*( x.*(1-2*x).*(y.^2).*(1-y)) ...  % Actually zero 
     + (abs(x-1)<EPSILON).*(-x.*(1-2*x).*(y.^2).*(1-y)) ...
     + (abs(y-0)<EPSILON).*( (x.^2).*(1-x).*y.*(1-2*y)) ...  % Actually zero 
     + (abs(y-1)<EPSILON).*(-(x.^2).*(1-x).*y.*(1-2*y));
end

% Known exact solution for "pressure" 
function p = fxnp(pt)
  x = pt(:,1);  y = pt(:,2);
  p = x.*(1-x).*y.*(1-y);
end

% Known exact solution for "pressure" gradient
function pg = fxnpg(pt) 
  x = pt(:,1);  y = pt(:,2);
  pg = [(1-2*x).*y.*(1-y), x.*(1-x).*(1-2*y)];
end

% Known exact solution for "velocity"
function u = fxnu(pt)
  x = pt(:,1);  y = pt(:,2);
  u1 = -x.*(1-2*x).*(y.^2).*(1-y);
  u2 = -(x.^2).*(1-x).*y.*(1-2*y);
  u = [u1,u2];
end

end