function EqnBC = EqnBC_Poisson_WheeYotovSINUM2006_Ex2
%% EqnBC: Poisson: Adopted from WheeYotov_SINUM_2006, Ex.2, p.2103
% James Liu, ColoState; 2012/07--2017/02

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

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

% The right-hand sied function in the Poisson/Darcy equation
function f = fxnf(pt)
  x = pt(:,1);  y = pt(:,2);
  px = 3*pi * sin(6*pi*x) .* sin(3*pi*y).^2;
  py = 3*pi * sin(3*pi*x).^2 .* sin(6*pi*y);
  pxx = 18*pi^2 * cos(6*pi*x) .* sin(3*pi*y).^2;
  pyy = 18*pi^2 * sin(3*pi*x).^2 .* cos(6*pi*y);
  pxy =  9*pi^2 * sin(6*pi*x) .* sin(6*pi*y);
  A = 4 + (x+2).^2 + y.^2;
  B = 1 + sin(x.*y);
  D = 2;
  Ax = 2*(x+2);
  Bx = y .* cos(x.*y);
  By = x .* cos(x.*y);
  f = -(Ax.*px + A.*pxx + Bx.*py + B.*pxy + By.*px + B.*pxy + D*pyy);
end

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

% Neumann boundary condition
% None 

% Known exact "pressure" solution 
function p = fxnp(pt)
  x = pt(:,1);  y = pt(:,2);
  p = sin(3*pi*x).^2 .* sin(3*pi*y).^2;
end

% Known gradient of the exact "pressure" solution 
function pg = fxnpg(pt) 
  x = pt(:,1);  y = pt(:,2);
  px = 3*pi * sin(6*pi*x) .* sin(3*pi*y).^2;
  py = 3*pi * sin(3*pi*x).^2 .* sin(6*pi*y);
  pg = [px,py];
end

% Known "velocity" function  
function u = fxnu(pt) 
  K = fxnK(pt);
  pg = fxnpg(pt);
  u = zeros(size(pt));
  u(:,1) = -( K(:,1,1).*pg(:,1) + K(:,1,2).*pg(:,2) );
  u(:,2) = -( K(:,2,1).*pg(:,1) + K(:,2,2).*pg(:,2) );  
end

end