gf_asm

Purpose

General assembly function.

Synopsis

vec F = gf_asm('volumic source', mesh_im mim, mesh_fem mf_u, mesh_fem mf_d, vec F)
vec F = gf_asm('boundary source',int boundary_num, mesh_im mim, mesh_fem mf_u, mesh_fem mf_d,vec G)
spmat M = gf_asm('mass matrix', mesh_im mim, mesh_fem mf1[, mesh_fem mf2])
spmat M = gf_asm('laplacian', mesh_im mim, mesh_fem mf_u, mesh_fem mf_d, vec A)
spmat K = gf_asm('linear elasticity', mesh_im mim, mesh_fem mf_u, mesh_fem mf_d, vec lambda_d, vec mu_d)
[spmat K,B] = gf_asm('stokes', mesh_im mim, mesh_fem mf_u, mesh_fem mf_p, mesh_fem mf_d, vec visc)
[spmat H,vec R] = gf_asm('dirichlet', int boundary_num, mesh_im mim, mesh_fem mf_u, 
mesh_fem mf_d, mat Hd, vec Rd)
M = gf_asm('boundary qu term', int boundary_num, mesh_im mim, mesh_fem mf_u, mesh_fem mf_d, mat Q)
[spmat Q, vec G,spmat H,vec R,vec F]=gf_asm('pdetool boundary conditions',
mesh_im mim, mesh_fem mf_u, mesh_fem mf_d, mat b, mat e[, string f_expr])
[...] = gf_asm('volumic'[, CVLST], string expr, mesh_im mim.., [mesh_fem mf1[, mf2,..]][,mat data...])
[...] = gf_asm('boundary', int bnum, string expr, mesh_im mim.., [mesh_fem mf1[, mf2,..]][,mat data...])
M = gf_asm('interpolation matrix', mesh_fem mf1, mesh_fem mf2)
M = gf_asm('extrapolation matrix', mesh_fem mf1, mesh_fem mf2)

Description

These assembly procedures all take an mf_u argument, which is the mesh_fem descriptor for the main unknown of the PDE. They usually take an mf_d argument, which describes the data FEM (i.e. Lamé coefficients for linear elasticity, fluid viscosity for stokes equation, etc...). Data mesh_fem are always expected to be scalar (i.e. Qdim==1), if they are used to describe a vector field V, then it is expected to have Q rows (and gf_mesh_fem_get(mf_d,'nbdof') columns). If you are not using exact integration methods, please make sure that the integration has a sufficiently high order (don't forget to take into account the degree of the geometrical transformation).

Fv=gf_asm('volumic source', mim, mf_u, mf_d, F) : assemble a volumic source term, on mf_u, using the data vector F defined on the data mesh_fem mf_d:

Fv = ∫Ωφi(x)F(x) dx with F(x)=∑j Fjψj(x)

Fb=gf_asm('boundary source', bnum, mim, mf_u, mf_d, F) is very similar, except that the integral is evaluated on the boundary bnum instead of the whole domain Ω.

M=gf_asm('mass matrix', mim, mf1 [, mf2]) : build the mass matrix

∫Ωφi(x).ψj(x) dx.

M=gf_asm('laplacian', mim, mf_u, mf_d, A) : do the assembly of elementary matrices for the Laplacian ∇.(a(x)∇u(x)):

∫a(x)(∇φu(x).∇φu(x)) with a(x)=∑Aiψi(x)

gf_asm('linear elasticity', mim, mf_u, mf_d, lambda_d, mu_d) : return the linear elasticity stiffness matrix: ∇.σ(x), where the stress tensor σ is σ(x)=Cijrsεrs and the strain tensor ε is εrs(u)=(∂rus+∂sur)/2 and Cijrs=λδijδrs + µ(δirδjs+δisδjr) (λ and µ are the Lamé coefficients). The mesh_fem mfu is expected to be such that gf_mesh_fem_get(mf_u,'Qdim') == gf_mesh_get(mf_u,'dim').

[K,B]=gf_asm('stokes', mim, mf_u, mf_p, mf_d, visc) : do the assembly of elementary matrices for the Stokes equation (viscous incompressible fluid) νdiv (ε(u))Δu - grad p=0, div u=0. On output, B is a sparse matrix corresponding to

∫Ωp(x).div v(x) dx
, and K is the linear elasticity stiffness matrix for λ=0 and 2µ=visc.

[H,R]=gf_asm('dirichlet', bnum, mim, mf_u, mf_d, Hd, Rd) : assemble conditions of type h(x).u(x) = r(x) where h is a small square matrix (of any rank) whose size is equal to gf_mesh_fem_get(mf_u,'Qdim'). This matrix is stored in Hd, one column per dof in mf_d, each column containing the values of the matrix h stored in Fortran order: for example Hd(:,j) = [h11(xj) h21(xj) h12(xj) h22(xj)] if u is a 2D vector field.

Of course, if the unknown u is a scalar field, Hd is just a row vector You may wonder why assembling Dirichlet conditions: these are usually expressed on a convenient mesh_fem (i.e. a Lagrangian one), while the mesh_fem of u might be more complex (i.e. non Lagrangian). Hence we need to project the constraints on the mesh_fem of u. This is basically identical to 
H = gf_asm('boundary qu term',bnum, mim, mf_u, mf_d, Hd);
R = gf_asm('boundary source',bnum, mim, mf_u, mf_d, Rd);
except that this function is smarter, in the sense that it tries to produce a "better" (more diagonal) constraints matrix HH (when possible): if it was not the case, H would be (in the general case) tridiagonal on 2D meshes when Hd is diagonal. Note that the rank of H still needs to be determined: [N,U0]=gf_spmat_get(H'dirichlet nullspace', R) does this. It solves the conditions HU=R, returning a solution U0 which has a minimum L2-norm. The sparse matrix N contains an orthogonal basis of the kernel of the constraints matrix H (hence, the PDE linear system should be solved on this subspace):
the initial problem
KU = B with constraints HU=R
is replaced by
(NKNT)V = N*B
U = N*V + U0

M=gf_asm('boundary qu term', boundary_num, mim, mf_u, mf_d, Q) : assemble the term ∫Γ (Q(x)φ(x)).ψ(x) dx where Q is a square matrix of size Qdim×Qdim, Qdim being the dimension of the unknown u (that is set when creating the mesh_fem). This is a kind of general boundary mass matrix.

The supplied argument Q should be a (Qdim2)×N array, where N is the number of degree of freedom of mf_d. Each column of Q contains the coefficients stored in the Fortran (and matlab order), for example if Qdim =2, Q(:,i) = q11, q21, q12, q22.

[Q,G,H,R,F]=gf_asm('pdetool boundary conditions', mim, mf_u, mf_d, b, e[, string f_expr]) is an easy way to assemble boundary conditions obtained from pdetool: b is the boundary matrix exported by the pdetool, and e is the edges array. f_expr is an optional expression (or vector) for the volumic term. On return Q,G,H,R,F contain the assembled boundary conditions (Q and H are matrices), similar to the ones returned by the function assemb from pdetool, and the solution U satisfies (K+Q)U=F+G under the constraints HU=R (K is the stiffness matrix of the PDE considered).

[...]=gf_asm({ 'volumic'[,CVLST] | 'boundary',bnum }, expr, mim1,..,[mf1[, mf2,..]][, data...]) is the generic assembly procedure for volumic and boundary assembly. The expression expr is evaluated over the mesh_fem listed in the arguments (with optional data) and assigned to the output arguments. For details about the syntax of assembly expressions, please refer to the getfem user manual (or look at the file getfem_assembling.h in the getfem++ sources).

For example, the L2 norm of a field can be computed with gf_compute(mf,U,'L2 norm') or with: 
gf_asm('volumic','u=data(#1); V()+=u(i).u(j).comp(Base(#1).Base(#1))(i,j)',mim,mf,U)
The Laplacian stiffness matrix can be evaluated with gf_asm('Laplacian',mim, mf, A) or equivalently with: 
gf_asm('volumic',['a=data(#2); ',...
       'M(#1,#1)+=sym(comp(Grad(#1).Grad(#1).Base(#2))(:,i,:,i,j).a(j))'], mim, mf, A);

gf_asm('interpolation matrix', mf1, mf2) : build the interpolation matrix from a mesh_fem onto another one (assumed to be Lagrangian). The returned sparse matrix M is such that V=M*U=gf_compute(mf1,U,'interpolate_on', mf1). This might be useful for repeated interpolations.

gf_asm('extrapolation matrix', mf1, mf2) is similar, but performs "light" extrapolation:,if some degrees of freedom of mf2 are slightly outside mf1, their value will be extrapolated from the values of the nearest D.o.F. of mf1.

See Also

gf_solve, gf_compute('interpolate on').