function EqnBC = EqnBC_Poisson_Ex64Aniso
%% EqnBC: Poisson: Strongly discontinuous anisotropic diffusion, 4 pieces 
% Diagonally opposite source & sink: See EdwZheng_JCP_2008 Case4 
% James Liu, ColoState; 2012/07--2017/02 

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

% Permeability as full tensor (2x2 matrix) 
function K = fxnK(pt)
  x = pt(:,1);  y = pt(:,2);
  NumPts = size(pt,1);
  a = 0.246436002E+4;  b = 0.1148683643E+4;  d = 0.05366399794E+4;
  K1 = zeros(NumPts,2,2);  K2 = zeros(NumPts,2,2);  K = zeros(NumPts,2,2);
  K1(:,1,1) = a;  K1(:,1,2) =  b;  K1(:,2,1) =  b;  K1(:,2,2) = d;  
  K2(:,1,1) = a;  K2(:,1,2) = -b;  K2(:,2,1) = -b;  K2(:,2,2) = d;  
  Rgn1 = ((x-0.5).*(y-0.5)>0);  Rgn2 = ~Rgn1;
  K = repmat(Rgn1,2,2).*K1 + repmat(Rgn2,2,2).*K2;
end

% The right-hand side function in the Poisson/Darcy equation
function f = fxnf(pt)
  x = pt(:,1);  y = pt(:,2);
  xc1 = 0.25;  yc1 = 0.25;  sigma2 = 0.005;
  xc2 = 0.75;  yc2 = 0.75;  
  tmp1 = (x-xc1).^2 + (y-yc1).^2;
  tmp2 = (x-xc2).^2 + (y-yc2).^2;
  f1 =  exp(-0.5/sigma2*tmp1);  % source 
  f2 = -exp(-0.5/sigma2*tmp2);  % sink 
  f = 1e5 * (f1+f2);
end

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

% Neumann boundary condition 
% None

% No known exact "pressure" 
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