# the example is based on the one described here for the 
# bim pkg:
# http://wiki.octave.org/Bim_package#3D_Time_dependent_problem

pkg load msh fpl fem-fenics

problem = 'Advection_Diffusion';
import_ufl_Problem (problem);

problem = 'Reaction';
import_ufl_BilinearForm (problem);
import_ufl_FunctionSpace (problem);

x = linspace (0, 1, 40);
y = z = linspace (0, 1, 20);
msho = msh3m_structured_mesh (x, y, z, 1, 1:6);
mshd = Mesh (msho);

x = msho.p(1, :).';
y = msho.p(2, :).';
z = msho.p(3, :).';
 
x0 = .2; sx = .1;
y0 = .2; sy = .1;
z0 = .8; sz = .1;
 
u = exp (- ((x-x0)/(2*sx)) .^2 - ((y-y0)/(2*sy)) .^2 - ((z-z0)/(2*sz)) .^2);

V  = FunctionSpace ('Reaction', mshd);

f = Expression ('phi', @(x,y,z) x + y -z);
g = Expression ('mu', @(x,y,z) 0.01);

a = BilinearForm ('Advection_Diffusion', V, f, g);
m = BilinearForm ('Reaction', V);
A = assemble (a);
M = assemble (m);

# Lumped reaction matrix
ML = diag (sum (M, 1));

function du = f (u, t, A, M)
  du = - M \ (A * u);
endfunction 
 
time = linspace (0, 1, 30);
lsode_options ("integration method", "adams");
U = lsode (@(u, t) f(u, t, A, ML), u, time);
 
for ii = 1:1:numel (time)
  name = sprintf ("u_%3.3d", ii);
  delete ([name ".vtu"]);
  fpl_vtk_write_field (name, msho, {U(ii,:)', 'u'}, {}, 1);
endfor