function EqnBC = EqnBC_Poisson_KelloggIntfc
%% EqnBC: Poisson: Kellogg interface problem 
% Example on the square [-1,1]^2 
% 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));
  idx = (x.*y >0);
  Ks(idx) = 161.4476387975881;
end

% Diffusion coeff. or permeability as 2x2 SPD matrix 
function K = fxnK(pt)
  Ks = fxnKs(pt);
  NumPts = size(pt,1);
  K = zeros(NumPts,2,2);
  K(:,1,1) = Ks(:);
  K(:,2,2) = Ks(:);
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);
  gamma = 0.1;
  sigma = -14.9225565104455152;
  rho = pi/4;
  r = sqrt(x.^2+y.^2);
  theta = atan2(y,x);
  theta = (theta>=0).*theta + (theta<0).*(theta+2*pi);
  mu = (theta>=0 & theta<pi/2).*cos((pi/2-sigma)*gamma).*cos((theta-pi/2+rho)*gamma)...
    + (theta>=pi/2 & theta<pi).*cos(rho*gamma).*cos((theta-pi+sigma)*gamma)...
    + (theta>=pi & theta<1.5*pi).*cos(sigma*gamma).*cos((theta-pi-rho)*gamma)...
    + (theta>=1.5*pi & theta<2*pi).*cos((pi/2-rho)*gamma).*cos((theta-1.5*pi-sigma)*gamma);
  p = -(r.^gamma).*mu;
end

% Known exact "pressure" gradient 
function pg = fxnpg(pt)
  x = pt(:,1);  y = pt(:,2);
  beta = 0.1;
  sigma = -14.9225565104455152;
  rho = pi/4;
  r = sqrt(x.^2+y.^2);
  theta = atan2(y,x);
  theta = (theta>=0).*theta + (theta<0).*(theta+2*pi);
  gx = (theta>=0 & theta<pi/2).*...
    (cos((pi/2-sigma)*beta).*(beta*x.*cos(beta*(theta-pi/2+rho))... 
      + beta*y.*sin(beta*(theta-pi/2+rho))))...
    + (theta>=pi/2 & theta<pi).*...
      (cos((rho)*beta).*(beta*x.*cos(beta*(theta-pi+sigma))...
         + beta*y.*sin(beta*(theta-pi+sigma))))...
    + (theta>=pi & theta<1.5*pi).*...
      (cos((sigma)*beta).*(beta*x.*cos(beta*(theta-pi-rho))...
        + beta*y.*sin(beta*(theta-pi-rho))))...
    + (theta>=1.5*pi & theta<2*pi).*...
      (cos((pi/2-rho).*beta)*(beta*x.*cos(beta*(theta-3*pi/2-sigma))... 
        + beta*y.*sin(beta*(theta-3*pi/2-sigma))));
  gx = r.^(beta-2).*gx;
  gy = (theta>=0 & theta<pi/2).*...
    (cos((pi/2-sigma)*beta).*(beta*y.*cos(beta*(theta-pi/2+rho))...
      - beta*x.*sin(beta*(theta-pi/2+rho))))...
    + (theta>=pi/2 & theta<pi).*...
      (cos((rho)*beta).*(beta*y.*cos(beta*(theta-pi+sigma))...
        - beta*x.*sin(beta*(theta-pi+sigma))))...
    + (theta>=pi & theta<1.5*pi).*...
      (cos((sigma)*beta).*(beta*y.*cos(beta*(theta-pi-rho))... 
        - beta*x.*sin(beta*(theta-pi-rho))))...
    + (theta>=1.5*pi & theta<2*pi).*...
      (cos((pi/2-rho).*beta)*(beta*y.*cos(beta*(theta-3*pi/2-sigma))...
        - beta*x.*sin(beta*(theta-3*pi/2-sigma))));
  gy= r.^(beta-2).*gy;
  pg = -[gx,gy];
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

end