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