--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/devel/example/Ficticious_Domain/Unsteady/L2_penalization/PressureUpdate.ufl	Sun Nov 10 19:33:44 2013 +0000
@@ -0,0 +1,30 @@
+## 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/>.
+
+
+V = VectorElement("CG", triangle, 2)
+Q = FiniteElement("CG", triangle, 1)
+
+# Define trial and test functions
+p = TrialFunction(Q)
+q = TestFunction(Q)
+
+# Define coefficients
+k  = Constant(triangle)
+u1 = Coefficient(V)
+
+# Define bilinear and linear forms
+a = inner(grad(p), grad(q))*dx
+L = -(1/k)*div(u1)*q*dx

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/devel/example/Ficticious_Domain/Unsteady/L2_penalization/TentativeVelocity.ufl	Sun Nov 10 19:33:44 2013 +0000
@@ -0,0 +1,37 @@
+## 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/>.
+
+
+V = VectorElement("CG", triangle, 2)
+Q = FiniteElement("CG", triangle, 1)
+P0 = FiniteElement("DG", triangle, 0)
+# Define trial and test functions
+u = TrialFunction(V)
+v = TestFunction(V)
+
+# Define coefficients
+k  = Constant(triangle)
+u0 = Coefficient(V)
+p0 = Coefficient(Q)
+f  = Coefficient(V)
+sig = Coefficient(P0)
+nu = Coefficient(P0)
+
+# Define bilinear and linear forms
+eq = (1/k)*inner(u - u0, v)*dx + inner(grad(u)*u0, v)*dx + \
+    nu*inner(grad(u), grad(v))*dx + sig * (inner (u, v))*dx \
+    - inner(f, v)*dx - div(v)*p0*dx
+a  = lhs(eq)
+L  = rhs(eq)

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/devel/example/Ficticious_Domain/Unsteady/L2_penalization/UNS_with_L2_penalization.m	Sun Nov 10 19:33:44 2013 +0000
@@ -0,0 +1,128 @@
+## 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 unsteady NS equation for a flow around a square 
+# cylinder in a channel. In this example, the no-slip BC are computed using
+# a L2 penalization technique.
+# We use the differential Chorin-Temam algorithm for each time step.
+# The solution shows the typical von-karman street pattern after a transition 
+# period.
+
+pkg load fem-fenics msh
+
+x = linspace (0, 10, 80*4+1);
+y = linspace (-2, 2, 80+1);
+msho = msh2m_structured_mesh (x, y, 1, 1:4);
+mesh = Mesh (msho);
+plot (mesh);
+import_ufl_Problem ("TentativeVelocity");
+import_ufl_Problem ("VelocityUpdate");
+import_ufl_Problem ("PressureUpdate");
+
+TH1 = FunctionSpace ('TentativeVelocity', mesh);
+TH2 = FunctionSpace ('PressureUpdate', mesh);
+
+bc1 = DirichletBC (TH1, @(x, y, z, n) [0, 0], [1, 3]);
+bc2 = DirichletBC (TH1, @(x, y) [1.5*4*(2 - y)*(2 + y)/16 , 0], 4);
+bc3 = DirichletBC (TH2, @(x, y) 0, 2);
+bcu = {bc1, bc2};
+bcp = bc3;
+## Set parameter values
+dt = 0.1;
+T = 50.;
+u0 = Expression ('u0', @(x, y) [0; 0]);
+p0 = Expression ('p0', @(x, y) 0);
+k = Constant ('k', dt);
+f = Constant ('f', [0; 0]);
+
+nu_1 = 5.e-3; 
+nu_0 = 5.e-3;
+r = 0.25;
+#ficticious domain
+
+## cylinder
+# dom = @(x,y) (((x - 3).^2 + (y).^2) < (r).^2);
+
+# square
+dom = @(x, y) (x <= 3+r)*(x >= 3-r)*(y >= -r)*(y <= r);
+nu = Expression ('nu', @(x, y) dom(x, y) * nu_0 + ...
+                            (1 - dom(x, y)) * nu_1);
+
+sig_1 = 0; 
+sig_0 = 1e8;
+sig = Expression ('sig',@(x, y) dom(x, y) * sig_0 + ... 
+                            (1 - dom(x, y)) * sig_1);
+
+
+a2 = BilinearForm ('PressureUpdate', TH2, TH2);
+A2 = assemble (a2, bcp);
+
+a3 = BilinearForm ('VelocityUpdate', TH1, TH1);
+A3 = assemble(a3);
+
+## Time-stepping
+t = dt; i = 0;
+pold = zeros (26001, 1);
+while t < T
+
+  # Compute tentative velocity step
+  a1 = BilinearForm ('TentativeVelocity', TH1, TH1, k, u0, sig, nu);
+  L1 = LinearForm ('TentativeVelocity', TH1, f, p0, k, u0);
+  [A1, b1] = assemble_system (a1, L1, bcu{:});
+  utmp = A1 \ b1;
+  u1 = Function ('u1', TH1, utmp);
+
+  # Pressure correction
+  L2 = LinearForm ('PressureUpdate', TH2, k, u1);
+  b2 = assemble (L2, bcp);
+  ptmp = A2 \ b2;
+  p1 = Function ('p1', TH2, ptmp);
+
+  # Velocity correction
+  L3 = LinearForm ('VelocityUpdate', TH1, k, u1, p1);
+  b3 = assemble (L3);
+  ut = A3 \ b3;
+  u1 = Function ('u0', TH1, ut);
+
+  #Update p
+  p = ptmp + pold;
+  p0 = Function ('p0', TH2, p);
+
+  # Save to file
+  save (p0, 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
+  pold = p;
+  u0 = u1;
+  t += dt;
+
+end
+
+plot (u0);
+
+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 (u0, [x, y]);
+    norm_err += err_L2(1).^2 + err_L2(2).^2;
+  endif
+endfor
+
+error_on_bc = sqrt (norm_err)

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/devel/example/Ficticious_Domain/Unsteady/L2_penalization/VelocityUpdate.ufl	Sun Nov 10 19:33:44 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/>.
+
+# Define function spaces (P2-P1)
+V = VectorElement("CG", triangle, 2)
+Q = FiniteElement("CG", triangle, 1)
+
+# Define trial and test functions
+u = TrialFunction(V)
+v = TestFunction(V)
+
+# Define coefficients
+k  = Constant(triangle)
+u1 = Coefficient(V)
+p1 = Coefficient(Q)
+
+# Define bilinear and linear forms
+a = inner(u, v)*dx
+L = inner(u1, v)*dx - k*inner(grad(p1), v)*dx