function EqnBC = EqnBC_Poisson_Ex60Aniso
%% EqnBC: Poisson: Anisotropy in diffusion, Self-made 
% Purpose: Checking positivity and boundedness, i.e., in [0,1] 
% Known exact pressure:  p(x,y) = exp(-x/\epsilon) sin(\pi y) 
% James Liu, ColoState; 2012/07--2017/02 

% This parameter can be changed 
epsilon = 1e-6;

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

% Diffusion coeff. or permeability as order-2 diagoal matrix 
function K = fxnK(pt)
  % x = pt(:,1);  y = pt(:,2);
  NumPts = size(pt,1);
  K = zeros(NumPts,2,2);
  K(:,1,1) = epsilon;
  K(:,2,2) = 1;
end

% The right-hand side function in the Poisson/Darcy equation
function f = fxnf(pt)
  x = pt(:,1);  y = pt(:,2);
  f = (-1/epsilon + pi*pi) * exp(-x/epsilon) .* sin(pi*y);
end

% Dirichlet boundary condition
function pD = fxnpD(pt) 
  pD = fxnp(pt);
end

% Neumann boundary condition 
% None 

% Known exact "pressure" 
function p = fxnp(pt)
  x = pt(:,1);  y = pt(:,2);
  epsilon = 1e-6;  % It could be changed 
  p = exp(-x/epsilon) .* sin(pi*y);
end

% Known exact "pressure" gradient 
function pg = fxnpg(pt)
  x = pt(:,1);  y = pt(:,2);
  px = -1/epsilon * exp(-x/epsilon) .* sin(pi*y);
  py = pi * exp(-x/epsilon) .* cos(pi*y);
  pg = [px,py];
end

% Known exact "velocity" 
function u = fxnu(pt)
  x = pt(:,1);  y = pt(:,2);
  u(:,1) = exp(-x/epsilon) .* sin(pi*y);  % JL20151031: This is p(x,y)
  u(:,2) = -pi * exp(-x/epsilon) .* cos(pi*y);
end

end