Mercurial > fem-fenics-eugenio
changeset 102:d880c19f691b
Add comment to the script and pot the mesh.
* NS.m: added comment
* lshape-domain.m: plot the mesh
author | gedeone-octave <marcovass89@hotmail.it> |
---|---|
date | Wed, 14 Aug 2013 11:56:36 +0200 |
parents | 6d05e6c85724 |
children | 2d90592f73ab |
files | example/NavierStokes/NS.m example/NavierStokes/lshape-domain.m |
diffstat | 2 files changed, 10 insertions(+), 3 deletions(-) [+] |
line wrap: on
line diff
--- a/example/NavierStokes/NS.m Mon Aug 12 23:57:03 2013 +0200 +++ b/example/NavierStokes/NS.m Wed Aug 14 11:56:36 2013 +0200 @@ -7,13 +7,16 @@ fem_create_all ("VelocityUpdate"); fem_create_all ("PressureUpdate"); +# We use the fem pkg and the script above to generate the L-shape domain run lshape-domain.m; mshd = Mesh (msho); V = VelocityUpdate_FunctionSpace (mshd); Q = PressureUpdate_FunctionSpace (mshd); +# no slip BC inside bcnoslip = DirichletBC (V, @(x,y) [0;0], [3,4]); +# constant pressure at the outflow bcout = DirichletBC (Q, @(x,y) 0, 2); t = 0; @@ -31,19 +34,21 @@ A3 = assemble (a3, bcnoslip); # solve the problem for each time step -# We need to update only the lhs +# We use the Chorin-Temam algorithm i = 0; while t < T t += dt; + # time dependent pressure at the inflow bcin = DirichletBC (Q, @(x,y) sin(3.0*t), 1); - L1 = TentativeVelocity_LinearForm (V, k, u0, f); b1 = assemble (L1, bcnoslip); utmp = A1 \ b1; u1 = Function ('u1', V, utmp); + + # we need to assemble the system at each time-step because of time dependent BC a2 = PressureUpdate_BilinearForm (Q); L2 = PressureUpdate_LinearForm (Q,u1, k); [A2, b2] = assemble_system (a2, L2, bcin, bcout); @@ -57,6 +62,8 @@ ut = A3 \ b3; u0 = Function ('u0', V, ut); + # We use the fem_save () function to store the solution in vtu format name = sprintf ("u_%3.3d", ++i); fem_save (u0, name); + delete ([name ".pvd"]); end
--- a/example/NavierStokes/lshape-domain.m Mon Aug 12 23:57:03 2013 +0200 +++ b/example/NavierStokes/lshape-domain.m Wed Aug 14 11:56:36 2013 +0200 @@ -16,5 +16,5 @@ fputs (fid,"Plane Surface(8) = {7};\n"); fclose (fid); msho = msh2m_gmsh (canonicalize_file_name (name)(1:end-4),"scale", 1,"clscale", .2); -# trimesh (msho.t(1:3,:)', msho.p(1,:)', msho.p(2,:)'); + trimesh (msho.t(1:3,:)', msho.p(1,:)', msho.p(2,:)'); unlink (canonicalize_file_name (name));