changeset 63:3231f888f25d

New examples for the fem-fenics pkg * Laplace: solve the laplace problem described on dolfin website * Biharmonic: solve the biharmonic equation described on dolfin * Heat: solve the heat evolutionary equation
author gedeone-octave <marco.vassallo@outlook.com>
date Fri, 26 Jul 2013 11:29:29 +0200
parents a00f95c0b048
children 55ffedcc59d4
files example/Biharmonic.ufl example/Heat.ufl example/Laplace.ufl example/fem_func_space.m example/test_biharmonic.m example/test_heat.m example/test_laplace.m obsolete/fem_func_space.m
diffstat 8 files changed, 290 insertions(+), 75 deletions(-) [+]
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/example/Biharmonic.ufl	Fri Jul 26 11:29:29 2013 +0200
@@ -0,0 +1,24 @@
+# Elements
+element = FiniteElement("Lagrange", triangle, 2)
+
+# Trial and test functions
+u = TrialFunction(element)
+v = TestFunction(element)
+f = Coefficient(element)
+
+# Normal component, mesh size and right-hand side
+n  = element.cell().n
+h = 2.0*triangle.circumradius
+h_avg = (h('+') + h('-'))/2
+
+# Parameters
+alpha = Constant(triangle)
+
+# Bilinear form
+a = inner(div(grad(u)), div(grad(v)))*dx \
+  - inner(avg(div(grad(u))), jump(grad(v), n))*dS \
+  - inner(jump(grad(u), n), avg(div(grad(v))))*dS \
+  + alpha('+')/h_avg*inner(jump(grad(u), n), jump(grad(v),n))*dS
+
+# Linear form
+L = f*v*dx
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/example/Heat.ufl	Fri Jul 26 11:29:29 2013 +0200
@@ -0,0 +1,13 @@
+element = FiniteElement("Lagrange", triangle, 1)
+
+u = TrialFunction(element)
+v = TestFunction(element)
+
+u0 = Coefficient(element)
+f = Coefficient(element)
+g = Coefficient(element)
+k = Constant(triangle)
+
+eq = (1/k)*(u-u0)*v*dx + inner(grad(u), grad(v))*dx -f*v*dx - g*v*ds
+a = rhs (eq)
+L = lhs (eq)
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/example/Laplace.ufl	Fri Jul 26 11:29:29 2013 +0200
@@ -0,0 +1,11 @@
+element = FiniteElement("Lagrange", triangle, 1)
+
+
+u = TrialFunction(element)
+v = TestFunction(element)
+
+f = Coefficient(element)
+g = Coefficient(element)
+
+a = inner(grad(u), grad(v))*dx
+L = f*v*dx + g*v*ds
--- a/example/fem_func_space.m	Fri Jul 26 10:26:16 2013 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,75 +0,0 @@
-## 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 2 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/>.
-
-
-## -*- texinfo -*-
-## @deftypefn {Function File} {} = fem_ffc (myproblem.ufl)
-## This function take as input a .ufl file with a Variationa Problem
-## @example
-## 
-## @end example
-## @seealso{} 
-## @end deftypefn
-
-function res = fem_func_space (mesh, name)
-
-  if nargin != 2
-    error ("fem_ffc: wrong number of input parameters.");
-  elseif ! isstr (name)
-    error ("fem_ffc: first argument is not a valid string");
-  endif
-
-  # fem_fs.cc is the file where I should include name.h
-  # I use name.h at line 35, where I call
-  # "new FunctionSpace (mshd);" 
-  # because FunctionSpace is defined inside of name.h
-  # (correct call would be "name::FunctionSpace (mshd)" but 
-  # I'm using the namespace "name")
-  command = "cp ../src/fem_fs.cc fem_fs.cc";
-  [output, text] = system (command);
-  if output != 0
-    display (text);
-    error ("Compilation failed");
-  endif
-
-  # add line to fem_fs.cc
-  filename = "fem_fs.cc";
-  fid = fopen (filename, "r+");
-  fseek (fid, 671, SEEK_SET ());
-  fprintf (fid, "%c%s.h%c", 34, name, 34);
-  fseek (fid, 721, SEEK_SET ());
-  fprintf (fid, "%s%c", name, 59);
-  fclose (fid);
-
-  #now compile the program
-  command = "mkoctfile -c fem_fs.cc $(LDFLAGS) -o fem_fs.o -I../src/";
-  [output, text] = system (command);
-  if output != 0
-    display (text);
-    error ("Compilation failed");
-  endif
-
-  command = "mkoctfile -s ../src/mesh.o ../src/fem_init_env.o fem_fs.o -o fem_fs.oct $(LDFLAGS) -ldolfin";
-  [output, text] = system (command);
-  if output != 0
-    display (text);
-    error ("Compilation failed");
-  endif
-
-  #run the program and get the result
-  res = fem_fs (mesh);
-
-endfunction
-
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/example/test_biharmonic.m	Fri Jul 26 11:29:29 2013 +0200
@@ -0,0 +1,49 @@
+# We follow the dolfin example for the Biharmonic problem
+# See (http://fenicsproject.org/documentation/dolfin/1.2.0/cpp/demo/pde/biharmonic/cpp/documentation.html#index-0)
+
+# initialize the environment and create function related to the problem
+# defined inside Biharmonic.ufl
+pkg load msh
+fem_init_env ();
+problem = 'Biharmonic';
+fem_create_fs (problem);
+fem_create_rhs (problem);
+fem_create_lhs (problem);
+
+# create a unit square mesh using msh: labels for the boundary sides are 1,2,3,4
+x = y = linspace (0, 1, 32);
+msho = msh2m_structured_mesh (x, y, 1, 1:4);
+mshd = fem_init_mesh (msho);
+
+V  = fem_fs_Biharmonic (mshd);
+
+# fem_bc takes as input the functionspace V, a function handler f,
+# and the sides where we want to apply the condition
+# The value on each point of the boundary is computed using
+# the eval method available inside expression.h
+f = @(x,y) 0;
+bc = fem_bc (V, f, 1:4);
+
+# fem_coeff takes as input a string and a function handler
+# and is used below to set the value of the coefficient of the rhs
+
+ff = @(x,y) 4.0*pi^4.*sin(pi.*x).*sin(pi.*y);
+f = fem_coeff ('f', ff);
+
+gg = @(x,y) 8;
+alpha = fem_coeff ('alpha', gg);
+
+# fem_rhs and fem_lhs takes as input the functionspace V, and the
+# boundarycondition bc and solve the Poisson problem with
+# the velues specified inside f and alpha;
+A = fem_rhs_Biharmonic (V, alpha, bc);
+
+b = fem_lhs_Biharmonic (V, f, bc);
+
+u = A \ b;
+
+func = fem_func (V, u)
+
+fem_plot (func);
+
+fem_save (func, problem);
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/example/test_heat.m	Fri Jul 26 11:29:29 2013 +0200
@@ -0,0 +1,60 @@
+# We follow the evolutionary heat equation specified in the file Heat.ufl
+
+pkg load msh
+fem_init_env ();
+problem = 'Heat';
+fem_create_fs (problem);
+fem_create_rhs (problem);
+fem_create_lhs (problem);
+
+# create a unit square mesh using msh: labels for the boundary sides are 1,2,3,4
+msho = msh2m_structured_mesh (0:0.05:1, 0:0.05:1, 1, 1:4);
+mshd = fem_init_mesh (msho);
+
+V  = fem_fs_Heat (mshd);
+
+# fem_bc takes as input the functionspace V, a function handler f,
+# and the sides where we want to apply the condition
+# The value on each point of the boundary is computed using
+# the eval method available inside expression.h
+# if a side is not specified, Neumann conditions are applied
+# with g specified below
+f = @(x,y) 0;
+bc = fem_bc (V, f, 1:4);
+
+# fem_coeff takes as input a string and a function handler
+# and is used below to set the value of the coefficient of the rhs
+
+ff = @(x,y) 0;
+f = fem_coeff ('f', ff);
+
+gg = @(x,y) 0;
+g = fem_coeff ('g', gg);
+
+# fem_rhs and fem_lhs takes as input the functionspace V, and the
+# boundarycondition bc and solve the Poisson problem with
+# the velues specified inside f and g;
+t = 0;
+dt = 0.05;
+T = 3;
+
+kk = @(x,y) dt;
+k = fem_coeff ('k', kk);
+
+A = fem_rhs_Heat (V, bc, k);
+
+uu = @(x,y) 10*exp(-((x - 0.5).^2 + (y - 0.5).^2) / 0.02);
+u0 = fem_coeff ('u0', uu);
+
+# solve the problem for each time step
+# We need to update only the lhs
+while t < T
+  t += dt;
+  b = fem_lhs_Heat (V, f, g, k, u0, bc);
+
+  u = A \ b;
+
+  u0 = fem_func ('u0', V, u);
+  fem_plot (u0);
+#  fem_save (u0, num2str(t));
+end
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/example/test_laplace.m	Fri Jul 26 11:29:29 2013 +0200
@@ -0,0 +1,58 @@
+# We follow the dolfin example for the Poisson problem
+#    -div ( grad (u) ) = f      on omega
+#                    u = h      on gamma_d;
+#                du/dn = g      on gamma_n;
+# See (http://fenicsproject.org/documentation/dolfin/1.2.0/cpp/demo/pde/poisson/cpp/documentation.html#index-0)
+
+
+pkg load msh
+pkg load fem-fenics
+
+# initialize the environment and create function related to the problem
+# defined inside Laplace.ufl
+fem_init_env ();
+problem = 'Laplace';
+fem_create_fs (problem);
+fem_create_rhs (problem);
+fem_create_lhs (problem);
+
+# create a unit square mesh using msh: labels for the boundary sides are 1,2,3,4
+msho = msh2m_structured_mesh (0:0.05:1, 0:0.05:1, 1, 1:4);
+mshd = fem_init_mesh (msho);
+
+V  = fem_fs_Laplace (mshd);
+
+# fem_bc takes as input the functionspace V, a function handler f,
+# and the sides where we want to apply the condition
+# The value on each point of the boundary is computed using
+# the eval method available inside expression.h
+# if a side is not specified, Neumann conditions are applied
+# with g specified below
+f = @(x,y) 0;
+bc = fem_bc (V, f, [2, 4]);
+
+# fem_coeff takes as input a string and a function handler
+# and is used below to set the value of the coefficient of the rhs
+
+ff = @(x,y) 10*exp(-((x - 0.5).^2 + (y - 0.5).^2) / 0.02);
+f = fem_coeff ('f', ff);
+
+gg = @(x,y) sin (5.0 * x);
+g = fem_coeff ('g', gg);
+
+# fem_rhs and fem_lhs takes as input the functionspace V, and the
+# boundarycondition bc and solve the Lapalce problem with
+# the velues specified inside f and g;
+A = fem_rhs_Laplace (V, bc);
+
+b = fem_lhs_Laplace (V, f, g, bc);
+
+u = A \ b;
+
+# fem_func create a function object which can be plotted, saved or
+#  used to set the coefficient of the lhs/rhs
+func = fem_func ('u', V, u)
+
+fem_plot (func);
+
+fem_save (func, problem);
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/obsolete/fem_func_space.m	Fri Jul 26 11:29:29 2013 +0200
@@ -0,0 +1,75 @@
+## 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 2 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/>.
+
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {} = fem_ffc (myproblem.ufl)
+## This function take as input a .ufl file with a Variationa Problem
+## @example
+## 
+## @end example
+## @seealso{} 
+## @end deftypefn
+
+function res = fem_func_space (mesh, name)
+
+  if nargin != 2
+    error ("fem_ffc: wrong number of input parameters.");
+  elseif ! isstr (name)
+    error ("fem_ffc: first argument is not a valid string");
+  endif
+
+  # fem_fs.cc is the file where I should include name.h
+  # I use name.h at line 35, where I call
+  # "new FunctionSpace (mshd);" 
+  # because FunctionSpace is defined inside of name.h
+  # (correct call would be "name::FunctionSpace (mshd)" but 
+  # I'm using the namespace "name")
+  command = "cp ../src/fem_fs.cc fem_fs.cc";
+  [output, text] = system (command);
+  if output != 0
+    display (text);
+    error ("Compilation failed");
+  endif
+
+  # add line to fem_fs.cc
+  filename = "fem_fs.cc";
+  fid = fopen (filename, "r+");
+  fseek (fid, 671, SEEK_SET ());
+  fprintf (fid, "%c%s.h%c", 34, name, 34);
+  fseek (fid, 721, SEEK_SET ());
+  fprintf (fid, "%s%c", name, 59);
+  fclose (fid);
+
+  #now compile the program
+  command = "mkoctfile -c fem_fs.cc $(LDFLAGS) -o fem_fs.o -I../src/";
+  [output, text] = system (command);
+  if output != 0
+    display (text);
+    error ("Compilation failed");
+  endif
+
+  command = "mkoctfile -s ../src/mesh.o ../src/fem_init_env.o fem_fs.o -o fem_fs.oct $(LDFLAGS) -ldolfin";
+  [output, text] = system (command);
+  if output != 0
+    display (text);
+    error ("Compilation failed");
+  endif
+
+  #run the program and get the result
+  res = fem_fs (mesh);
+
+endfunction
+