Assembly examples in Python

This is the use with Python interface of the C++ assembly examples in Compute arbitrary terms - high-level generic assembly procedures - Generic Weak-Form Language (GWFL).

[1]:
import getfem as gf
import numpy as np

elements_degree = 1
# elements_degree = 2

mesh = gf.Mesh("cartesian", np.arange(0.0, 2.0, 1.0))

mf is supposed to be an already declared gf.MeshFem object and mim a already declared gf.MeshIm object on the same mesh.

[2]:
mf = gf.MeshFem(mesh, 1)
mf.set_classical_fem(elements_degree)
print(mf)

BEGIN MESH_FEM

QDIM 1
 CONVEX 0 'FEM_PK(1,1)'
 BEGIN DOF_ENUMERATION
  0:  0 1
 END DOF_ENUMERATION
END MESH_FEM

[3]:
mim = gf.MeshIm(mesh, elements_degree*2)
print(mim)

BEGIN MESH_IM

 CONVEX 0 'IM_GAUSS1D(2)'
END MESH_IM

As a first example, if one needs to perform the assembly of a Poisson problem

\[-\mathrm{div} \nabla u=f,\mathrm{in} \Omega ,\]

the stiffness matrix is given

\[K_{i,j} =\int _{\Omega } \nabla \varphi _{i} \cdot \nabla \varphi _{j} dx,\]

and will be assembled by the following code:

[4]:
md = gf.Model("real")
md.add_fem_variable("u", mf)
md.add_nonlinear_term(mim, "Grad_u.Grad_Test_u")
md.assembly(option="build_matrix")
K = md.tangent_matrix()
print(K)
matrix(2, 2)
( (r0, 1) (r1, -1) )
( (r0, -1) (r1, 1) )

Note that the value of the variable do not really intervene because of the linearity of the problem. This allows to pass "u" as the value of the variable which will not be used. Note also that two other possible expressions for exactly the same result for the assembly string are "Grad_Test2_u.Grad_Test_u" (i.e. an order 2 expression) or "Norm_sqr(Grad_u)/2" (i.e. a potential). In fact other possible assembly string will give the same result such as "Grad_u.Grad_u/2" or "[Grad_u(1), Grad_u(2)].[Grad_Test_u(1), Grad_Test_u(2)]" for two-dimensional problems. However, the recommendation is preferably to give an order 1 expression (weak formulation) if there is no particular reason to prefer an order 0 or an order 2 expression.

[5]:
md = gf.Model("real")
md.add_fem_variable("u", mf)
md.add_nonlinear_term(mim, "Grad_u.Grad_u/2")
md.assembly(option="build_matrix")
K = md.tangent_matrix()
print(K)
matrix(2, 2)
( (r0, 1) (r1, -1) )
( (r0, -1) (r1, 1) )

[6]:
# for two-dimensional problems
# md = gf.Model("real")
# md.add_fem_variable("u", mf)
# md.add_nonlinear_term(mim, "[Grad_u(1), Grad_u(2)].[Grad_Test_u(1), Grad_Test_u(2)]")
# md.assembly(option="build_matrix")
# K = md.tangent_matrix()
# print(K)

As a second example, let us consider a coupled problem, for instance the mixed problem of incompressible elasticity given by the equations

\[\begin{split}-\mathrm{div}\left( \mu \left( \nabla u+( \nabla u)^{T}\right) -pI_{d}\right) =f,\mathrm{in} \thinspace \Omega , \\ -\mathrm{div} \thinspace u=0.\end{split}\]

where u is the vector valued displacement and p the pressure. The assembly of the matrix for the whole coupled system can be performed as follows:

[7]:
epsilon = 1.; E = 21E6; nu = 0.3;
clambda = E*nu/((1+nu)*(1-2*nu));
cmu = E/(2*(1+nu));

mf_u = gf.MeshFem(mesh, 1)
mf_u.set_classical_fem(elements_degree)
mf_p = gf.MeshFem(mesh, 1)
mf_p.set_classical_fem(elements_degree)

md = gf.Model("real")
md.add_fem_variable("u", mf_u)
md.add_fem_variable("p", mf_p)
md.add_initialized_data("mu", cmu)
md.add_nonlinear_term(mim, "2*mu*Sym(Grad_u):Grad_Test_u"
                     "- p*Trace(Grad_Test_u) - Test_p*Trace(Grad_u)")
md.assembly(option="build_matrix")
K = md.tangent_matrix()
print(K)
matrix(4, 4)
( (r2, 0.5) (r3, -0.5) )
( (r2, 0.5) (r3, -0.5) )
( (r0, 0.5) (r1, 0.5) (r2, 1.61538e+07) (r3, -1.61538e+07) )
( (r0, -0.5) (r1, -0.5) (r2, -1.61538e+07) (r3, 1.61538e+07) )

where, here, mf_u and mf_p are supposed to be some already declared getfem::mesh_fem objects defined on the same mesh, mim a already declared getfem::mesh_im object and mu is the Lame coefficient. It is also possible to perform the assembly of the sub-matrix of this system separately.

Let us see now how to perform the assembly of a source term. The weak formulation of a volumic source term is

Let us see now how to perform the assembly of a source term. The weak formulation of a volumic source term is

\[\int _{\Omega } fvdx\]

where \(f\) is the source term and \(v\) the test function. The corresponding assembly can be written:

[8]:
F = np.ones(mf_u.nbdof())

md = gf.Model("real")
md.add_fem_variable("u", mf_u)
md.add_initialized_fem_data("f", mf_u, F)
md.add_nonlinear_term(mim, "f*Test_u")
md.assembly("build_rhs")
rhs = md.rhs()
print(rhs)
[-0.5 -0.5]

if the source term is describe on a finite element mf_data and the corresponding vector of degrees of freedom F. Explicit source terms are also possible. For instance:

[9]:
region = -1 # ALL
md = gf.Model("real")
md.add_fem_variable("u", mf_u)
md.add_nonlinear_term(mim, "sin(X(1))*Test_u", region)
# for two-dimensional problems
# md.add_nonlinear_term(mim, "sin(X(1)+X(2))*Test_u", region)
md.assembly("build_rhs")
rhs = md.rhs()
print(rhs)
[-0.15767352 -0.30191429]

where region is the mesh region number.

As another example, let us describe a simple nonlinear elasticity problem. Assume that we consider a Saint-Venant Kirchhoff constitutive law which means that we consider the following elastic energy on a body of reference configuration \(\Omega\):

\[\int _{\Omega }\dfrac{\lambda }{2}\left(\mathrm{tr}( E)^{2} +\mu \mathrm{tr}\left( E^{2}\right)\right) dx\]

where \(\lambda\), \(\mu\) are the Lamé coefficients and \(E\) is the strain tensor given by \(E=\left( \nabla u+( \nabla u)^{T} +( \nabla u)^{T} \nabla u\right) /2\).

This is possible to perform the assembly of the corresponding tangent problem as follows:

[10]:
md = gf.Model("real")

md.add_fem_variable("u", mf_u)
md.add_initialized_data("lambda", clambda)
md.add_initialized_data("mu", cmu)
md.add_nonlinear_term(mim, "lambda*sqr(Trace(Grad_u+Grad_u'+Grad_u'*Grad_u))"
                      "+ mu*Trace((Grad_u+Grad_u'+Grad_u'*Grad_u)"
                      "*(Grad_u+Grad_u'+Grad_u'*Grad_u))")
md.assembly("build_rhs")
rhs = md.rhs()
print(rhs)
md.assembly("build_matrix")
print(K)
[0. 0.]
matrix(4, 4)
( (r2, 0.5) (r3, -0.5) )
( (r2, 0.5) (r3, -0.5) )
( (r0, 0.5) (r1, 0.5) (r2, 1.61538e+07) (r3, -1.61538e+07) )
( (r0, -0.5) (r1, -0.5) (r2, -1.61538e+07) (r3, 1.61538e+07) )

and to adapt a Newton-Raphson algorithm to solve that nonlinear problem. Of course the expression is rather repetitive and it would be preferable to define some intermediate nonlinear operators. However, note that repeated expressions are automatically detected and computed only once in the assembly.

The last example is the assembly of the stiffness matrix of an order four problem, the Kirchhoff-Love plate problem:

[11]:
h = 1.0 # mm
D = 2.0*h**3*E/(3*(1-nu**2))

md = gf.Model("real")

md.add_fem_variable("u", mf)
md.add_initialized_data("D", D)
md.add_initialized_data("nu", nu)
md.add_nonlinear_term(mim, "D*(1-nu)*(Hess_u:Hess_Test_u) - "
                           "D*nu*Trace(Hess_u)*Trace(Hess_Test_u)")
md.assembly(option="build_all")
K = md.tangent_matrix()

with \(D\) the flexion modulus and \(\nu\) the Poisson ratio.

[ ]: