function EqnBC = EqnBC_Poisson_ArboWheeYotovSINUM1997_Ex2
%% EqnBC: Poisson: Adopted from ArboWheeYotov_SINUM_1997, Ex.2, p.848
% James Liu, ColoState; 2012/07--2018/12 

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

% 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) = 11;
  K(:,1,2) =  9;
  K(:,2,1) =  9;
  K(:,2,2) = 13;
end

% The right-hand sied function in the Poisson/Darcy equation
function f = fxnf(pt)
  x = pt(:,1);  y = pt(:,2);
  f11 = 11*(-2) * y.*(1-y);
  f12 = 9 * (1-2*x) .* (1-2*y);
  f21 = f12;
  f22 = 13 * x.*(1-x) * (-2);
  f = -(f11+f12+f21+f22);
end

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

% Neumann boundary condition
function uN = fxnuN(pt) 
  x = pt(:,1);  y = pt(:,2);
  ux = 11*(1-2*x).*y.*(1-y) + 9*x.*(1-x).*(1-2*y);
  ux = -ux;
  uy = 9*(1-2*x).*y.*(1-y) + 13*x.*(1-x).*(1-2*y);
  uy = -uy;
  uN = (x==0) .* (-ux)...
     + (x==1) .*   ux ...
     + (y==0) .* (-uy)...
     + (y==1) .*   uy;
end

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

% Known gradient of the exact "pressure" solution 
function pg = fxnpg(pt) 
  x = pt(:,1);  y = pt(:,2);
  px = (1-2*x) .* y.*(1-y);
  py = x.*(1-x) .* (1-2*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

% Known velocity divergence function: divu = f (source term)  
function divu = fxndivu(pt)
  x = pt(:,1);  y = pt(:,2);
  f11 = 11*(-2) * y.*(1-y);
  f12 = 9 * (1-2*x) .* (1-2*y);
  f21 = f12;
  f22 = 13 * x.*(1-x) * (-2);
  f = -(f11+f12+f21+f22);
  divu = f;
end

end