function EqnBC = EqnBC_SteadyCD_BookOUP2005ElmanSilverWathen_Ex312
%% EqnBC: Steady conv.-diff.: Varying velocity, Chara.bndry.layers 
% Example on square [-1,1]^2 
% Adopted from ElmanSilverWathen_BookOUP_2005 (p.117): Ex.3.1.2 
% James Liu, ColoState; 2012/07--2017/02 

D = 2e-2;  % Diffusion coefficient 

EqnBC = struct('fxnv',@fxnv, 'D',D, 'fxnf',@fxnf, ...
  'fxncD',@fxncD, 'fxnuN',@fxnuN, 'fxnc',@fxnc);

% Known velocity field v=(0,1+((x+1)^2)/4) 
function v = fxnv(pt) 
  x = pt(:,1);  % y = pt(:,2);
  v = [zeros(size(x)), 1+((x+1).^2)/4];
end

% The right-hand function in steady convection-diffusion equation
function f = fxnf(pt)
  x = pt(:,1);  % y = pt(:,2);
  f = ones(size(x));
end

% Dirichlet boundary condition for the concentration 
function cD = fxncD(pt) 
  x = pt(:,1);  y = pt(:,2);
  cD = zeros(size(x));
  k = find(y==-1);  cD(k) = 1;  % Bottom: Inflow/entrance
  k = find(x==-1);  cD(k) = (1-(1+y(k))/2).^3;  % Left-wall: Decr. cubic. 
  k = find(x== 1);  cD(k) = (1-(1+y(k))/2).^2;  % Right wall: Decr. quad. 
end

% Neumann boundary condition for the concentration 
function uN = fxnuN(pt) 
  x = pt(:,1);  % y = pt(:,2);
  uN = zeros(size(x));
end

% ? No known exact concentration, so set to 0
function c = fxnc(pt)
  x = pt(:,1);  % y = pt(:,2);
  c = zeros(size(x));
  % For x: c = x - (exp(x/D)-1)/(exp(1/D)-1);
  % For y: ?
end

end