function EqnBC = EqnBC_Poisson_Ex31Ldom
%% EqnBC: Poisson: An L-shape problem, singularity at the corner 
% On L-shape domain [-1,1]^2\([0,1\times [-1,0]) 
% James Liu, ColoState; 2012/07--2017/02 

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 = zeros(size(x));
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);
  r = sqrt(x.^2+y.^2);
  theta = atan2(y,x);
  theta = (theta>=0).*theta + (theta<0).*(theta+2*pi);
  p = r.^(2/3).*sin(2*theta/3);
end

% Known exact "pressure" gradient 
function pg = fxnpg(pt)
  x = pt(:,1);  y = pt(:,2);
  r = sqrt(x.^2+y.^2);
  theta = atan2(y,x);
  theta = (theta>=0).*theta + (theta<0).*(theta+2*pi);
  pg(:,1) = 2/3*r.^(-1/3).*sin(2*theta/3).*x./r ...
          - 2/3*r.^(2/3).*cos(2*theta/3).*y./r.^2;
  pg(:,2) = 2/3*r.^(-1/3).*sin(2*theta/3).*y./r ...
          + 2/3*r.^(2/3).*cos(2*theta/3).*x./r.^2;
end

% Exact "velocity" UNKNOWN or just set to zero
function u = fxnu(pt) 
  u = zeros(size(pt));
end

end