--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/inst/femfenics_examples.m	Fri Sep 13 12:45:22 2013 +0200
@@ -0,0 +1,372 @@
+## 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/>.
+
+%# -*- texinfo -*-
+%# @deftypefn {Function File} {[@var{}] =} femfenics_examples (@var{})
+%# Open the PDE examples menu and allow the user to select a demo that will be evaluated.
+%# @end deftypefn
+
+function [] = femfenics_examples ()
+
+  vode = 1; while (vode > 0)
+    clc;
+    fprintf (1, ...
+      ['FEM-FENICS examples menu:\n', ...
+       '==================\n', ...
+       '\n', ...
+       '   (1) Solve the Poisson equation\n', ...
+       '   (2) Solve the Poisson equation with a Mixed formulation\n', ...
+       '   (3) Solve the Elasticity problem for an HyperElastic material\n', ...
+       '   (4) Solve the Navier-Stokes equation with Chorin-Teman algorithm\n', ...
+       '   (5) Solve the Biharmonic equation\n', ...
+       '\n', ...
+       '   Note: There are further FEM-FeNICS examples available with the pkg\n', ...
+       '         testsuite functions.\n', ...
+       '\n', ...
+       '   If you have another interesting example that you would like\n', ...
+       '   to share then please modify this file, create a patch and push\n', ...
+       '   your patch with your added example to the Fem-Fenics repository.\n', ...
+       '\n' ]);
+    vode = input ('Please choose a number from above or press <Enter> to return: ');
+    clc; if (vode > 0 && vode < 6)
+      %# We can't use the function 'demo' directly here because it does
+      %# not allow to run other functions within a demo.
+      vexa = example (mfilename (), vode);
+      disp (vexa);
+      eval (vexa);
+      input ('Press <Enter> to continue: ');
+    end %# if (vode > 0)
+  end %# while (vode > 0)
+
+%!demo
+%! # In this example the Poisson equation is solved using P1 element
+%! # on a uniform squared mesh.
+%! # Read about the Poisson equation at
+%! #   http://en.wikipedia.org/wiki/Poisson%27s_equation.
+%!
+%! n = length (mfilename ("fullpath")) - length (mfilename());
+%! path = strtrunc(mfilename ("fullpath"), n);
+%! private = fullfile (path, "example/Poisson/Poisson.ufl");
+%! [status, msg, msgid] = copyfile  (private, pwd, 'f');
+%! if (status != 1)
+%!   error ("couldn't copy file: %s", msg);
+%! endif
+%!
+%! pkg load fem-fenics msh
+%! import_ufl_Problem ('Poisson')
+%! 
+%! # Create mesh and define function space
+%! x = y = linspace (0, 1, 33);
+%! mesh = Mesh(msh2m_structured_mesh (x, y, 1, 1:4));
+%! 
+%! V = FunctionSpace('Poisson', mesh);
+%! 
+%! # Define boundary condition
+%! bc = DirichletBC(V, @(x, y) 0.0, [2;4]);
+%! 
+%! f = Expression ('f', @(x,y) 10*exp(-((x - 0.5)^2 + (y - 0.5)^2) / 0.02));
+%! g = Expression ('g', @(x,y) sin (5.0 * x));
+%! 
+%! a = BilinearForm ('Poisson', V);
+%! L = LinearForm ('Poisson', V, f, g);
+%! 
+%! # Compute solution
+%! [A, b] = assemble_system (a, L, bc);
+%! sol = A \ b;
+%! u = Function ('u', V, sol);
+%! 
+%! # Save solution in VTK format
+%! save (u, 'poisson')
+%! 
+%! # Plot solution
+%! plot (u);
+
+%!demo
+%! # In this example the Poisson equation is solved using a Mixed Approach.
+%! # We use the stable FE space obtained using Brezzi-Douglas-Marini polynomial 
+%! # of order 1 and Dicontinuos element of order 0.
+%! # Read about the Poisson equation at
+%! #   http://en.wikipedia.org/wiki/Poisson%27s_equation.
+%!
+%! n = length (mfilename ("fullpath")) - length (mfilename());
+%! path = strtrunc(mfilename ("fullpath"), n);
+%! private = fullfile (path, "example/Mixed-Poisson/MixedPoisson.ufl");
+%! [status, msg, msgid] = copyfile  (private, pwd, 'f');
+%! if (status != 1)
+%!   error ("couldn't copy file: %s", msg);
+%! endif
+%!
+%! pkg load fem-fenics msh
+%! import_ufl_Problem ('MixedPoisson')
+%! 
+%! # Create mesh
+%! x = y = linspace (0, 1, 33);
+%! mesh = Mesh(msh2m_structured_mesh (x, y, 1, 1:4));
+%! 
+%! V = FunctionSpace('MixedPoisson', mesh);
+%! 
+%! f = Expression ('f', @(x,y) 10*exp(-((x - 0.5)^2 + (y - 0.5)^2) / 0.02));
+%! 
+%! a = BilinearForm ('MixedPoisson', V);
+%! L = LinearForm ('MixedPoisson', V, f);
+%! 
+%! # Define essential boundary
+%! bc1 = DirichletBC (SubSpace (V, 1), @(x,y) [0; -sin(5.0*x)], 1);
+%! bc2 = DirichletBC (SubSpace (V, 1), @(x,y) [0;  sin(5.0*x)], 3);
+%! 
+%! # Compute solution
+%! [A, b] = assemble_system (a, L, bc1, bc2);
+%! sol = A \ b;
+%! func = Function ('func', V, sol);
+%! 
+%! sigma = Function ('sigma', func, 1);
+%! u = Function ('u', func, 2);
+%! 
+%! # Plot solution
+%! plot (sigma);
+%! plot (u);
+
+%!demo
+%! # In this example we solve the equation for an Hyperelastic material
+%! # minimizing the potential energy of the material.
+%! # Read about the Hyperelastic materials at
+%! #   http://en.wikipedia.org/wiki/Hyperelastic_material
+%!
+%! n = length (mfilename ("fullpath")) - length (mfilename());
+%! path = strtrunc(mfilename ("fullpath"), n);
+%! private = fullfile (path, "example/Elasticity/HyperElasticity.ufl");
+%! [status, msg, msgid] = copyfile  (private, pwd, 'f');
+%! if (status != 1)
+%!   error ("couldn't copy file: %s", msg);
+%! endif
+%!
+%! pkg load fem-fenics msh
+%! problem = 'HyperElasticity';
+%! import_ufl_Problem (problem);
+%! 
+%! Rotation = @(x,y,z) ...
+%!    [0; ...
+%!    0.5*(0.5 + (y - 0.5)*cos(pi/3) - (z-0.5)*sin(pi/3) - y);...
+%!    0.5*(0.5 + (y - 0.5)*sin(pi/3) + (z-0.5)*cos(pi/3) - z)];
+%! 
+%! #Create mesh and define function space
+%! x = y = z = linspace (0, 1, 17);
+%! mshd = Mesh (msh3m_structured_mesh (x, y, z, 1, 1:6));
+%! V  = FunctionSpace (problem, mshd);
+%! 
+%! # Create Dirichlet boundary conditions
+%! bcl = DirichletBC (V, @(x,y,z) [0; 0; 0], 1);
+%! bcr = DirichletBC (V, Rotation, 2);
+%! bcs = {bcl, bcr};
+%!  
+%! # Define source and boundary traction functions
+%! B = Constant ('B', [0.0; -0.5; 0.0]);
+%! T = Constant ('T', [0.1; 0.0; 0.0]);
+%! 
+%! # Set material parameters
+%! E = 10.0;
+%! nu = 0.3;
+%! mu = Constant ('mu', E./(2*(1+nu)));
+%! lmbda = Constant ('lmbda', E*nu./((1+nu)*(1-2*nu)));
+%! u = Expression ('u', @(x,y,z) [0; 0; 0]);
+%! 
+%! # Create (linear) form defining (nonlinear) variational problem
+%! L = ResidualForm (problem, V, mu, lmbda, B, T, u);
+%! 
+%! # Solve nonlinear variational problem F(u; v) = 0
+%! 
+%! # Create function for the resolution of the NL problem
+%!   function [y, jac] = f (problem, xx, V, bc1, bc2, B, T, mu, lmbda)
+%!     u = Function ('u', V, xx);
+%!     a = JacobianForm (problem, V, mu, lmbda, u);
+%!     L = ResidualForm (problem, V, mu, lmbda, B, T, u);
+%!     if (nargout == 1)
+%!       [y, xx] = assemble (L, xx, bc1, bc2);
+%!     elseif (nargout == 2)
+%!       [jac, y, xx] = assemble_system (a, L, xx, bc1, bc2);
+%!     endif
+%!   endfunction
+%! 
+%! u0 = assemble (L, bcs{:});
+%! fs = @(xx) f (problem, xx, V, bcl, bcr, B, T, mu, lmbda);
+%! [x, fval, info] = fsolve (fs, u0, optimset ("jacobian", "on"));
+%! func = Function ('u', V, x);
+%!  
+%! # Save solution in VTK format
+%! save (func, 'displacement');
+%! 
+%! # Plot solution
+%! plot (func);
+
+%!demo
+%! # In this example we solve the Navier-Stokes equation in a L-shape domain.
+%! # We impose a sinusoidal pressure at the inflow and a constant pressure at
+%! # the output. On the other sides, no-slip BC are applied.
+%! # We solve the problem with the Chorin-Temam algorithm.
+%! # The pressure is discretized with P1 element while the velocity with P2.
+%! # Read about the NS equation at
+%! #   http://en.wikipedia.org/wiki/Navier%E2%80%93Stokes_equations
+%!
+%! n = length (mfilename ("fullpath")) - length (mfilename());
+%! path = strtrunc(mfilename ("fullpath"), n);
+%! private = fullfile (path, "example/NavierStokes/*.ufl");
+%! [status, msg, msgid] = copyfile  (private, pwd, 'f');
+%! if (status != 1)
+%!   error ("couldn't copy file: %s", msg);
+%! endif
+%!
+%! pkg load fem-fenics msh
+%! import_ufl_Problem ("TentativeVelocity");
+%! import_ufl_Problem ("VelocityUpdate");
+%! import_ufl_Problem ("PressureUpdate");
+%!  
+%! # We use the msh pkg to generate the L-shape domain
+%! 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", .4);
+%! unlink (canonicalize_file_name (name));
+%!  
+%! mesh = Mesh (msho);
+%!  
+%! # Define function spaces (P2-P1). From ufl file
+%! V = FunctionSpace ('VelocityUpdate', mesh);
+%! Q = FunctionSpace ('PressureUpdate', mesh);
+%!  
+%! # Set parameter values
+%! dt = 0.01;
+%! T = 3.;
+%!  
+%! # Define boundary conditions
+%! noslip = DirichletBC (V, @(x,y) [0; 0], [3, 4]);
+%! outflow = DirichletBC (Q, @(x,y) 0, 2);
+%!  
+%! # Create functions
+%! u0 = Expression ('u0', @(x,y) [0; 0]);
+%!  
+%! # Define coefficients
+%! k = Constant ('k', dt);
+%! f = Constant ('f', [0; 0]);
+%!  
+%! # Tentative velocity step.
+%! a1 = BilinearForm ('TentativeVelocity', V, k);
+%!  
+%! # Pressure update.
+%! a2 = BilinearForm ('PressureUpdate', Q);
+%!  
+%! # Velocity update
+%! a3 = BilinearForm ('VelocityUpdate', V);
+%!  
+%! # Assemble matrices
+%! A1 = assemble (a1, noslip);
+%! A3 = assemble (a3, noslip);
+%!  
+%! # Time-stepping
+%! t = dt; i = 0;
+%! while t < T
+%!  
+%!   # Update pressure boundary condition
+%!   inflow = DirichletBC (Q, @(x,y) sin(3.0*t), 1);
+%!  
+%!   # Compute tentative velocity step
+%!   L1 = LinearForm ('TentativeVelocity', V, k, u0, f);
+%!   b1 = assemble (L1, noslip);
+%!   utmp = A1 \ b1;
+%!   u1 = Function ('u1', V, utmp);
+%!  
+%!   # Pressure correction
+%!   L2 = LinearForm ('PressureUpdate', Q, u1, k);
+%!   [A2, b2] = assemble_system (a2, L2, inflow, outflow);
+%!   ptmp = A2 \ b2;
+%!   p1 = Function ('p1', Q, ptmp);
+%!  
+%!   # Velocity correction
+%!   L3 = LinearForm ('VelocityUpdate', V, k, u1, p1);
+%!   b3 = assemble (L3, noslip);
+%!   ut = A3 \ b3;
+%!   u1 = Function ('u0', V, ut);
+%!  
+%!   # Plot solution
+%!   plot (p1);
+%!   plot (u1);
+%!  
+%!   # Save to file
+%!   save (p1, sprintf ("p_%3.3d", ++i));
+%!   delete (sprintf ("p_%3.3d.pvd", i));
+%!   save (u1, sprintf ("u_%3.3d", i));
+%!   delete (sprintf ("u_%3.3d.pvd", i));
+%!  
+%!   # Move to next time step
+%!   u0 = u1;
+%!   t += dt
+%!  
+%! end
+
+%!demo
+%! # In this example the Biharmonic equation is solved using P2 element
+%! # on a uniform squared mesh.
+%! # Read more about the Biharmonic equation at
+%! #   http://en.wikipedia.org/wiki/Biharmonic
+%!
+%! n = length (mfilename ("fullpath")) - length (mfilename());
+%! path = strtrunc(mfilename ("fullpath"), n);
+%! private = fullfile (path, "example/Biharmonic/Biharmonic.ufl");
+%! [status, msg, msgid] = copyfile  (private, pwd, 'f');
+%! if (status != 1)
+%!   error ("couldn't copy file: %s", msg);
+%! endif
+%! 
+%! pkg load fem-fenics msh
+%!  
+%! problem = 'Biharmonic';
+%! import_ufl_Problem (problem);
+%!  
+%! # Create mesh and define function space
+%! x = y = linspace (0, 1, 33);
+%! mesh = Mesh(msh2m_structured_mesh (x, y, 1, 1:4));
+%!  
+%! V = FunctionSpace(problem, mesh);
+%!  
+%! bc = DirichletBC (V, @(x,y) 0, 1:4);
+%!  
+%! f = Expression ('f', @(x,y) 4.0*pi^4.*sin(pi.*x).*sin(pi.*y));
+%! g = Expression ('alpha', @(x,y) 8);
+%!  
+%! a = BilinearForm (problem, V, g);
+%! L = LinearForm (problem, V, f);
+%!  
+%! [A, b] = assemble_system (a, L, bc);
+%! u = A \ b;
+%!  
+%! func = Function ('u', V, u);
+%! plot (func);
+
+%# Local Variables: ***
+%# mode: octave ***
+%# End: ***