## 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)