FElupe documentation#
FElupe is a Python 3.8+ ๐ finite element analysis package ๐ฆ focussing on the formulation and numerical solution of nonlinear problems in continuum mechanics ๐ง of solid bodies ๐. Its name is a combination of FE (finite element) and the german word Lupe ๐ (magnifying glass) as a synonym for getting an insight ๐ how a finite element analysis code ๐งฎ looks like under the hood ๐ณ๏ธ.
Highlights
highlevel finiteelementanalysis API
flexible building blocks for finite element assembly
hyperelastic
integral (weak) forms
straightforward definition of
mixedfields
Installation#
Install Python, fire up ๐ฅ a terminal and run ๐
pip install felupe[all]
where [all]
installs all optional dependencies. FElupe has minimal requirements, all available at PyPI supporting all platforms.
numpy for array operations
scipy for sparse matrices
tensortrax for automatic differentiation
In order to make use of all features of FElupe ๐๐ฐ๐๐๐, it is suggested to install all optional dependencies.
Getting Started#
This tutorial covers the essential highlevel parts of creating and solving problems with FElupe. As an introductory example ๐จโ๐ซ, a quarter model of a solid cube
with hyperelastic material behaviour is subjected to a uniaxial()
elongation applied at a clamped endface.
First, letโs import FElupe and create a meshed cube
out of hexahedron
cells with a given number of points per axis. A numeric region
, predefined for hexahedrons, is created on the mesh. A vectorvalued displacement field
is initiated on the region. Next, a field container
is created on top of this field.
A uniaxial()
load case is applied on the displacement field
stored inside the field container
. This involves setting up symmetry()
planes as well as the absolute value of the prescribed displacement at the meshpoints on the rightend face of the cube. The rightend face is clamped ๐ ๏ธ: only displacements in direction x are allowed. The dict of boundary
conditions for this predefined load case are returned as boundaries
and the partitioned degrees of freedom as well as the external displacements are stored within the returned dict loadcase
.
An isotropic pseudoelastic OgdenRoxburgh
Mullinssoftening model formulation in combination with an isotropic hyperelastic NeoHookean
material formulation is applied on the displacement field
of a nearlyincompressible solid body
.
A step
generates the consecutive substepmovements of a given boundary
condition. The step
is further added to a list of steps of a job
๐ฉโ๐ป (here, a characteristic curve
๐ job is used). During evaluation
โณ, each substep of each step
is solved by an iterative NewtonRhapson
procedure โ๏ธ. The solution
is exported after each completed substep as a timeseries โ XDMF file. Finally, the result of the last completed substep is plotted.
import felupe as fem
mesh = fem.Cube(n=6)
region = fem.RegionHexahedron(mesh)
field = fem.FieldContainer([fem.Field(region, dim=3)])
boundaries, loadcase = fem.dof.uniaxial(field, clamped=True)
umat = fem.OgdenRoxburgh(material=fem.NeoHooke(mu=1), r=3, m=1, beta=0)
solid = fem.SolidBodyNearlyIncompressible(umat, field, bulk=5000)
move = fem.math.linsteps([0, 1, 0, 1, 2, 1], num=5)
step = fem.Step(items=[solid], ramp={boundaries["move"]: move}, boundaries=boundaries)
job = fem.CharacteristicCurve(steps=[step], boundary=boundaries["move"])
job.evaluate(filename="result.xdmf")
fig, ax = job.plot(
xlabel="Displacement $u$ in mm $\longrightarrow$",
ylabel="Normal Force $F$ in N $\longrightarrow$",
)
ax2 = solid.imshow("Principal Values of Cauchy Stress", theme="paraview")
Extension Packages#
The capabilities of FElupe may be enhanced with extension packages created by the community.
Package 
Description 

Constitutive hyperelastic material formulations 

Material Definition with Automatic Differentiation (AD) 

Math on HyperDual Tensors with Trailing Axes (bundled with FElupe) 

A visualization tool for FElupe 
Performance#
This is a simple benchmark to compare assembly times for linear elasticity and hyperelasticity on tetrahedrons.
Analysis 
DOF/s 

LinearElastic 
84523 +/ 9541 
Hyperelastic 
49771 +/ 6781 
Tested on: KDE Neon (Ubuntu 22.04), Python 3.10, Intelยฎ Coreโข i76650U CPU @ 2.20GHz, 8GB RAM (Surface Pro 4).
from timeit import timeit
import matplotlib.pyplot as plt
import numpy as np
import felupe as fem
def pre_linear_elastic(n, **kwargs):
mesh = fem.Cube(n=n).triangulate()
region = fem.RegionTetra(mesh)
field = fem.FieldContainer([fem.Field(region, dim=3)])
umat = fem.LinearElastic(E=1, nu=0.3)
solid = fem.SolidBody(umat, field)
return mesh, solid
def pre_hyperelastic(n, **kwargs):
mesh = fem.Cube(n=n).triangulate()
region = fem.RegionTetra(mesh)
field = fem.FieldContainer([fem.Field(region, dim=3)])
umat = fem.NeoHooke(mu=1.0, bulk=2.0)
solid = fem.SolidBody(umat, field)
return mesh, solid
print("# Assembly Runtimes")
print("")
print(" DOF  LinearElastic in s  Hyperelastic in s ")
print("      ")
points_per_axis = np.round((np.logspace(3, 5, 6) / 3)**(1 / 3)).astype(int)
number = 3
parallel = False
runtimes = np.zeros((len(points_per_axis), 2))
for i, n in enumerate(points_per_axis):
mesh, solid = pre_linear_elastic(n)
matrix = solid.assemble.matrix(parallel=parallel)
time_linear_elastic = (
timeit(lambda: solid.assemble.matrix(parallel=parallel), number=number) / number
)
mesh, solid = pre_hyperelastic(n)
matrix = solid.assemble.matrix(parallel=parallel)
time_hyperelastic = (
timeit(lambda: solid.assemble.matrix(parallel=parallel), number=number) / number
)
runtimes[i] = time_linear_elastic, time_hyperelastic
print(
f" {mesh.points.size:7d}  {runtimes[i][0]:19.2f}  {runtimes[i][1]:17.2f} "
)
dofs_le = points_per_axis ** 3 * 3 / runtimes[:, 0]
dofs_he = points_per_axis ** 3 * 3 / runtimes[:, 1]
print("")
print(" Analysis  DOF/s ")
print("    ")
print(
f" LinearElastic  {np.mean(dofs_le):5.0f} +/{np.std(dofs_le):5.0f} "
)
print(f" Hyperelastic  {np.mean(dofs_he):5.0f} +/{np.std(dofs_he):5.0f} ")
plt.figure()
plt.loglog(
points_per_axis ** 3 * 3,
runtimes[:, 1],
"C0",
label=r"Stiffness Matrix (Hyperelastic)",
)
plt.loglog(
points_per_axis ** 3 * 3,
runtimes[:, 0],
"C1",
label=r"Stiffness Matrix (LinearElastic)",
)
plt.xlabel(r"Number of degrees of freedom $\longrightarrow$")
plt.ylabel(r"Runtime in s $\longrightarrow$")
plt.legend()
plt.tight_layout()
plt.savefig("benchmark.png")
License#
FElupe  Finite Element Analysis (C) 20212024 Andreas Dutzler, Graz (Austria).
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 https://www.gnu.org/licenses/.