# HG changeset patch # User gedeone-octave # Date 1374830969 -7200 # Node ID 3231f888f25d73931e4bbcf339941e228437b69f # Parent a00f95c0b04868812f5e5edb15febb5e03b64b69 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 diff -r a00f95c0b048 -r 3231f888f25d example/Biharmonic.ufl --- /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 diff -r a00f95c0b048 -r 3231f888f25d example/Heat.ufl --- /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) diff -r a00f95c0b048 -r 3231f888f25d example/Laplace.ufl --- /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 diff -r a00f95c0b048 -r 3231f888f25d example/fem_func_space.m --- 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 . - - -## -*- 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 - diff -r a00f95c0b048 -r 3231f888f25d example/test_biharmonic.m --- /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); diff -r a00f95c0b048 -r 3231f888f25d example/test_heat.m --- /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 diff -r a00f95c0b048 -r 3231f888f25d example/test_laplace.m --- /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); diff -r a00f95c0b048 -r 3231f888f25d obsolete/fem_func_space.m --- /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 . + + +## -*- 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 +