function EqnBC = EqnBC_DarcyGen_GiraultVassiYotovNumerMath2014_Ex 
%% Darcy eqn. in a general form in domain (0,1)x(0,0.5)  
% Adopted from Girault,Vassilev,Yotov, Numer.Math., 2014 Ex.(p.145) 
% By Graham Harper, James Liu, ColoState; 2017/07--2020/12 

%% Parameters 
kappa = 1;  % For Darcy (corresponding to K=1 in the paper) 
mu = 0.1;  alpha = 1;
G = sqrt(mu*kappa)/alpha; 
xi = (1-G)/(2*(1+G));
chi = -(30*xi+17)/48; 
omega = 6;

EqnBC = struct('fxnK',@fxnK, 'fxns',@fxns, 'fxng',@fxng, ...
  'fxnpD',@fxnpD, 'fxnuN',@fxnuN, ...
  'fxnp',@fxnp, 'fxnu',@fxnu);

%% Darcy: Permeability as elementwise const 2x2 SPD matrix 
function K = fxnK(pt)
  x = pt(:,1);  % y = pt(:,2);
  K(:,1,1) = (kappa/mu) * ones(size(x));
  K(:,2,2) = (kappa/mu) * ones(size(x));
  K(:,1,2) = zeros(size(x));
  K(:,2,1) = zeros(size(x));
end

%% Darcy: Source term 
function s = fxns(pt)
  x = pt(:,1);  y = pt(:,2);
  s = (-omega^2)*sin(omega*x).*y + chi;
end

%% Darcy: An additional term, e.g., gravity 
function g = fxng(pt) 
  x = pt(:,1);  y = pt(:,2);
  u1 = omega * cos(omega*x) .* y;
  u2 = sin(omega*x) + chi*(y+0.5);
  g = (1 - 1/mu) * [u1,u2];
end

%% Darcy: Dirichlet boundary conditions
function pD = fxnpD(pt)
  x = pt(:,1);  y = pt(:,2);
  pD = -sin(omega*x).*y/kappa - chi/(2*kappa)*(y+0.5).^2;
end

%% Darcy: Neumann boundary condition
function uN = fxnuN(pt)
  x = pt(:,1);  y = pt(:,2);
  u = fxnu(pt);
  u1 = u(:,1);  u2 = u(:,2);
  uN = (-u1).*(x==0) + u1.*(x==1) + (-u2).*(y==0) + u2.*(y==0.5);
end

%% Darcy: Known exact solution for pressure 
function p = fxnp(pt)
  x = pt(:,1);  y = pt(:,2);  
  p = -sin(omega*x).*y/kappa - chi/(2*kappa)*(y+0.5).^2;
end

%% Darcy: Known exact solution for velocity 
function u = fxnu(pt)
  x = pt(:,1);  y = pt(:,2);
  u1 = omega * cos(omega*x) .* y;
  u2 = sin(omega*x) + chi*(y+0.5);
  u = [u1, u2];
end

end