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));