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