# HG changeset patch # User gedeone-octave # Date 1384111842 0 # Node ID 501e5eb8be335f939ab2d5d11ee47f6e611a198b # Parent 98f598376b3a45459b6acf0191b993e3907fd0bd Steady NS with H1 penalization. diff -r 98f598376b3a -r 501e5eb8be33 devel/example/Ficticious_Domain/Steady_state/H1_penalization/A_H1.ufl --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/devel/example/Ficticious_Domain/Steady_state/H1_penalization/A_H1.ufl Sun Nov 10 19:30:42 2013 +0000 @@ -0,0 +1,31 @@ +## Copyright (C) 2013 Marco Vassallo +## +## 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 . + +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 * (1 + sig) * inner (grad (u), grad (v)) + sig * (inner (u, v)) + inner(grad(u)*u0, v)) * dx +L = (inner(f, v))*dx diff -r 98f598376b3a -r 501e5eb8be33 devel/example/Ficticious_Domain/Steady_state/H1_penalization/B_H1.ufl --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/devel/example/Ficticious_Domain/Steady_state/H1_penalization/B_H1.ufl Sun Nov 10 19:30:42 2013 +0000 @@ -0,0 +1,23 @@ +## Copyright (C) 2013 Marco Vassallo +## +## 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 . + +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 diff -r 98f598376b3a -r 501e5eb8be33 devel/example/Ficticious_Domain/Steady_state/H1_penalization/C_H1.ufl --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/devel/example/Ficticious_Domain/Steady_state/H1_penalization/C_H1.ufl Sun Nov 10 19:30:42 2013 +0000 @@ -0,0 +1,24 @@ +## Copyright (C) 2013 Marco Vassallo +## +## 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 . + +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 diff -r 98f598376b3a -r 501e5eb8be33 devel/example/Ficticious_Domain/Steady_state/H1_penalization/Err_p.ufl --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/devel/example/Ficticious_Domain/Steady_state/H1_penalization/Err_p.ufl Sun Nov 10 19:30:42 2013 +0000 @@ -0,0 +1,21 @@ +## Copyright (C) 2013 Marco Vassallo +## +## 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 . + +P0 = FiniteElement("DG", triangle, 0) + +p = Coefficient(P0) +p0 = Coefficient(P0) + +M = ((p-p0)*(p-p0))*dx diff -r 98f598376b3a -r 501e5eb8be33 devel/example/Ficticious_Domain/Steady_state/H1_penalization/Err_u.ufl --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/devel/example/Ficticious_Domain/Steady_state/H1_penalization/Err_u.ufl Sun Nov 10 19:30:42 2013 +0000 @@ -0,0 +1,22 @@ +## Copyright (C) 2013 Marco Vassallo +## +## 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 . + + +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 diff -r 98f598376b3a -r 501e5eb8be33 devel/example/Ficticious_Domain/Steady_state/H1_penalization/NS_with_H1_penalization.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/devel/example/Ficticious_Domain/Steady_state/H1_penalization/NS_with_H1_penalization.m Sun Nov 10 19:30:42 2013 +0000 @@ -0,0 +1,151 @@ +## Copyright (C) 2013 Marco Vassallo +## +## 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 . + +# In this example we solve the steady 2D NS equation for a flow around a square +# cylinder in a channel. In this example, the no-slip BC are applied using +# a H1 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); +plot (mesh); + +import_ufl_Problem ('A_H1'); +import_ufl_BilinearForm ('B_H1'); +import_ufl_BilinearForm ('C_H1'); +import_ufl_Functional ('Err_u'); +import_ufl_Functional ('Err_p'); +import_ufl_FunctionSpace ('C_H1'); +TH1 = FunctionSpace ('A_H1', mesh); +TH2 = FunctionSpace ('C_H1', 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_H1', TH1, TH1, nu, sig, u0); +L = LinearForm ('A_H1', TH1, f); +[A, ff] = assemble_system (a, L, bc{:}); + +b = BilinearForm ('B_H1', TH1, TH2); +B = assemble(b); + +m = BilinearForm ('C_H1', 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)))); + +#Iteration +err = 10; +tol = 1e-4; +maxit = 100; +i = 1; + +while (err > tol && i < maxit) + + a = BilinearForm ('A_H1', 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)