Learn how easy it is to sync an existing GitHub or Google Code repo to a SourceForge project! See Demo

Close

Diff of /inst/example/NavierStokes/NS.m [000000] .. [d03d99] Maximize Restore

  Switch to side-by-side view

--- a
+++ b/inst/example/NavierStokes/NS.m
@@ -0,0 +1,102 @@
+pkg load fem-fenics msh
+import_ufl_Problem ("TentativeVelocity");
+import_ufl_Problem ("VelocityUpdate");
+import_ufl_Problem ("PressureUpdate");
+
+# We can either load the mesh from the file as in Dolfin but we can also use the msh pkg to generate the L-shape domain
+name = [tmpnam ".geo"];
+fid = fopen (name, "w");
+fputs (fid,"Point (1)  = {0, 0, 0, 0.1};\n");
+fputs (fid,"Point (2)  = {1, 0, 0, 0.1};\n");
+fputs (fid,"Point (3)  = {1, 0.5, 0, 0.1};\n");
+fputs (fid,"Point (4)  = {0.5, 0.5, 0, 0.1};\n");
+fputs (fid,"Point (5) = {0.5, 1, 0, 0.1};\n");
+fputs (fid,"Point (6) = {0, 1, 0,0.1};\n");
+
+fputs (fid,"Line (1)  = {5, 6};\n");
+fputs (fid,"Line (2) = {2, 3};\n");
+
+fputs (fid,"Line(3) = {6,1,2};\n");
+fputs (fid,"Line(4) = {5,4,3};\n");
+fputs (fid,"Line Loop(7) = {3,2,-4,1};\n");
+fputs (fid,"Plane Surface(8) = {7};\n");
+fclose (fid);
+msho = msh2m_gmsh (canonicalize_file_name (name)(1:end-4),"scale", 1,"clscale", .2);
+unlink (canonicalize_file_name (name));
+
+mesh = Mesh (msho);
+
+# Define function spaces (P2-P1). From ufl file
+V = FunctionSpace ('VelocityUpdate', mesh);
+Q = FunctionSpace ('PressureUpdate', mesh);
+
+# Set parameter values
+dt = 0.01;
+T = 3.;
+
+# Define boundary conditions
+noslip = DirichletBC (V, @(x,y) [0; 0], [3, 4]);
+outflow = DirichletBC (Q, @(x,y) 0, 2);
+
+# Create functions
+u0 = Expression ('u0', @(x,y) [0; 0]);
+
+# Define coefficients
+k = Constant ('k', dt);
+f = Constant ('f', [0; 0]);
+
+# Tentative velocity step.
+a1 = BilinearForm ('TentativeVelocity', V, k);
+
+# Pressure update.
+a2 = BilinearForm ('PressureUpdate', Q);
+
+# Velocity update
+a3 = BilinearForm ('VelocityUpdate', V);
+
+# Assemble matrices
+A1 = assemble (a1, noslip);
+A3 = assemble (a3, noslip);
+
+# Time-stepping
+t = dt; i = 0;
+while t < T
+
+  # Update pressure boundary condition
+  inflow = DirichletBC (Q, @(x,y) sin(3.0*t), 1);
+
+  # Compute tentative velocity step
+  L1 = LinearForm ('TentativeVelocity', V, k, u0, f);
+  b1 = assemble (L1, noslip);
+  utmp = A1 \ b1;
+  u1 = Function ('u1', V, utmp);
+
+  # Pressure correction
+  L2 = LinearForm ('PressureUpdate', Q, u1, k);
+  [A2, b2] = assemble_system (a2, L2, inflow, outflow);
+  ptmp = A2 \ b2;
+  p1 = Function ('p1', Q, ptmp);
+
+  # Velocity correction
+  L3 = LinearForm ('VelocityUpdate', V, k, u1, p1);
+  b3 = assemble (L3, noslip);
+  ut = A3 \ b3;
+  u1 = Function ('u0', V, ut);
+
+  # Plot solution
+  plot (p1);
+  plot (u1);
+
+  # Save to file
+  save (p1, sprintf ("p_%3.3d", ++i));
+  delete (sprintf ("p_%3.3d.pvd", i));
+  save (u1, sprintf ("u_%3.3d", i));
+  delete (sprintf ("u_%3.3d.pvd", i));
+
+  # Move to next time step
+  u0 = u1;
+  t += dt
+
+end
+
+