```--- a
@@ -0,0 +1,146 @@
+## Copyright (C) 2013 Marco Vassallo <gedeone-octave@users.sourceforge.net>
+##
+## This program is free software; you can redistribute it and/or modify it under
+## Foundation; either version 3 of the License, or (at your option) any later
+## version.
+##
+## This program is distributed in the hope that it will be useful, but WITHOUT
+## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
+## details.
+##
+## You should have received a copy of the GNU General Public License along with
+## this program; if not, see <http://www.gnu.org/licenses/>.
+
+# In this example we solve the 2D NS equation for a the flow around a square
+# cylinder in a chennel. In this example, homogeneous no-slip BC are imposed
+# on the square sides.
+
+# WARNING: read the warning inside generate_mesh.m
+
+
+generate_mesh;
+mesh = Mesh (msh19);
+plot (mesh);
+
+import_ufl_Problem ('A');
+import_ufl_BilinearForm ('B');
+import_ufl_BilinearForm ('C');
+import_ufl_Functional ('Err_u');
+import_ufl_Functional ('Err_p');
+import_ufl_FunctionSpace ('C');
+
+TH1 = FunctionSpace ('A', mesh);
+TH2 = FunctionSpace ('C', mesh);
+
+bc1 = DirichletBC (TH1, @(x, y, z, n) [0, 0],
+                        [1, 5, 8, 19, 23, 26, 7, 16, 11, 21, 27, 30]);
+bc2 = DirichletBC (TH1, @(x, y) [4*(1 - y)*(y) , 0], [4, 13, 20, 28]);
+bc = {bc1, bc2};
+
+u0 = Expression ('u0', @(x, y) [0; 0]);
+f = Expression ('f', @(x, y) [0; 0]);
+
+nu = 1/40;
+nu = Expression ('nu', @(x, y) nu);
+sig = Expression ('sig', @(x, y) 0);
+
+a = BilinearForm ('A', TH1, TH1, nu, sig, u0);
+L = LinearForm ('A', TH1, f);
+[A, ff] = assemble_system (a, L, bc{:});
+
+b = BilinearForm ('B', TH1, TH2);
+B = assemble(b);
+
+m = BilinearForm ('C', TH2, TH2, nu);
+M = assemble(m);
+
+[x1, y1, v1] = find (A);
+[x2, y2, v2] = find (B');
+y2 += size (A, 1);
+[x3, y3, v3] = find (B);
+x3 += size (A, 1);
+[x4, y4, v4] = find (M);
+x4 += size (A, 1);
+y4 += size (A, 1);
+
+C = sparse ([x1; x2; x3],[y1; y2; y3],[v1; v2; v3],
+            (size (A,1) + size (B,1)), (size (A, 1) + size (B, 1)));
+
+P = sparse ([x1; x4],[y1; y4],[v1; v4],
+            (size (A,1) + size (B,1)), (size (A, 1) + size (B, 1)));
+
+F = [ff; (zeros (size (B, 1), 1))];
+
+[sol, flag, relres, iter, resvec] = gmres (C, F, 200, 1e-6, 1000, P);
+fprintf('Gmres converges in %d Iteration\n',iter (2));
+
+u = Function ('u', TH1, sol(1: (size(A,1))));
+p = Function ('p', TH2, sol((size(A,1))+1 : end));
+
+#Compute the initial norm
+p0 = Expression ('p0', @(x, y) 0);
+
+E1 = Functional ('Err_u', mesh, u, u0);
+normu0 = sqrt (assemble(E1));
+
+E2 = Functional ('Err_p', mesh, p, p0);
+normp0 = sqrt (assemble(E2));
+
+u0 = Function('u0',TH1, sol(1: (size(A,1))));
+err = 10;
+tol = 1e-4;
+maxit = 100;
+i = 1;
+
+while (err > tol && i < maxit)
+
+  a = BilinearForm ('A', TH1, TH1, nu, sig, u0);
+  [A, ff] = assemble_system (a, L, bc{:});
+
+  [x1, y1, v1] = find (A);
+  C = sparse ([x1; x2; x3],[y1; y2; y3],[v1; v2; v3],
+              (size (A,1) + size (B,1)), (size (A, 1) + size (B, 1)));
+  F = [ff; (zeros (size (B, 1), 1))];
+  [sol, flag, relres, iter, resvec] = gmres (C, F, 100, 1e-6, 2000, P);
+  fprintf('iteration %d: Gmres converges in %d Iteration\n',i, iter (2));
+
+  u = Function ('u', TH1, sol(1: (size(A,1))));
+  p = Function ('p', TH2, sol((size(A,1))+1 : end));
+
+  E1 = Functional ('Err_u', mesh, u, u0);
+  normu = sqrt (assemble(E1));
+
+  E2 = Functional ('Err_p', mesh, p, p0);
+  normp = sqrt (assemble(E2));
+
+  err = normu/normu0 + normp/normp0;
+  fprintf('the error is %f \n',err);
+  pause (0);
+  i++;
+  u0 = Function ('u0',TH1, sol(1: (size(A,1))));
+  p0 = Function ('p0', TH2, sol((size(A,1))+1 : end));
+  save (u, 'velocity');
+  save (p, 'pressure');
+
+end
+
+
+plot (u);
+save (u, 'velocity');
+save (p, 'pressure');
+
+dom = @(x, y) (x <= 1+r)*(x >= 1)*(y >= (0.5 - r/2))*(y <= (0.5 + r/2));
+norm_err = 0;
+for i = 1:size(msh19.p, 2)
+  x = msh19.p (1, i);
+  y = msh19.p (2, i);
+  if (dom (x, y) == true)
+    err_L2 = feval (u, [x, y]);
+    norm_err += err_L2(1).^2 + err_L2(2).^2;
+  endif
+endfor
+
+error_on_bc = sqrt (norm_err)
```