Mercurial > fem-fenics-eugenio
changeset 191:0c748179f6d4
Steady NS with no slip BC.
author | gedeone-octave <marcovass89@hotmail.it> |
---|---|
date | Sun, 10 Nov 2013 19:32:47 +0000 |
parents | 501e5eb8be33 |
children | 945c19831a16 |
files | devel/example/Ficticious_Domain/Steady_state/no_slip_bc/A.ufl devel/example/Ficticious_Domain/Steady_state/no_slip_bc/B.ufl devel/example/Ficticious_Domain/Steady_state/no_slip_bc/C.ufl devel/example/Ficticious_Domain/Steady_state/no_slip_bc/Err_p.ufl devel/example/Ficticious_Domain/Steady_state/no_slip_bc/Err_u.ufl devel/example/Ficticious_Domain/Steady_state/no_slip_bc/NS_with_no_slip_bc.m devel/example/Ficticious_Domain/Steady_state/no_slip_bc/generate_mesh.m |
diffstat | 7 files changed, 311 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/devel/example/Ficticious_Domain/Steady_state/no_slip_bc/A.ufl Sun Nov 10 19:32:47 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/no_slip_bc/B.ufl Sun Nov 10 19:32:47 2013 +0000 @@ -0,0 +1,8 @@ +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/no_slip_bc/C.ufl Sun Nov 10 19:32:47 2013 +0000 @@ -0,0 +1,9 @@ +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/no_slip_bc/Err_p.ufl Sun Nov 10 19:32:47 2013 +0000 @@ -0,0 +1,6 @@ +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/no_slip_bc/Err_u.ufl Sun Nov 10 19:32:47 2013 +0000 @@ -0,0 +1,6 @@ +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/no_slip_bc/NS_with_no_slip_bc.m Sun Nov 10 19:32:47 2013 +0000 @@ -0,0 +1,146 @@ +## 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 the flow around a square +# cylinder in a chennel. In this example, homogeneous no-slip BC are imposed +# on the square sides. + +# WARNING: read the warning inside generate_mesh.m + +pkg load fem-fenics msh + +generate_mesh; +mesh = Mesh (msh19); +plot (mesh); + +import_ufl_Problem ('A'); +import_ufl_BilinearForm ('B'); +import_ufl_BilinearForm ('C'); +import_ufl_Functional ('Err_u'); +import_ufl_Functional ('Err_p'); +import_ufl_FunctionSpace ('C'); + +TH1 = FunctionSpace ('A', mesh); +TH2 = FunctionSpace ('C', mesh); + +bc1 = DirichletBC (TH1, @(x, y, z, n) [0, 0], + [1, 5, 8, 19, 23, 26, 7, 16, 11, 21, 27, 30]); +bc2 = DirichletBC (TH1, @(x, y) [4*(1 - y)*(y) , 0], [4, 13, 20, 28]); +bc = {bc1, bc2}; + +u0 = Expression ('u0', @(x, y) [0; 0]); +f = Expression ('f', @(x, y) [0; 0]); + +nu = 1/40; +nu = Expression ('nu', @(x, y) nu); +sig = Expression ('sig', @(x, y) 0); + +a = BilinearForm ('A', TH1, TH1, nu, sig, u0); +L = LinearForm ('A', TH1, f); +[A, ff] = assemble_system (a, L, bc{:}); + +b = BilinearForm ('B', TH1, TH2); +B = assemble(b); + +m = BilinearForm ('C', 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, 200, 1e-6, 1000, 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)); + +#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', 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); +save (u, 'velocity'); +save (p, 'pressure'); + +dom = @(x, y) (x <= 1+r)*(x >= 1)*(y >= (0.5 - r/2))*(y <= (0.5 + r/2)); +norm_err = 0; +for i = 1:size(msh19.p, 2) + x = msh19.p (1, i); + y = msh19.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)
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/devel/example/Ficticious_Domain/Steady_state/no_slip_bc/generate_mesh.m Sun Nov 10 19:32:47 2013 +0000 @@ -0,0 +1,105 @@ +## 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/>. + + +# we generate a structured mesh of size [0, 4]x[0, 1] +# with a square hole at [1, 1.25]x[0.375, 0.625]. + +# WARNING: the following code is an attempt to create a structured mesh with +# an hole using the pkg msh. +# Unfortunately at the moment I've not yet been able to do it an a good way, +# and the following code will raise an error with the version 1.0.8 of msh. +# To avoid this error, you should comment the line number 132 in the file +# msh2m_join_structured_mesh.m. + +pkg load msh + +r = 0.25; +nel = 2; +n1 = nel*8 + 1; +n2 = nel*2 + 1; +n3 = nel*22 + 1; +n4 = nel*3 + 1; +n5 = nel*1 + 1; + +x = linspace (0, 1, n1); +y = linspace (0, 0.5 - r/2, n4); +msh1 = msh2m_structured_mesh (x, y, 1, 1:4); + +x = linspace (1, 1+r, n2); +msh2 = msh2m_structured_mesh (x, y, 1, 5:8); + +x = linspace (1+r, 4, n3); +msh3 = msh2m_structured_mesh (x, y, 1, 9:12); + +msh12 = msh2m_join_structured_mesh (msh1, msh2, 2, 8); +msh13 = msh2m_join_structured_mesh (msh12, msh3, 6, 12); + + +#### +x = linspace (0, 1, n1); +y = linspace (0.5 + r/2, 1, n4); +msh7 = msh2m_structured_mesh (x, y, 1, 25:28); + +x = linspace (1, 1+r, n2); +msh8 = msh2m_structured_mesh (x, y, 1, 29:32); +msh78 = msh2m_join_structured_mesh (msh7, msh8, 26, 32); + +x = linspace (1+r, 4, n3); +msh9 = msh2m_structured_mesh (x, y, 1, 33:36); +msh79 = msh2m_join_structured_mesh (msh78, msh9, 30, 36); + + +### +x = linspace (0, 1, n1); +y = linspace (0.5 - r/2, 0.5, n5); +msh4 = msh2m_structured_mesh (x, y, 1, 13:16); +msh14 = msh2m_join_structured_mesh (msh13, msh4, 3, 13); + +x = linspace (1, 1+r, n2); +msh5 = msh2m_structured_mesh (x, y, 1, 17:20); + +x = linspace (1+r, 4, n3); +msh6 = msh2m_structured_mesh (x, y, 1, 21:24); +msh16 = msh2m_join_structured_mesh (msh14, msh6, 10, 21); + +for i = 1:size(msh16.e (5, :),2) + if (msh16.e(5, i) == 15) + msh16.e(5, i) = 12; + endif +endfor + +#### +x = linspace (0, 1, n1); +y = linspace (0.5, 0.5 + r/2, n5); +msh4 = msh2m_structured_mesh (x, y, 1, 37:40); +msh74 = msh2m_join_structured_mesh (msh79, msh4, 25, 39); + +x = linspace (1, 1+r, n2); +msh5 = msh2m_structured_mesh (x, y, 1, 41:44); + +x = linspace (1+r, 4, n3); +msh6 = msh2m_structured_mesh (x, y, 1, 45:48); +msh76 = msh2m_join_structured_mesh (msh74, msh6, 32, 47); + +for i = 1:size(msh76.e (5, :),2) + if (msh76.e(5, i) == 38) + msh76.e(5, i) = 35; + endif +endfor + +########### +msh19 = msh2m_join_structured_mesh (msh16, msh76, 12, 35); +