function EqnBC = EqnBC_Poisson_WihlerRivieJSC2011_Ex
%% EqnBC: Poisson: On [0,1]^2; Corner singularity at the origin 
% Cf. Wihler,Riviere, J.Sci.Comput., 2011, 
% James Liu, ColoState; 2012/07--2017/02 

% This parameter a (0<a<=1) can be changed 
a = 0.4;

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 = ones(size(x));
end

% Diffusion coeff. or permeability as 2x2 SPD matrix 
function K = fxnK(pt)
  NumPts = size(pt,1);
  K = zeros(NumPts,2,2);
  K(:,1,1) = 1;
  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 = -(x.^2+y.^2).^(-2+0.5*a).*(...
    2*x.^4-2*x.^3+(-a^2+6-2*a)*y.*x.^2 + (a^2-8+4*a)*y.^2.*x.^2 ...
    + (a^2-4)*y.*x + (-a^2+6-2*a)*y.^2.*x + 2*y.^4-2*y.^3 );
end

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

% Neumann boundary condition
% None 

% Known exact "pressure" 
function p = fxnp(pt)
  x = pt(:,1);  y = pt(:,2);
  r = sqrt(x.^2+y.^2);
  p = x.*(1-x).*y.*(1-y).*r.^(-2+a);
end

% Known exact "pressure" gradient
function pg = fxnpg(pt) 
  x = pt(:,1);  y = pt(:,2);
  r = sqrt(x.^2+y.^2);
  px = (1-2*x).*(r.^(-2+a)) + (x.^2).*(1-x).*(r.^(-4+a))*(-2+a);
  py = (1-2*y).*(r.^(-2+a)) + (y.^2).*(1-y).*(r.^(-4+a))*(-2+a);
  px = y.*(1-y).*px;  
  py = x.*(1-x).*py;
  pg = [px,py];
end

% Known exact "Darcy velocity" 
function u = fxnu(pt)
  NumPts = size(pt,1);
  u = zeros(NumPts,2);
  K = fxnK(pt);
  pg = fxnpg(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