### 55 lines (40 with data), 1.3 kB

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54``` ```# 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_BilinearForm (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 ```