Mercurial > fem-fenics-eugenio
view inst/example/NavierStokes/NS.m @ 164:7190852cf57f
Update the examples giving two input parameters to the BilinearForm.
author | gedeone-octave <marcovass89@hotmail.it> |
---|---|
date | Sat, 05 Oct 2013 10:40:06 +0100 |
parents | d03d998e64af |
children | e85390ed620d |
line wrap: on
line source
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, V, k); # Pressure update. a2 = BilinearForm ('PressureUpdate', Q, Q); # Velocity update a3 = BilinearForm ('VelocityUpdate', V, 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