gf_solve

Purpose

Solve PDEs. THIS FUNCTION IS DEPRECATED, USE THE MODEL BRICKS INSTEAD - the model bricks are much more powerful and fast, however they act as black-boxes, so for now the gf_solve function is left in the getfem-interface, for educational purposes.

Synopsis

U[,pde]=gf_solve(pde)

Description

The aim of this function is not to provide a general fast solver for all kinds of PDEs, but to serve as an example of use of the previous functions (especially assembly routines), and to provide an easy way to solve some basic PDEs. There are currently three PDEs handled by gf_solve:
  • the Laplacian: div (a(x)grad u(x)) + f = 0;
  • the linear elasticity: div σ(u) + f = 0, with σij=λεττ+2µεij;
  • and the Stokes equation: νΔu-∇p+f = 0.
The argument pde is a structure describing the PDE that is to be solved. The member pde.type can be 'laplacian', 'linear elasticity' or 'stokes'. The member pde.mf_u and pde.mf_d must be mesh_fem handles to the chosen mesh_fem (if the stokes solver is to be used, then one has to set also pde.mf_p, for the pressure). The coefficient (dependent of the PDE) must also be set. For the Laplacian, it is pde.lambda. For the linear elasticity, it is pde.lambda and pde.mu (Lamé coefficients). For Stokes, it is pde.viscos. These coefficients can be expressed in various forms: 
pde.viscos = { 1 };                        % constant coefficient
pde.viscos = { 'x.^2+y.^2' }                 % string expression
f=inline('x.^2+y.^2'); pde.viscos = { @f }   % function handle
pde.viscos = ones(1,gf_mesh_fem_get(1,pde.mf_d))  % dof values
Note the use of braces, this allows to express non-scalar coefficients with heterogeneous expressions such as { 1, 'x.*y' }. The volumic term must be set in pde.F. If the boundary condition were obtained from pdetool, one just has to set 
pde.pdetool.b = b;       
pde.pdetool.e = e;
and gf_solve will set the boundary numbers itself. For the general case, you will have to express the boundary conditions yourself, i.e. define the boundaries with gf_mesh_set(m,'boundary'), and fill the array pde.bound. For example, if the Qdim of mf_u is equal to 2, 
% Dirichlet condition HU=R on the boundary number 1
pde.bound(1).type = 'Dirichlet';
pde.bound(1).R = { 0, 0 };
pde.bound(1).H = { 1, 0; 0, 1 } % optional, if not set H will be eye(Qdim)

% Neumann condition (+optional boundary mass matrix) on the boundary number 2
pde.bound(1).type = 'Neumann';
pde.bound(1).G = { 0,0 };
pde.bound(1).Q = { 1, 0; 0, 1 } % optional, if not set Q will be zeros(Qdim)
    
% Mixed condition
pde.bound(1).type = 'Mixed';
pde.bound(1).R = { 'x', 0 };
pde.bound(1).H = { 1, 1; 0, 0 } % hence we impose $u_x+u_y=x$
pde.bound(1).G = { 0, 1 }
On output, the solution of the pde is returned, an the pde structure can also been returned (filled with assembled matrices and vectors in its field pde.asm). Note that if this structure is passed again as an argument to gf_solve, nothing will be computed, since gf_solve does the assembly of elements which are not found in the structure pde.asm. The solver itself is (for the moment) the slash operator of matlab (i.e. LU-factorization), except for the stokes problem where a conjugate gradient is used for the pressure.

Examples

Solving the stokes equation, using boundary condition and mesh from the pdetool: 
pde.type = 'stokes';
pde.viscos=1.0;
pde.pdetool.b = b;       % b and e were exported from the pdetool
pde.pdetool.e = e;
pde.F = { 0, 0 };   % volumic source term
m=gf_mesh('pt2D',p,t);   % mesh creation from the p and t arrays exported by  pdetool 
pde.mf_u=gf_mesh_fem(m,2);   % the displacement u is a vector field
pde.mf_p=gf_mesh_fem(m,1);   % the pressure is a scalar field
pde.mf_d=gf_mesh_fem(m,2); 
pde.mim=gf_mesh_im(m, gf_integ('IM_EXACT_SIMPLEX(2)'));
% we set the FEMs
gf_mesh_fem_set(pde.mf_u,'fem',gf_fem('FEM_PK(2,3)'));
gf_mesh_fem_set(pde.mf_d,'fem',gf_fem('FEM_PK(2,3)'));
gf_mesh_fem_set(pde.mf_p,'fem',gf_fem('FEM_PK_DISCONTINUOUS(2,1)'));

% and now we let the solver do its job
[U,P]=gf_solve(pde);

See Also

gf_asm, introduction Laplacian example *.