Mercurial > fem-fenics-eugenio

annotate inst/example/NavierStokes/NS.m @ 164:7190852cf57f

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
145
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
1 pkg load fem-fenics msh
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
2 import_ufl_Problem ("TentativeVelocity");
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
3 import_ufl_Problem ("VelocityUpdate");
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
4 import_ufl_Problem ("PressureUpdate");
76
6cb2a3e3be68 New examples which solve vector field problems
gedeone-octave <marcovass89@hotmail.it>
parents:
diff changeset
5
145
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
6 # 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
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
7 name = [tmpnam ".geo"];
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
8 fid = fopen (name, "w");
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
9 fputs (fid,"Point (1) = {0, 0, 0, 0.1};\n");
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
10 fputs (fid,"Point (2) = {1, 0, 0, 0.1};\n");
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
11 fputs (fid,"Point (3) = {1, 0.5, 0, 0.1};\n");
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
12 fputs (fid,"Point (4) = {0.5, 0.5, 0, 0.1};\n");
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
13 fputs (fid,"Point (5) = {0.5, 1, 0, 0.1};\n");
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
14 fputs (fid,"Point (6) = {0, 1, 0,0.1};\n");
77
ff95326e6f13 Make the example compliant with the fem-fenics pkg
gedeone-octave <marcovass89@hotmail.it>
parents: 76
diff changeset
15
145
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
16 fputs (fid,"Line (1) = {5, 6};\n");
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
17 fputs (fid,"Line (2) = {2, 3};\n");
76
6cb2a3e3be68 New examples which solve vector field problems
gedeone-octave <marcovass89@hotmail.it>
parents:
diff changeset
18
145
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
19 fputs (fid,"Line(3) = {6,1,2};\n");
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
20 fputs (fid,"Line(4) = {5,4,3};\n");
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
21 fputs (fid,"Line Loop(7) = {3,2,-4,1};\n");
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
22 fputs (fid,"Plane Surface(8) = {7};\n");
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
23 fclose (fid);
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
24 msho = msh2m_gmsh (canonicalize_file_name (name)(1:end-4),"scale", 1,"clscale", .2);
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
25 unlink (canonicalize_file_name (name));
76
6cb2a3e3be68 New examples which solve vector field problems
gedeone-octave <marcovass89@hotmail.it>
parents:
diff changeset
26
145
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
27 mesh = Mesh (msho);
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
28
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
29 # Define function spaces (P2-P1). From ufl file
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
30 V = FunctionSpace ('VelocityUpdate', mesh);
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
31 Q = FunctionSpace ('PressureUpdate', mesh);
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
32
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
33 # Set parameter values
76
6cb2a3e3be68 New examples which solve vector field problems
gedeone-octave <marcovass89@hotmail.it>
parents:
diff changeset
34 dt = 0.01;
6cb2a3e3be68 New examples which solve vector field problems
gedeone-octave <marcovass89@hotmail.it>
parents:
diff changeset
35 T = 3.;
6cb2a3e3be68 New examples which solve vector field problems
gedeone-octave <marcovass89@hotmail.it>
parents:
diff changeset
36
145
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
37 # Define boundary conditions
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
38 noslip = DirichletBC (V, @(x,y) [0; 0], [3, 4]);
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
39 outflow = DirichletBC (Q, @(x,y) 0, 2);
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
40
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
41 # Create functions
89
d06b423169fa The examples now follow the new naming convention.
gedeone-octave <marcovass89@hotmail.it>
parents: 77
diff changeset
42 u0 = Expression ('u0', @(x,y) [0; 0]);
95
77aabcc087dc Update the examples to the new naming convention.
gedeone-octave <marcovass89@hotmail.it>
parents: 89
diff changeset
43
145
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
44 # Define coefficients
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
45 k = Constant ('k', dt);
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
46 f = Constant ('f', [0; 0]);
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
47
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
48 # Tentative velocity step.
164
7190852cf57f Update the examples giving two input parameters to the BilinearForm.
gedeone-octave <marcovass89@hotmail.it>
parents: 159
diff changeset
49 a1 = BilinearForm ('TentativeVelocity', V, V, k);
145
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
50
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
51 # Pressure update.
164
7190852cf57f Update the examples giving two input parameters to the BilinearForm.
gedeone-octave <marcovass89@hotmail.it>
parents: 159
diff changeset
52 a2 = BilinearForm ('PressureUpdate', Q, Q);
95
77aabcc087dc Update the examples to the new naming convention.
gedeone-octave <marcovass89@hotmail.it>
parents: 89
diff changeset
53
145
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
54 # Velocity update
164
7190852cf57f Update the examples giving two input parameters to the BilinearForm.
gedeone-octave <marcovass89@hotmail.it>
parents: 159
diff changeset
55 a3 = BilinearForm ('VelocityUpdate', V, V);
145
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
56
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
57 # Assemble matrices
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
58 A1 = assemble (a1, noslip);
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
59 A3 = assemble (a3, noslip);
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
60
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
61 # Time-stepping
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
62 t = dt; i = 0;
76
6cb2a3e3be68 New examples which solve vector field problems
gedeone-octave <marcovass89@hotmail.it>
parents:
diff changeset
63 while t < T
145
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
64
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
65 # Update pressure boundary condition
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
66 inflow = DirichletBC (Q, @(x,y) sin(3.0*t), 1);
76
6cb2a3e3be68 New examples which solve vector field problems
gedeone-octave <marcovass89@hotmail.it>
parents:
diff changeset
67
145
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
68 # Compute tentative velocity step
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
69 L1 = LinearForm ('TentativeVelocity', V, k, u0, f);
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
70 b1 = assemble (L1, noslip);
76
6cb2a3e3be68 New examples which solve vector field problems
gedeone-octave <marcovass89@hotmail.it>
parents:
diff changeset
71 utmp = A1 \ b1;
89
d06b423169fa The examples now follow the new naming convention.
gedeone-octave <marcovass89@hotmail.it>
parents: 77
diff changeset
72 u1 = Function ('u1', V, utmp);
76
6cb2a3e3be68 New examples which solve vector field problems
gedeone-octave <marcovass89@hotmail.it>
parents:
diff changeset
73
145
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
74 # Pressure correction
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
75 L2 = LinearForm ('PressureUpdate', Q, u1, k);
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
76 [A2, b2] = assemble_system (a2, L2, inflow, outflow);
76
6cb2a3e3be68 New examples which solve vector field problems
gedeone-octave <marcovass89@hotmail.it>
parents:
diff changeset
77 ptmp = A2 \ b2;
89
d06b423169fa The examples now follow the new naming convention.
gedeone-octave <marcovass89@hotmail.it>
parents: 77
diff changeset
78 p1 = Function ('p1', Q, ptmp);
76
6cb2a3e3be68 New examples which solve vector field problems
gedeone-octave <marcovass89@hotmail.it>
parents:
diff changeset
79
145
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
80 # Velocity correction
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
81 L3 = LinearForm ('VelocityUpdate', V, k, u1, p1);
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
82 b3 = assemble (L3, noslip);
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
83 ut = A3 \ b3;
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
84 u1 = Function ('u0', V, ut);
76
6cb2a3e3be68 New examples which solve vector field problems
gedeone-octave <marcovass89@hotmail.it>
parents:
diff changeset
85
145
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
86 # Plot solution
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
87 plot (p1);
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
88 plot (u1);
76
6cb2a3e3be68 New examples which solve vector field problems
gedeone-octave <marcovass89@hotmail.it>
parents:
diff changeset
89
145
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
90 # Save to file
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
91 save (p1, sprintf ("p_%3.3d", ++i));
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
92 delete (sprintf ("p_%3.3d.pvd", i));
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
93 save (u1, sprintf ("u_%3.3d", i));
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
94 delete (sprintf ("u_%3.3d.pvd", i));
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
95
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
96 # Move to next time step
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
97 u0 = u1;
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
98 t += dt
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
99
100
5c5ef0f7498d save the solution instead of plotting it
gedeone-octave <marcovass89@hotmail.it>
parents: 95
diff changeset
100 end
145
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
101
93a4ee13aa75 Update the examples.
gedeone-octave <marcovass89@hotmail.it>
parents: 102
diff changeset
102