Mercurial > fem-fenics-eugenio
view devel/example/Advection-Diffusion/AD_time.m @ 161:7e922376b3e8
Devel folder for examples which need to be improved
author | gedeone-octave <marcovass89@hotmail.it> |
---|---|
date | Thu, 12 Sep 2013 15:23:48 +0200 |
parents | example/Advection-Diffusion/AD_time.m@17358a4eb648 |
children |
line wrap: on
line source
# 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