Mercurial > fem-fenics-eugenio
changeset 76:6cb2a3e3be68
New examples which solve vector field problems
* NS: solve the incompressible NS equations using the chorin-temam algorithm
* MixedPoisson: solve the Poisson equation using amixed formulation using DG and
Brezzi-Douglas-marini finite elements
author | gedeone-octave <marcovass89@hotmail.it> |
---|---|
date | Fri, 02 Aug 2013 23:01:39 +0200 |
parents | d6df5cc8ef53 |
children | ff95326e6f13 |
files | example/MixedPoisson.m example/MixedPoisson.ufl example/NavierStokes/NS.m example/NavierStokes/PressureUpdate.ufl example/NavierStokes/TentativeVelocity.ufl example/NavierStokes/VelocityUpdate.ufl example/NavierStokes/lshape-domain.m |
diffstat | 7 files changed, 234 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/example/MixedPoisson.m Fri Aug 02 23:01:39 2013 +0200 @@ -0,0 +1,25 @@ +pkg load msh + +fem_init_env (); + +x = y = linspace (0, 1, 32); +msho = msh2m_structured_mesh (x, y, 1, 1:4); +mshd = fem_init_mesh (msho); + +V = fem_fs_MixedPoisson (mshd); + +bc1 = fem_bc (V, @(x,y) [0; -sin(5.0*x); 0], 1); +bc2 = fem_bc (V, @(x,y) [0; sin(5.0*x); 0], 3); + + +f = fem_coeff ('f', @(x,y) 10*exp(-((x - 0.5)^2 + (y - 0.5)^2) / 0.02)); + +A = fem_rhs_MixedPoisson (V, bc1, bc2); +b = fem_lhs_MixedPoisson (V, f, bc1, bc2); + +u = A \ b; + + +func = fem_func ('u', V, u); + +fem_plot (func);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/example/MixedPoisson.ufl Fri Aug 02 23:01:39 2013 +0200 @@ -0,0 +1,12 @@ +BDM = FiniteElement("BDM", triangle, 1) +DG = FiniteElement("DG", triangle, 0) +W = BDM * DG + +(sigma, u) = TrialFunctions(W) +(tau, v) = TestFunctions(W) + +CG = FiniteElement("CG", triangle, 1) +f = Coefficient(CG) + +a = (dot(sigma, tau) + div(tau)*u + div(sigma)*v)*dx +L = - f*v*dx
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/example/NavierStokes/NS.m Fri Aug 02 23:01:39 2013 +0200 @@ -0,0 +1,49 @@ +pkg load msh + +fem_init_env (); + +run lshape-domain.m; +mshd = fem_init_mesh (msho); + +V = fem_fs_VelocityUpdate (mshd); +Q = fem_fs_PressureUpdate (mshd); + +bcnoslip = fem_bc (V, @(x,y) [0;0], [3,4]); +bcout = fem_bc (Q, @(x,y) 0, 2); + +t = 0; +dt = 0.01; +T = 3.; + +k = fem_coeff ('k', @(x,y) dt); +f = fem_coeff ('f', @(x,y) [0; 0]); +u0 = fem_coeff ('u0', @(x,y) [0; 0]); +p0 = fem_coeff ('p1', @(x,y) 0); +A1 = fem_rhs_TentativeVelocity (V, k, bcnoslip); +A3 = fem_rhs_VelocityUpdate (V, bcnoslip); +# solve the problem for each time step +# We need to update only the lhs +while t < T + t += dt; + bcin = fem_bc (Q, @(x,y) sin(3.0*t), 1); + + b1 = fem_lhs_TentativeVelocity (V, k, u0, f, bcnoslip); + + utmp = A1 \ b1; + + u1 = fem_func ('u1', V, utmp); + + A2 = fem_rhs_PressureUpdate (Q, bcin, bcout); + b2 = fem_lhs_PressureUpdate (Q, bcin, bcout, u1, k); + + ptmp = A2 \ b2; + p1 = fem_func ('p1', Q, ptmp); + + b3 = fem_lhs_VelocityUpdate (V, k, u1, p1, bcnoslip); + + ut = A3 \ b3; + u0 = fem_func ('u0', V, ut); + + #press Q to make the plot continue + fem_plot (u0); +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/example/NavierStokes/PressureUpdate.ufl Fri Aug 02 23:01:39 2013 +0200 @@ -0,0 +1,41 @@ +# Copyright (C) 2010 Anders Logg +# +# This file is part of DOLFIN. +# +# DOLFIN is free software: you can redistribute it and/or modify +# it under the terms of the GNU Lesser General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# DOLFIN 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 Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public License +# along with DOLFIN. If not, see <http://www.gnu.org/licenses/>. +# +# First added: 2010-08-30 +# Last changed: 2010-09-01 +# +# The bilinear form a(u, v) and linear form L(v) for the pressure +# update step in Chorin's method for the incompressible Navier-Stokes +# equations. +# +# Compile this form with FFC: ffc -l dolfin PressureUpdate.ufl. + +# Define function spaces (P2-P1) +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/example/NavierStokes/TentativeVelocity.ufl Fri Aug 02 23:01:39 2013 +0200 @@ -0,0 +1,45 @@ +# Copyright (C) 2010 Anders Logg +# +# This file is part of DOLFIN. +# +# DOLFIN is free software: you can redistribute it and/or modify +# it under the terms of the GNU Lesser General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# DOLFIN 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 Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public License +# along with DOLFIN. If not, see <http://www.gnu.org/licenses/>. +# +# First added: 2010-08-30 +# Last changed: 2011-06-30 +# +# The bilinear form a(u, v) and linear form L(v) for the tentative +# velocity step in Chorin's method for the incompressible +# Navier-Stokes equations. +# +# Compile this form with FFC: ffc -l dolfin TentativeVelocity.ufl. + +# 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) +u0 = Coefficient(V) +f = Coefficient(V) +nu = 0.01 + +# Define bilinear and linear forms +eq = (1/k)*inner(u - u0, v)*dx + inner(grad(u0)*u0, v)*dx + \ + nu*inner(grad(u), grad(v))*dx - inner(f, v)*dx +a = lhs(eq) +L = rhs(eq)
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/example/NavierStokes/VelocityUpdate.ufl Fri Aug 02 23:01:39 2013 +0200 @@ -0,0 +1,42 @@ +# Copyright (C) 2010 Anders Logg +# +# This file is part of DOLFIN. +# +# DOLFIN is free software: you can redistribute it and/or modify +# it under the terms of the GNU Lesser General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# DOLFIN 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 Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public License +# along with DOLFIN. If not, see <http://www.gnu.org/licenses/>. +# +# First added: 2010-08-30 +# Last changed: 2010-09-01 +# +# The bilinear form a(u, v) and linear form L(v) for the velocity +# update step in Chorin's method for the incompressible Navier-Stokes +# equations. +# +# Compile this form with FFC: ffc -l dolfin TentativeVelocity.ufl. + +# 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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/example/NavierStokes/lshape-domain.m Fri Aug 02 23:01:39 2013 +0200 @@ -0,0 +1,20 @@ +name = [tmpnam ".geo"]; +fid = fopen (name, "w"); +fputs (fid,"Point (1) = {0, 0, 0, 0.1};\n"); +fputs (fid,"Point (2) = {1, 0, 0, 0.1};\n"); +fputs (fid,"Point (3) = {1, 0.5, 0, 0.1};\n"); +fputs (fid,"Point (4) = {0.5, 0.5, 0, 0.1};\n"); +fputs (fid,"Point (5) = {0.5, 1, 0, 0.1};\n"); +fputs (fid,"Point (6) = {0, 1, 0,0.1};\n"); + +fputs (fid,"Line (1) = {5, 6};\n"); +fputs (fid,"Line (2) = {2, 3};\n"); + +fputs (fid,"Line(3) = {6,1,2};\n"); +fputs (fid,"Line(4) = {5,4,3};\n"); +fputs (fid,"Line Loop(7) = {3,2,-4,1};\n"); +fputs (fid,"Plane Surface(8) = {7};\n"); + fclose (fid); + msho = msh2m_gmsh (canonicalize_file_name (name)(1:end-4),"scale", 1,"clscale", .2); +# trimesh (msho.t(1:3,:)', msho.p(1,:)', msho.p(2,:)'); + unlink (canonicalize_file_name (name));