Note

Go to the end to download the full example code.

Building Blocks#

Start setting up a problem in FElupe by the creation of a numeric Region with a geometry Mesh, a finite Element formulation Hexahedron and a Quadrature rule GaussLegendre.

import felupe as fem

mesh = fem.Cube(n=6)
element = fem.Hexahedron()
quadrature = fem.GaussLegendre(order=1, dim=3)

mesh.plot().show()

Static Scene

Interactive Scene

An overview of the mesh-metadata is available in the console.

print(mesh)

<felupe Mesh object>
  Number of points: 216
  Number of cells:
    hexahedron: 125

Region#

A Region essentially pre-calculates element shape/ansatz/basis functions and derivatives evaluated at every quadrature point of every cell w.r.t. the undeformed coordinates (as attribute dhdX).

Note

By using a template region like RegionHexahedron, only the mesh has to be created. The element formulation and the appropriate quadrature scheme are chosen automatically.

region = fem.Region(mesh, element, quadrature)
# region = fem.RegionHexahedron(mesh)
print(region)

<felupe Region object>
  Element formulation: Hexahedron
  Quadrature rule: GaussLegendre
  Gradient evaluated: True
  Hessian evaluated: False

The scheme of the region, i.e. the element with its point labels and the integration points of the quadrature rule may be plotted.

region.plot().show()

Static Scene

Interactive Scene

An array containing products of quadrature weights multiplied by the determinants of the (geometric) jacobians is stored as the array of (undeformed) differential volumes dV. The sum of all differential volumes gives the total (undeformed) volume of the region.

dV = region.dV
V = dV.sum()

Field#

In a next step, a displacement Field is added to the Region. This may be either a scalar- or a vector-valed field. The values at mesh-points are obtained with the attribute values. Interpolated field values at quadrature points are calculated with the interpolate() method. Additionally, the displacement gradient w.r.t. the undeformed coordinates is calculated for every quadrature point of every cell in the region with the field method grad(). A generalized extraction method extract(grad=True, add_identity=True, sym=False) allows several arguments to be passed. This involves whether the gradient or the values are extracted. If the gradient is extracted, the identity matrix may be added to the gradient (useful for the calculation of the deformation gradient). Optionally, the symmetric part is returned (small strain tensor).

displacement = fem.Field(region, dim=3)

u = displacement.values
v = displacement.interpolate()
dudX = displacement.grad()

Next, the field is added to a FieldContainer, which holds one or more fields. Like a Field, the FieldContainer also provides the extract(grad=True, add_identity=True, sym=False) method, returning a list of interpolated field values or gradients.

field = fem.FieldContainer([displacement])
print(field)

<felupe FieldContainer object>
  Number of fields: 1
  Dimension of fields:
    Field: 3

The deformation gradient is obtained by a sum of the identity and the displacement gradient. The extract()-method generates a list with extracted values at the integration points for each field in the field container. To evaluate the deformation gradient tensor of the first field of a container, the identity matrix is added to the gradient of the field.

F = field.extract(grad=True, sym=False, add_identity=True)

Methods to evaluate the deformation gradient as well as strain measures are provided in FieldContainer.evaluate.

import numpy as np

log_strain = field.evaluate.strain(fun=lambda stretch: np.log(stretch), tensor=True)
principal_stretches = field.evaluate.strain(fun=lambda stretch: stretch, tensor=False)

Constitution#

The material behavior has to be provided by the first Piola-Kirchhoff stress tensor as a function of the deformation gradient. FElupe provides a comprehensive constitutive library (see Constitution) including NeoHooke, LinearElastic and a generalized Hu-Washizu (u,p,J) ThreeFieldVariation. By alternative, an isotropic material formulation Hyperelastic is defined by a strain energy density function where both variation (stress) and linearization (elasticity) are carried out by automatic differentiation. The latter one is demonstrated here with a nearly-incompressible version of the Neo-Hookean material model formulation.

Note

FElupe supports different backends for automatic differentiation and uses tensortrax by default. It is important to use only differentiable math-functions from the backend, e.g. NumPy-like modules tensortrax.math, tensortrax.math.linalg or tensortrax.math.special. tensortrax is a dependency of FElupe and hence, is installed along with FElupe.

\[\psi = \frac{\mu}{2} \left( J^{-2/3}\ \text{tr}(\boldsymbol{C}) - 3 \right) + \frac{K}{2} \left( J - 1 \right)^2\]

import tensortrax.math as tm

import felupe.constitution.tensortrax as mat

def W(C, mu, bulk):
    "Isotropic hyperelastic Neo-Hookean material formulation."

    J = tm.sqrt(tm.linalg.det(C))

    return mu / 2 * (J ** (-2 / 3) * tm.trace(C) - 3) + bulk * (J - 1) ** 2 / 2

umat = mat.Hyperelastic(W, mu=1.0, bulk=2.0)

P = umat.gradient
A = umat.hessian

Boundary Conditions#

Next we enforce boundary conditions on the displacement field. Boundary conditions are stored as a dictionary of multiple Boundary instances. First, the left end of the cube is fixed. Displacements on the right end are fixed in directions y and z whereas displacements in direction x are prescribed with a user-defined value=0.5. A Boundary holds useful attributes like points or dof.

f0 = lambda x: np.isclose(x, 0)
f1 = lambda x: np.isclose(x, 1)

boundaries = fem.BoundaryDict()

boundaries["left"] = fem.Boundary(displacement, fx=f0)
boundaries["right"] = fem.Boundary(displacement, fx=f1, skip=(1, 0, 0))
boundaries["move"] = fem.Boundary(displacement, fx=f1, skip=(0, 1, 1), value=0.5)

boundaries.plot().show()

Static Scene

Interactive Scene

Partition of deegrees of freedom#

The separation of active and inactive degrees of freedom is performed by a so-called partition(). External values of prescribed displacement degrees of freedom are obtained by the application of the boundary values on the displacement field.

dof0, dof1 = fem.dof.partition(field, boundaries)
ext0 = fem.dof.apply(field, boundaries, dof0)

Integral forms of equilibrium equations#

The integral (or weak) forms of equilibrium equations are defined by the IntegralForm class, see Eq. (1). The pre-evaluated function of interest has to be passed as the fun argument and the test field as the v argument. By setting grad_v=[True] (default), FElupe passes the gradient of the first test field to the integral form. FElupe assumes a linear form if u=None (default) and creates a bilinear form if a field is passed to the trial field argument u.

(1)#\[\int_V P_{iJ} : \frac{\partial \delta u_i}{\partial X_J} \ dV \qquad \text{and} \qquad \int_V \frac{\partial \delta u_i} {\partial X_J} : \mathbb{A}_{iJkL} : \frac{\partial u_k}{\partial X_L} \ dV\]

linearform = fem.IntegralForm(P(F)[:-1], field, dV, grad_v=[True])
bilinearform = fem.IntegralForm(A(F), field, dV, u=field, grad_v=[True], grad_u=[True])

The assembly of both forms lead to the (point-based) internal force vector and the (sparse) stiffness matrix.

r = linearform.assemble()
K = bilinearform.assemble()

Prepare (partition) and solve the linearized equation system#

In order to solve the linearized equation system a partition() into active and inactive degrees of freedom has to be performed, see Eqs. (2). This system may then be passed to the (sparse direct) solver. Given a set of nonlinear equilibrium equations \(\boldsymbol{g}\) the unknowns \(\boldsymbol{u}\) are found by linearization at a valid initial state of equilibrium and an iterative Newton-Rhapson solution prodecure. The incremental values of inactive degrees of freedom are given as the difference of external prescribed and current values of unknowns. The (linear) solution is equal to the first result of a Newton-Rhapson iterative solution procedure. The solution du is finally added to the displacement field.

(2)#\[ \begin{align}\begin{aligned}\boldsymbol{g}_1(\boldsymbol{u}) &= -\boldsymbol{r}_1(\boldsymbol{u}) + \boldsymbol{f}_1\\\boldsymbol{g}_1(\boldsymbol{u} + d\boldsymbol{u}) &\approx -\boldsymbol{r}_1 + \boldsymbol{f}_1 - \frac{\partial \boldsymbol{r}_1}{\partial \boldsymbol{u}_1} \ d\boldsymbol{u}_1 - \frac{ \partial \boldsymbol{r}_1}{\partial \boldsymbol{u}_0} \ d\boldsymbol{u}_0 = \boldsymbol{0}\\d\boldsymbol{u}_0 &= \boldsymbol{u}_0^{(ext)} - \boldsymbol{u}_0\\\text{solve} \qquad \boldsymbol{K}_{11}\ d\boldsymbol{u}_1 &= \boldsymbol{g}_1 - \boldsymbol{K}_{10}\ d\boldsymbol{u}_{0}\\\boldsymbol{u}_0 &+= d\boldsymbol{u}_0\\\boldsymbol{u}_1 &+= d\boldsymbol{u}_1\end{aligned}\end{align} \]

The default solver of FElupe is SuperLU provided by the sparse package of SciPy. A significantly faster alternative is pypardiso which may be installed from PyPI with pip install pypardiso (not included with FElupe). The optional argument solver of solve() accepts any user-defined solver.

from scipy.sparse.linalg import spsolve  # default

# from pypardiso import spsolve

system = fem.solve.partition(field, K, dof1, dof0, r)
dfield = fem.solve.solve(*system, ext0, solver=spsolve)  # .reshape(*u.shape)
# field += dfield

A very simple Newton-Rhapson code looks like this:

for iteration in range(8):
    F = field.extract()

    linearform = fem.IntegralForm(P(F)[:-1], field, dV)
    bilinearform = fem.IntegralForm(A(F), field, dV, field)

    r = linearform.assemble()
    K = bilinearform.assemble()

    system = fem.solve.partition(field, K, dof1, dof0, r)
    dfield = fem.solve.solve(*system, ext0, solver=spsolve)

    norm = np.linalg.norm(dfield)
    print(iteration, norm)
    field += dfield

    if norm < 1e-12:
        break

4.525274987383307
0.16865632784122678
0.012102266767787884
6.444069678837376e-05
3.407185882229353e-09
3.6816610046338393e-16

By alternative, one may also use the Newton-Rhapson function of FElupe.

field[0].fill(0)
solid = fem.SolidBody(umat, field)
loadcase = {"dof1": dof1, "dof0": dof0, "ext0": ext0}
res = fem.newtonrhapson(items=[solid], verbose=2, tol=1e-12, **loadcase)
field = res.x

Newton-Rhapson solver
=====================

| # | norm(fun) |  norm(dx) |
|---|-----------|-----------|
| 1 | 3.188e-01 | 4.525e+00 |
| 2 | 2.719e-02 | 1.687e-01 |
| 3 | 1.758e-04 | 1.210e-02 |
| 4 | 6.337e-09 | 6.444e-05 |
| 5 | 9.403e-16 | 3.407e-09 |

Converged in 5 iterations (Assembly: 0.07556 s, Solve: 0.0147 s).

All 3x3 components of the deformation gradient of integration point 2 of cell 1 (Python is 0-indexed) are obtained with

defgrad = field.evaluate.deformation_gradient()
print(defgrad[:, :, 1, 0])

[[ 1.47073968 -0.029986   -0.029986  ]
 [ 0.23485302  0.90111657  0.00477373]
 [ 0.23485302  0.00477373  0.90111657]]

Export and plot of results#

Results are exported as VTK or XDMF files using meshio.

fem.save(region, field, filename="result.vtk")

Any tensor at quadrature points, shifted or projected to and averaged at mesh-points, is evaluated for quad and hexahedron cell types by topoints() or project(), respectively. For example, the calculation of the Cauchy stress involves the conversion from the first Piola-Kirchhoff stress to the Cauchy stress followed by the shift or the projection. The stress results at mesh points are passed as a dictionary to the point_data argument.

from felupe.math import det, dot, tovoigt, transpose

s = dot(P(F)[0], transpose(F[0])) / det(F[0])
# s = solid.evaluate.cauchy_stress()

fem.save(
    region,
    field,
    filename="result_with_cauchy.vtk",
    point_data={
        "CauchyStressShifted": fem.topoints(s, region),
        "CauchyStressProjected": fem.project(tovoigt(s), region),
    },
)

Instead of writing the simulation results to a VTK-file, point- and cell-data related to a solid body may also be viewed with PyVista. Once again, additional point-data like our stress results at mesh points are passed as a dictionary to the point_data argument of the view()-method. Cauchy stresses on solid bodies are also pre-defined in FElupe and may be plotted as cell- or (shifted / projected) point-data.

view = solid.view(
    point_data={
        "Cauchy Stress (Shifted)": fem.topoints(s, region),
        "Cauchy Stress (Projected)": fem.project(tovoigt(s), region),
    }
)

view.plot("Cauchy Stress (Shifted)", component=0).show()
solid.plot("Cauchy Stress", component=0).show()
# solid.plot("Cauchy Stress", component=0, project=fem.topoints).show()
# solid.plot("Cauchy Stress", component=0, project=fem.project).show()

Static Scene

Interactive Scene

Static Scene

Interactive Scene

Total running time of the script: (0 minutes 2.489 seconds)

Gallery generated by Sphinx-Gallery

Building Blocks#

Region#

Field#

Constitution#

Boundary Conditions#

Partition of deegrees of freedom#

Integral forms of equilibrium equations#

Prepare (partition) and solve the linearized equation system#

Export and plot of results#

This Page