function  EqnBIC = EqnBIC_TransCDR_CodinaCMAME1998_Ex
%% EqnBIC: TransCDR: 
% Example on the unit square ...
% Assuming a constant reaction coefficient
% Adopted from Codina_CMAME_1998 (p.205)
% James Liu, ColoState; 2012/07--2017/02

% These parameters can be adjusted 
T = 1;  % Time period [0,T] 
a = 1;  % a = 0.5;  % Magnitude of a constant velocity 
theta = 3*pi/8;  % Direction of the consatnt velocity
D = 1e-1;  % Diffusion coefficient
R = 1;  % R = 1e-4;  % Reaction coefficient 

EqnBIC = struct('T',T, ...
  'fxnv',@fxnv, 'D',D, 'R',R, 'fxnf',@fxnf, ...
  'fxncD',@fxncD, 'fxncN',@fxncN, 'fxnc0',@fxnc0, ...
  'fxnc',@fxnc);

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

% (1st order linear) reaction coefficient function
function Reac = fxnR(pt)
  % x = pt(:,1);  % y = pt(:,2);
  Reac = R;
end

% The right-hand side in the CDR equation
function f = fxnf(pt,t)
  x = pt(:,1);  % y = pt(:,2);
  f = ones(size(x));
end

% Dirichlet bndry. cond. for concentration 
function cD = fxncD(pt,t) 
  x = pt(:,1);  % y = pt(:,2);
  cD = zeros(size(x));
end

% Neumann boundary condition for the concentration 
% --None 

% Known initial condition for the concentration 
function c0 = fxnc0(pt)
  x = pt(:,1);  % y = pt(:,2);
  c0 = zeros(size(x));
end

% ??? No known exact concentration or just set to zero
function c = fxnc(pt,t)
  x = pt(:,1);  % y = pt(:,2);
  c = zeros(size(x));
end

end