## 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, 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, 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, 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, V, k);
%!  
%! # Pressure update.
%! a2 = BilinearForm ('PressureUpdate', Q, Q);
%!  
%! # Velocity update
%! a3 = BilinearForm ('VelocityUpdate', V, 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, 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: ***