--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/devel/example/Ficticious_Domain/Steady_state/L2_penalization/A_L2.ufl	Sun Nov 10 19:30:10 2013 +0000
@@ -0,0 +1,31 @@
+## Copyright (C) 2013 Marco Vassallo <gedeone-octave@users.sourceforge.net>
+##
+## This program is free software; you can redistribute it and/or modify it under
+## the terms of the GNU General Public License as published by the Free Software
+## 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/>.
+
+P2  = VectorElement("CG", triangle, 2)
+P1  = FiniteElement("CG", triangle, 1)
+P0  = FiniteElement("DG", triangle, 0)
+TH  = P2
+
+
+u = TrialFunction(TH)
+v = TestFunction(TH)
+
+nu  = Coefficient(P0)
+sig = Coefficient(P0)
+u0 = Coefficient(P2)
+f = Coefficient(P2)
+
+a = (nu * inner (grad (u), grad (v)) + sig * (inner (u, v)) + inner(grad(u)*u0, v)) * dx
+L = (inner(f, v))*dx

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/devel/example/Ficticious_Domain/Steady_state/L2_penalization/B_L2.ufl	Sun Nov 10 19:30:10 2013 +0000
@@ -0,0 +1,23 @@
+## Copyright (C) 2013 Marco Vassallo <gedeone-octave@users.sourceforge.net>
+##
+## This program is free software; you can redistribute it and/or modify it under
+## the terms of the GNU General Public License as published by the Free Software
+## 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/>.
+
+P2  = VectorElement("CG", triangle, 2)
+P1  = FiniteElement("CG", triangle, 1)
+P0  = FiniteElement("DG", triangle, 0)
+
+u = TrialFunction(P2)
+q = TestFunction(P1)
+
+a = (- div(u) * q)* dx

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/devel/example/Ficticious_Domain/Steady_state/L2_penalization/C_L2.ufl	Sun Nov 10 19:30:10 2013 +0000
@@ -0,0 +1,24 @@
+## Copyright (C) 2013 Marco Vassallo <gedeone-octave@users.sourceforge.net>
+##
+## This program is free software; you can redistribute it and/or modify it under
+## the terms of the GNU General Public License as published by the Free Software
+## 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/>.
+
+P1  = FiniteElement("CG", triangle, 1)
+P0  = FiniteElement("DG", triangle, 0)
+
+p = TrialFunction(P1)
+q = TestFunction(P1)
+
+nu  = Coefficient(P0)
+
+a = ((1/nu)*p*q) * dx

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/devel/example/Ficticious_Domain/Steady_state/L2_penalization/Err_p.ufl	Sun Nov 10 19:30:10 2013 +0000
@@ -0,0 +1,21 @@
+## Copyright (C) 2013 Marco Vassallo <gedeone-octave@users.sourceforge.net>
+##
+## This program is free software; you can redistribute it and/or modify it under
+## the terms of the GNU General Public License as published by the Free Software
+## 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/>.
+
+P0  = FiniteElement("DG", triangle, 0)
+
+p = Coefficient(P0)
+p0 = Coefficient(P0)
+
+M = ((p-p0)*(p-p0))*dx

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/devel/example/Ficticious_Domain/Steady_state/L2_penalization/Err_u.ufl	Sun Nov 10 19:30:10 2013 +0000
@@ -0,0 +1,21 @@
+## Copyright (C) 2013 Marco Vassallo <gedeone-octave@users.sourceforge.net>
+##
+## This program is free software; you can redistribute it and/or modify it under
+## the terms of the GNU General Public License as published by the Free Software
+## 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/>.
+
+P2 = VectorElement("CG", triangle, 2)
+
+u = Coefficient(P2)
+u0 = Coefficient(P2)
+
+M = (inner(u-u0,u-u0) + inner(grad(u -u0), grad(u-u0)))*dx

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/devel/example/Ficticious_Domain/Steady_state/L2_penalization/NS_with_L2_penalization.m	Sun Nov 10 19:30:10 2013 +0000
@@ -0,0 +1,148 @@
+## Copyright (C) 2013 Marco Vassallo <gedeone-octave@users.sourceforge.net>
+##
+## This program is free software; you can redistribute it and/or modify it under
+## the terms of the GNU General Public License as published by the Free Software
+## 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 flow around a square 
+# cylinder in a channel. In this example, the no-slip BC are applied using
+# a L2 penalization technique.
+# We use the preconditioned gmres algorithm.
+pkg load fem-fenics msh
+
+x = linspace (0, 4, 33);
+y = linspace (0, 1, 9);
+msho = msh2m_structured_mesh (x, y, 1, 1:4);
+mesh = Mesh (msho);
+
+import_ufl_Problem ('A_L2');
+import_ufl_BilinearForm ('B_L2');
+import_ufl_BilinearForm ('C_L2');
+import_ufl_Functional ('Err_u');
+import_ufl_Functional ('Err_p');
+import_ufl_FunctionSpace ('C_L2');
+TH1 = FunctionSpace ('A_L2', mesh);
+TH2 = FunctionSpace ('C_L2', mesh);
+
+bc1 = DirichletBC (TH1, @(x, y, z, n) [0, 0], [1, 3]);
+bc2 = DirichletBC (TH1, @(x, y) [4*(1 - y)*(y) , 0], 4);
+
+bc = {bc1, bc2};
+
+u0 = Expression ('u0', @(x, y) [0; 0]);
+f = Expression ('f', @(x, y) [0; 0]);
+
+nu_1 = 1/40; 
+nu_0 = 1/40;
+r = 0.25;
+#ficticious domain
+dom = @(x, y) (x <= 1+r)*(x >= 1)*(y >= (0.5 - r/2))*(y <= (0.5 + r/2));
+nu = Expression ('nu', @(x, y) dom(x, y) * nu_0 + ... 
+                            (1 - dom(x, y)) * nu_1);
+
+sig_1 = 0; 
+sig_0 = 1e4;
+sig = Expression ('sig',@(x, y) dom(x, y) * sig_0 + ... 
+                            (1 - dom(x, y)) * sig_1);
+
+a = BilinearForm ('A_L2', TH1, TH1, nu, sig, u0);
+L = LinearForm ('A_L2', TH1, f);
+[A, ff] = assemble_system (a, L, bc{:});
+
+b = BilinearForm ('B_L2', TH1, TH2);
+B = assemble(b);
+
+m = BilinearForm ('C_L2', 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, [], 1e-6, 100, 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));
+save (u, 'velocity');
+save (p, 'pressure');
+
+#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_L2', 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);
+
+norm_err = 0;
+for i = 1:size(msho.p, 2)
+  x = msho.p (1, i);
+  y = msho.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)