adtzlr / felupe

:mag: finite element analysis for continuum mechanics of solid bodies

Finite element analysis for continuum mechanics of solid bodies.

FElupe is a Python 3.8+ 🐍 finite element analysis package 📦 focusing 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 🕳️.

Installation

Install Python, fire up 🔥 a terminal and run 🏃

pip install felupe[all]

where [all] is a combination of [io,parallel,plot,progress,view] and 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.

• einsumt for parallel (threaded) assembly
• h5py for writing XDMF result files
• matplotlib for plotting graphs
• meshio for mesh-related I/O
• pyvista for interactive visualizations
• tqdm for showing progress bars during job evaluations

Getting Started

This tutorial covers the essential high-level 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 end-face.

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, pre-defined for hexahedrons, is created on the mesh. A vector-valued 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 mesh-points on the right-end face of the cube. The right-end face is clamped 🛠️: only displacements in direction x are allowed. The dict of boundary conditions for this pre-defined 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 pseudo-elastic Ogden-Roxburgh Mullins-softening model formulation in combination with an isotropic hyperelastic Neo-Hookean material formulation is applied on the displacement field of a nearly-incompressible solid body.

A step generates the consecutive substep-movements 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 Newton-Rhapson procedure ⚖️. The solution is exported after each completed substep as a time-series ⌚ XDMF file. Finally, the result of the last completed substep is plotted.

For more details beside this high-level code snippet, please have a look at the 📝 documentation.

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.1)
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$",
)

solid.plot("Principal Values of Cauchy Stress").show()

Documentation

The documentation is located here.

Extension Packages

The capabilities of FElupe may be enhanced with extension packages created by the community.

Package	Description
hyperelastic	Constitutive hyperelastic material formulations
matadi	Material Definition with Automatic Differentiation (AD)
tensortrax	Differentiable Tensors based on NumPy Arrays
feplot	A visualization tool for FElupe

Testing

To run the FElupe unit tests, check out this repository and type

tox

Scientific Publications

A list of articles in which is involved. If you use FElupe in your scientific work, please star this repository, cite it and add your publication to this list.

1. A. Dutzler, C. Buzzi, and M. Leitner, "Nondimensional translational characteristics of elastomer components", Journal of Applied Engineering Design and Simulation, vol. 1, no. 1. UiTM Press, Universiti Teknologi MARA, Sep. 21, 2021. doi: 10.24191/jaeds.v1i1.20.

2. C. Buzzi, A. Dutzler, T. Faethe, J. Lassacher, M. Leitner, and F.-J. Weber, "Development of a tool for estimating the characteristic curves of rubber-metal parts", in Proceedings of the 12th International Conference on Railway Bogies and Running Gears, 2022 (ISBN 978-963-9058-46-0).

3. J. Torggler, A. Dutzler, B. Oberdorfer, T. Faethe, H. Müller, C. Buzzi, and M. Leitner, "Investigating Damage Mechanisms in Cord-Rubber Composite Air Spring Bellows of Rail Vehicles and Representative Specimen Design", Applied Composite Materials. Springer Science and Business Media LLC, Aug. 22, 2023. doi: 10.1007/s10443-023-10157-1. Simulation-related Python scripts are available on GitHub at adtzlr/fiberreinforcedrubber.

Changelog

All notable changes to this project will be documented in this file. The format is based on Keep a Changelog, and this project adheres to Semantic Versioning.

License

FElupe - finite element analysis (C) 2021-2024 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/.

felupe's Issues

IntegralFormMixed

numpy upper triangle indices are hardcoded for three-field-formulations. Change this to the number of passed fields...

self.i, self.j = np.triu_indices(3)

[v1.0.1] Tune import time

Guess: if numba is installed, felupe compiles the parallel version of the integrate method of a form on import - even if it is not used. This should be handled in a better way in the future.

tests

Felupe's core code base has now evolved in a way that it should not change too much in the future. Now it would be a good time to add tests and evaluate coverage...

Wrong basisprime calculation in Hex20 elements

in fe.element.Hex2s() the basis function derivative is wrong. This leads to an error.

wrong solution of NewtonRhapson

it seems that u[dof0] - u0ext is not updated...

Improve `tools.newtonrhapson`

Make this function actually usable for both single- and mixed-field problems.

• add default fun
• add default jac
• add default solve
• add default check
• add tests

[v1.0.1] Add inverse method for GaussLegendre quadrature

Add a new inv() method for GaussLegendre quadrature schemes.

add general purpose Newton-Rhapson template function

As the Newton-Rhapson algorithm is frequently used, it would be a great addition to felupe to have a flexible template function. It should be able to

• handle both

• single and
• mixed field formulations

• let the user pass a function and its jacobian along with a given state of equilibrium
• both the input variables (unknowns) and the output (equilibrium residuals) should be able to be pre- and postprocess-able
• the result should be stored in a Result class, similar to scipy.optimize

[v1.0.1] Wrong stress projection to points

Quadrature points are indexed differently compared to the cell connectivity (due to the usage of quadpy). Thus, stress projection delivers wrong results.

A collection of tasks:

• For linear quadratures the point ordering should be permuted in order to match cell ordering.
• add a tools.project function for the stress projection to mesh-points (basically the code-snippet above)
• argument to average results
• add tools.project() to Getting Started Example (see tests)
• check axisymmetric stress result
• add permute for Gauss-Legendre order=2 dim=3

Fix tetrahedral meshes

If a tet mesh is imported with node numbering resulting in negative volumes this leads to an error in felupe.

Tasks:

• generate a fix function which can be applied to both triangle and tetrahedral meshes
• add quadrature rules for linear/quadratic triangle/tetrahedrons

Rename `mesh.as_discontinous` to `mesh.disconnect`

Better express what this function does. It disconnects all elements.

Access field values by slice

A field cannot be accessed with a slice. However, this could sometimes be convenient. Take the following code:

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

# this is the current way to obtain "u1"
u1 = displacement.values.ravel()[dof1]

# this should be possible...
u1 = displacement[dof1]

[v1.0.1] Improve Documentation

Better explain the calculation of cauchy stresses and the usage of pypardiso in the getting started example. Add this example to the test suite.

add reference loadcase "uniaxial"

as felupe.doftools.uniaxial(field, right, move, clamped=True)

and use it as

boundaries, dof0, dof1, u0ext = fe.doftools.uniaxial(displacement, right=1.0, move=0.4, clamped=False)

Add mesh conversion: insert edge-based midnodes

...in order to convert a mesh for the usage with quadratic elements (Quad8, Tria6, Tetra10 and Hexahedron20)

Error in fe.constitution.MaterialFrom

It seems that C, F, J are not passed to the subsequent materials AsIsochoric, etc.

Fix utils.axito3d

• quadrature.order was removed and hence, the order for the 3d quadrature can't be evaluated.
• element class names are not available anymore (Quad1 -> Quad, Hex1 -> Hexahedron)

Add Triangle and Tetrahedron element (cell) classes

As discussed in here,

• add support for linear & quadratic Triangle (Tri1, Tri2) and Tetrahedron (Tet1, Tet2) elements

progress:

• Tri1
• Tri2
• Tet1
• Tet2

Status: experimental

Improve LinearElastic material

modify LinearElastic(E, nu).elasticity(strain) to LinearElastic(E, nu).elasticity(strain, stress=None) in order to calculate the initial stress contribution of the elasticity if the stress argument is passed.

Remove `tovoigt2()` function in `math.py`

I think this function is not used in felupe.

Update field values

Field values can be updated with field += dxbut not withfield = field + dx. This is because Field.add()returnsNone.`

Improve Element Namings

should be similar to vtk cell types...

• Line, Triangle, Quad, Tetrahedron and Hexahedron for the linear elements
• QuadraticTriangle, QuadraticTetrahedron and QuadraticTetrahedron for quadratic elements with edge-based midpoints
• TriQuadraticHexahedron for quadratic lagrange-type hexahedron with mid-faces and mid-volume points
• ArbitraryOrderLagrange for "VTK_LAGRANGE_QUAD" or "VTK_LAGRANGE_HEXAHEDRON"

Simplify Usage of IntegralForm

During the design of FElupe, the intention was to create a single function for both LinearForm and BilinearForm, called IntegralForm. Unfortunately, this introduced the separate creation of IntegralFormMixed, as things are slightly different here. The idea is to unify IntegralForm and IntegralFormMixed into a single IntegralFormGeneralized class.

IntegralForm(fun, v, dV, grad_v=False, u=None, grad_u=False)

Also, the grad_v argument should be placed in front of u, so that one can code a Lineaform as

IntegralForm(fun, v, dV, True)

...No, keep it as it is.

For mixed formulations, fun is a list whereas for single-field formulations fun is casted into a length-one list. The same applies for the test- and trial-fields v and u. If v is a list or a tuple, then grad_v must also be of the same length, i.e. grad_v=(True, False, False). Inside this new IntegralFormGeneralized the code is 💯% the same for both single and mixed field formulations.

Region Templates

The idea is to introduce Region templates, i.e. RegionHexahedron, RegionQuad, etc.

These template Regions have pre-defined elements and quadratures and are shorter to use.

AxisymmetricField: region.cells does not exist anymore

change to region.mesh.cells

felupe.mesh.sweep re-sorts nodes

Due to the usage of np.unique the nodes are getting sorted. There should be at least an option to disable the sorting (which should then be the default mode). As Numpy does not provide an option this will require some work. Pandas provide a unique function too but returns no inverse index array.

Try/except statement on modules using sparse

In order that FElupe works without all the optional dependencies it is necessary to wrap import sparse into a try/except statement.

felupe/felupe/dof/__init__.py

Lines 13 to 16 in 3bf3a15

	from ._multipoint import (
	MultiPointConstraint,
	MultiPointContact,
	)

and here

felupe/felupe/dof/_multipoint.py

Line 30 in 3bf3a15

import sparse

Axisymmetry in combination with mixed-field formulations

Try to implement mixed-field formulations for axisymmetric problems.

First Task would be to merge IntegralFormAxisymmetric into IntegralForm. This should be done by adding a kind="axi" attribute to FieldAxisymmetric and then perform the adopted assemblage.

• add kind attribute to Field
• check kind of field in IntegralForm and adopt integrate and assemble methods
  - [ ] remove IntegralFormAxisymmetric class

Add `strain(field)` function in felupe.tools

calculate symmetric part of displacement gradient in voigt notation

• dim = 2-> 3 strain components
• dim = 3 -> 6 strain components

[v1.0.1] Wrong tovoigt() method in math

There is a typo in converting stress (3,3) to stress (6)

ij = [(0, 0), (1, 1), (2, 2), (0, 0), (1, 2), (0, 2)]

should be

ij = [(0, 0), (1, 1), (2, 2), (0, 1), (1, 2), (0, 2)]

instead.

[v1.0.1] Drop `quadpy` dependency

As mentioned in #92 quadpy increases the import time drastically. We can drop that dependency and use built-in functions of numpy.

Define slow?

Hi @adtzlr

Thanks for open-sourcing felupe. Looks great. I have used both FEniCS and scikit-fem and completely agree with you on native-installation of FEniCS on Windows (I use it mostly via Docker) and the the convenience of scikit-fem. For the latter, with numba even though the forms are considerably longer, the timings have been reasonable so far, especially for a pure python code.

I will try felupe for myself sometime this week. But have you observed significant differences between scikit-fem and felupe timings? Is this documented somewhere? My use cases also reside in finite incompressible elasticity so it might be worthwhile knowing it.

In any case, yet another nice tool in pure python to have at disposal.

Cheers!

_Form files are missing

Region: Can't invert `dXdr` if `mesh.ndim != element.ndim`

Easy fix: don't calculate the inverse to allow the creation of such Regions.

Improve tools.py

Several tasks have to be done here. Currently tools.py (newer) and utils.py (older) both contain useful tools/utilities - some are outdated, some need to be tweaked, etc. In the end, only one file should be kept.

Tasks

• improve tools.py
• merge utils.py -> tools.py
• delete utils.py

Error in `IntegralForm.assemble()` if `Field.ndim` not equal `mesh.ndim`

E.g. if a one-dimensional field is applied on a 2d quad-mesh, there are Errors due to dimension mismatch.

Remove unused `mathsymmetric.py`

This file was created to seperate some linear algebra functions for symmetrc matrices. However, these (simpler) calculations are now in math.py, where some functions have an optional sym argument.

Naming convention / Nomenclature

Go through all files and unify namings.

Overview of tasks...

• mesh.py (see 810f82f)
• doftools.py (see 875d9e8)
• element.py (see 66552d4)
• field.py (see 6fa911b)
• forms.py (see 716ce21)
  - [x] math.py
  - [x] quadrature.py
• region.py (see 90a8838)
  - [x] solve.py
• utils.py (see b32ee30)

Document (docstrings) and comment code

Add docstrings and comments to code.

Overview

• doftools.py
• element.py
• field.py
• forms.py
• math.py
• quadrature.py
• solve.py
• utils.py
• region.py

Simplify mesh.py

Mesh generators are too complicated for now. There should be only 2d generators (rectangles, disks). Cubes or cylinders (with a hole) should only be expansions or revolutions from 2d sections. The scaled versions of Rectangle and Cube should be created with a scale function.

Implement Multi point constraints

as RBE2, Penalty Method, see here

update docs

according to #19 some wordings have to be adopted.

Error if phi < 360 in mesh.revolve()

must have happened in #12

Fix Tet10 element definition

Tet10 element implementation is currently broken

Extend default loadcases: provide uniaxal, biaxial and planar loadings

Also set clamped=False because this is the real homogenous loadcase. Of course, for these loading situations a finite element analysis doesn't make too much sense. They just serve as a reference implementation of some simple loadcases.

Unify attributes

Across felupe classes, several attributes have different names, e.g. ndim and dim. This issue serves as a collection of them.

Dimensions

attribute	description
dim, ndim	dimension of ...
npoints, ncells	number of ...

PyPardiso replaced by SuperLU - PyPardiso can't be installed via pip

PyPardiso was removed in 3b4806d because it can't be installed via pip. (dependancy mkl-service is missing, see haasad/PyPardiso#24).

If someone would like to use pypardiso, install it the recommended way and specify the solver argument in felupe.solve.solve.

import pypardiso

# user code

felupe.solve.solve(*system, solver=pypardiso.spsolve)

[v1.0.1] add mesh-method `discontinous`

Introduce a method for transforming a connected mesh into a disconnected mesh. Template code:

import numpy as np
from felupe import Mesh


def as_discontinous(self):
    points = self.points[self.cells].reshape(-1, self.dim)
    cells = np.arange(self.cells.size).reshape(*self.cells.shape)
    return Mesh(points, cells, cell_type=self.cell_type)

Gradient and Hessian of AD-materials for repeated values of principal stretches

casADi's implementation of eig_symbolic returns nan if a matrix contains repeated eigenvalues.

Simplify poisson example

The current poisson example should be simplified to a one-dimensional field. As this was not possible in the past, the example is based on a 2d field. However, this is possible since #37 and should be adopted now.

• add a laplace-helper function to the math module which takes a field as argument felupe.math.laplace(field)