function EqnBC = EqnBC_SteadyCD_Ex13Gauss
%% EqnBC: Steady-state convection-diffusion equation: Gaussian
% Domain [0,1]^2;  Concentration maximum 1;  Steep fronts in all directions
% James Liu, ColoState; 2012/07--2017/02 

D = 1;  % Diffusion coefficient
a = 1;  b = 1;  % Velocity vector v=[a,b]
xc = 0.25;  yc = 0.25;  sigma22 = 0.01;  % 2*sigma^2 

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

% Known velocity field 
function v = fxnv(pt) 
  x = pt(:,1);  % y = pt(:,2);
  v = [a*ones(size(x)), b*ones(size(x))];
end

% The right-hand function in the PDE 
function f = fxnf(pt)
  x = pt(:,1);  y = pt(:,2);
  tmp = ((x-xc).^2 +(y-yc).^2)/sigma22;
  tmp1 = (-a*(x-xc)-b*(y-yc)+2*D)*2/sigma22 - tmp*(4*D/sigma22);
  f = exp(-tmp) .* tmp1;  
end

% Dirichlet boundary condition for the concentration 
function cD = fxncD(pt) 
  % x = pt(:,1);  y = pt(:,2);
  % cD = zeros(size(x));
  cD = fxnc(pt);
end

% Neumann boundary condition for the concentration 
% --None 

% Known exact concentration 
function c = fxnc(pt)
  x = pt(:,1);  y = pt(:,2);
  tmp = (x-xc).^2 + (y-yc).^2;
  c = exp(-tmp/sigma22);
end

% Known exact concentration gradient 
function cg = fxncg(pt)
  x = pt(:,1);  y = pt(:,2);
  tmp = (x-xc).^2 + (y-yc).^2;
  cx = -(x-xc)*2/sigma22 .* exp(-tmp/sigma22);
  cy = -(y-yc)*2/sigma22 .* exp(-tmp/sigma22);
  cg = [cx,cy];
end

end