Simple finite element assemblers
scikit-fem
is a pure Python 3.8+ library for performing finite element
assembly. Its main
purpose is the transformation of bilinear forms into sparse matrices and linear
forms into vectors.
The library
If you use the library in your research, you can cite the following article:
@article{skfem2020,
doi = {10.21105/joss.02369},
year = {2020},
volume = {5},
number = {52},
pages = {2369},
author = {Tom Gustafsson and G. D. McBain},
title = {scikit-fem: A {P}ython package for finite element assembly},
journal = {Journal of Open Source Software}
}
The most recent release can be installed simply by
pip install scikit-fem[all]
Remove [all]
to not install the optional dependencies meshio
for mesh
input/output, and matplotlib
for creating simple visualizations.
The minimal dependencies are numpy
and scipy
.
You can also try the library in browser through Google Colab.
Solve the Poisson problem (see also ex01.py
):
from skfem import *
from skfem.helpers import dot, grad
# create the mesh
mesh = MeshTri().refined(4)
# or, with your own points and elements:
# mesh = MeshTri(points, elements)
basis = Basis(mesh, ElementTriP1())
@BilinearForm
def laplace(u, v, _):
return dot(grad(u), grad(v))
@LinearForm
def rhs(v, _):
return 1. * v
A = laplace.assemble(basis)
b = rhs.assemble(basis)
# Dirichlet boundary conditions
A, b = enforce(A, b, D=mesh.boundary_nodes())
# solve the linear system
x = solve(A, b)
# plot using matplotlib
mesh.plot(x, shading='gouraud', colorbar=True).show()
# or, save to external file:
# mesh.save('output.vtk', point_data={'solution': x})
Meshes can be initialized manually, loaded from external files using meshio, or created with the help of special constructors:
import numpy as np
from skfem import MeshLine, MeshTri, MeshTet
mesh = MeshLine(np.array([0., .5, 1.]))
mesh = MeshTri(
np.array([[0., 0.],
[1., 0.],
[0., 1.]]).T,
np.array([[0, 1, 2]]).T,
)
mesh = MeshTri.load("docs/examples/meshes/square.msh") # requires meshio
mesh = MeshTet.init_tensor(*((np.linspace(0, 1, 60),) * 3))
We support many common finite elements. Below the stiffness matrix is assembled using second-order tetrahedra:
from skfem import Basis, ElementTetP2
basis = Basis(mesh, ElementTetP2()) # quadratic tetrahedron
A = laplace.assemble(basis) # type: scipy.sparse.csr_matrix
More examples can be found in the gallery.
The following benchmark (docs/examples/performance.py
) demonstrates the time
spent on finite element assembly in comparison to the time spent on linear
solve. The given numbers were calculated using a ThinkPad X1 Carbon laptop (7th
gen). Note that the timings are only illustrative as they depend on, e.g., the
type of element used, the number of quadrature points used, the type of linear
solver, and the complexity of the forms. This benchmark solves the Laplace
equation using linear tetrahedral elements and the default direct sparse solver
of scipy.sparse.linalg.spsolve
.
Degrees-of-freedom | Assembly (s) | Linear solve (s) |
---|---|---|
4096 | 0.04805 | 0.04241 |
8000 | 0.09804 | 0.16269 |
15625 | 0.20347 | 0.87741 |
32768 | 0.46399 | 5.98163 |
64000 | 1.00143 | 36.47855 |
125000 | 2.05274 | nan |
262144 | 4.48825 | nan |
512000 | 8.82814 | nan |
1030301 | 18.25461 | nan |
The project is documented using Sphinx under docs/
.
Built version can be found from Read the Docs.
Here are direct links to additional resources:
If you encounter an issue you can use GitHub issue tracker. If you cannot find help from the documentation, you can use the GitHub Discussions to ask questions. Try to provide a snippet of code which fails and include also the version of the library you are using. The version can be found as follows:
import skfem; print(skfem.__version__)
The minimal dependencies for installing scikit-fem
are
numpy and scipy. In addition,
many
examples use
matplotlib for visualization and
meshio for loading/saving different mesh
file formats. Some examples demonstrate the use of other external packages;
see requirements.txt
for a list of test dependencies.
The tests are run by GitHub Actions. The Makefile
in the repository root has
targets for running the testing container locally using docker
. For example,
make test_py38
runs the tests using py38
branch from
kinnala/scikit-fem-docker-action.
The releases are tested in
kinnala/scikit-fem-release-tests.
The contents of skfem/
and the PyPI package scikit-fem
are licensed under
the 3-clause BSD license. Some examples under docs/examples/
or snippets
in the documentation may have a different license.
This project was started while working under a grant from the Finnish Cultural Foundation. Versions 2.0.0+ were prepared while working in a project funded by Academy of Finland (decisions nr. 324611 and 338341). The approach used in the finite element assembly has been inspired by the work of A. Hannukainen and M. Juntunen.
We are happy to welcome any contributions to the library. Reasonable projects for first timers include:
By contributing code to scikit-fem, you are agreeing to release it under BSD-3-Clause, see LICENSE.md.
The format is based on Keep a Changelog, and this project adheres to Semantic Versioning with respect to documented and/or tested features.
numpy==2.0rc1
MappingAffine
now uses lazy evaluation also for element
mappings, in addition to boundary mappingsMeshTet.init_tensor
uses significantly less memory for
large meshesMesh.load
uses less memory when loading and matching tagsBasis
has new optional disable_doflocs
kwarg which
can be set to True
to avoid computing Basis.doflocs
, to reduce memory
usageElementVector
works also for split_bases/split_indices in case mesh.dim() != elem.dim
MappingMortar
and MortarFacetBasis
in favor of
skfem.supermeshing
skfem.visuals.glvis
; current version is broken and no
replacement is being plannedMesh.load
supports new keyword arguments
ignore_orientation=True
and ignore_interior_facets=True
which
will both speed up the loading of larger three-dimensional meshes by
ignoring facet orientation and all tags not on the boundary,
respectively.skfem.supermeshing
(requires shapely>=2
) for creating quadrature
rules for interpolating between two 1D or 2D meshes.Mesh.remove_unused_nodes
Mesh.remove_duplicate_nodes
Mesh.remove_elements
now supports passing any subdomain
reference through Mesh.normalize_elements
; subdomains and
boundaries are also properly preservedMeshTet
uniform refine was reindexing subdomains incorrectlyMeshDG.draw
did not work; now calls Basis.draw
which
works for any mesh topologyFacetBasis
now works with MeshTri2
, MeshQuad2
,
MeshTet2
and MeshHex2
ElementGlobal
now uses outward normals to initialize DOFs
on boundary facetsElementTriHHJ0
and ElementTriHHJ1
matrix finite elements
for implementing the Hellan-Hermann-Johnson mixed method (see ex48.py
)ElementHexSkeleton0
, piecewise constant element on the skeleton
of a hexahedral meshElementHexC1
, C1-continuous hexahedral elementMesh.restrict
now preserves subdomains and boundariesskfem.utils.mpc
for setting multipoint constraintsBasis.get_dofs
now supports fetching DOFs at specific nodes
through kwarg nodes
DofsView.sort
for sorting a set of DOFs returned by
Basis.get_dofs
CellBasis.with_elements
as a shortcut for initiating a new Basis with subset of elementsMesh.refined
now preserves subdomains in adaptive mesh refinementMesh.refined
used to modify the element connectivity of the original meshBasis.find_dofs
method, see Basis.get_dofs
for a
replacementElementTetN0
to ElementTriN1
and added alias for backwards
compatibilityElementQuadN1
, first order H(curl) conforming quadrilateral elementElementTriN1
, first order H(curl) conforming triangle elementElementTriN2
, second order H(curl) conforming triangle elementElementTetSkeletonP0
, extension of ElementTriSkeletonP0
to
tetrahedral meshesMesh.restrict
which returns a new mesh given a subset of elements or subdomainMesh.trace
which turns facets into a trace meshskfem.utils.bmat
, a variant of scipy.sparse.bmat
which adds the
indices of the different blocks as an attribute out.blocks
skfem.models.elasticity
Basis.interpolate
and Basis.project
now support specifying dtype
and using complex fieldsBasis.intepolate
did not work properly with ElementComposite
when the basis was defined only for a subset of elementsBasis.split
worked incorrectly for ElementVector
and multiple DOFs
of same typeElementQuadP
basis for reused quadrature points
did not work correctlyMappingMortar
and MortarFacetBasis
in favor of
plain FacetBasis
and the upcoming skfem.experimental.supermeshing
(see ex04.py
)DiscreteField.is_zero
in the
helpers to skip the evaluation of zero components in ElementComposite
to
improve type stability with respect to the size of the underlying numpy
arrays; this is technically a backwards incompatible change and might affect
self-created helper functionsFacetBasis.trace
in favor of Basis.interpolator
and Basis.project
Basis.interpolator
supports trailing axes; can be now
passed to Basis.project
for (inexact) interpolation between meshesElementTriRT0
to ElementTriRT1
and added alias for
backwards compatibilityElementTetRT0
to ElementTetRT1
and added alias for
backwards compatibilityElementQuadRT0
to ElementQuadRT1
and added alias for
backwards compatibilityElementTriRT2
, the second order Raviart-Thomas elementElementHexRT1
, the first order Raviart-Thomas element for hexahedral meshesBasis.project
now better supports ElementComposite
solver_iter_cg
, a simple pure Python conjugate gradient solver for
environments that do not have sparse solver libraries (e.g., Pyodide)ElementTriP2B
and ElementTriP1B
, new aliases for ElementTriMini
and ElementTriCCR
ElementTriP1G
and ElementTriP2G
, variants of ElementTriP1
and
ElementTriP2
using ElementGlobal
so that second derivatives are available
(useful, e.g., for stabilized methods and the Stokes problem)Basis.plot3
, a wrapper to skfem.visuals.*.plot3
Basis.__repr__
was slow and incorrectElementHdiv
did not work together with FacetBasis
DiscreteField
is now a subclass of ndarray
instead of
NamedTuple
and, consequently, the components of DiscreteField
cannot no
more be indexed inside forms like u[1]
(use u.grad
instead)w['u']
and w.u
inside the form definition is now
equivalent (previously w.u == w['u'].value
)Mesh.draw
now uses matplotlib
by default, calling
Mesh.draw("vedo")
uses vedo
skfem.visuals.matplotlib
now uses jet
as the default colormapBoundaryFacetBasis
is now an alias of FacetBasis
instead of
other way aroundDiscreteField.value
remains for backwards-compatibility but is
now deprecated and can be droppedMesh.plot
, a wrapper to skfem.visuals.*.plot
Basis.plot
, a wrapper to skfem.visuals.*.plot
Basis.refinterp
now supports vectorial fieldsskfem.visuals.matplotlib.plot
now has a basic quiver plot for vector
fieldsMesh.facets_around
which constructs a set of facets around a
subdomainMesh.save
and load
now preserve the orientation of boundaries and
interfacesOrientedBoundary
which is a subclass of ndarray
for facet index
arrays with the orientation information (0 or 1 per facet) available as
OrientedBoundary.ori
FacetBasis
will use the facet orientations (if present) to calculate
traces and normal vectorsskfem.visuals.matplotlib.draw
will visualize the orientations if
boundaries=True
is givenMesh.facets_satisfying
allows specifying the keyword argument
normal
for orienting the resulting interfaceFacetBasis
constructor now has the keyword argument side
which
allows changing the side of the facet used to calculate the basis function
values and gradientsBasis.boundary
method to quickly initialize the corresponding
FacetBasis
asm
/assemble
keyword
argumentsCellBasis
did not properly support elements
argumentBasis.interpolate
did not properly interpolate all components of
ElementComposite
Basis.project
, a more general and easy to use alternative for
projection
Basis
and FacetBasis
kwargs elements
and facets
can now be a
string refering to subdomain and boundary tagsElementQuadRT0
, lowest-order quadrilateral Raviart-Thomas elementFunctional
returned only the first component for forms with
non-scalar outputskfem.helpers.mul
for matrix multiplicationBasis.split
will now also split ElementVector
into its componentsElementDG
was not included in the wildcard importMeshTri2
and MeshQuad2
in Jupyter
notebooks raised exceptionmeshio
is now an optional dependencyElementComposite
uses DiscreteField()
to represent zeroBasis.get_dofs
skfem.__version__
Mesh.boundaries
during uniform refinement of MeshLine1
,
MeshTri1
and MeshQuad1
MeshTet1
that failed to produce
a valid meshBasis.find_dofs
in favor of Basis.get_dofs
DofsView
objects via +
and |
because of surprising
behavior in some edge casesMappingIsoparametric
can now be pickledMesh.save
/Mesh.load
now exports/imports Mesh.subdomains
and
Mesh.boundaries
Mesh.load
now optionally writes any mesh data to a list passed via
the keyword argument out
, e.g., out=data
where data = ['point_data']
Mesh.load
(and skfem.io.meshio.from_file
) now supports the
additional keyword argument force_meshio_type
for loading mesh files that
have multiple element types written in the same file, one element type at
a timeasm
will now accept a list of bases, assemble the same form using
all of the bases and sum the result (useful for jump terms and mixed meshes,
see ex41.py
Mesh.with_boundaries
now allows the definition of internal boundaries/interfaces
via the flag boundaries_only=False
MeshTri1DG
, MeshQuad1DG
, MeshHex1DG
, MeshLine1DG
; new mesh
types for describing meshes with a discontinuous topology, e.g., periodic
meshes (see ex42.py
)ElementHexDG
for transforming hexahedral H1 elements to DG/L2 elements.ElementTriP1DG
, ElementQuad1DG
, ElementHex1DG
,
ElementLineP1DG
; shorthands for ElementTriDG(ElementTriP1())
etc.ElementTriSkeletonP0
and ElementTriSkeletonP1
for defining
Lagrange multipliers on the skeleton mesh (see ex40.py
)TrilinearForm
for assembling a sparse 3-tensor, e.g., when dealing
with unknown material dataMeshTri.oriented
for CCW oriented triangular meshes which can be
useful for debugging or interfacing to external toolsMeshWedge1
and ElementWedge1
, the lowest order
wedge mesh and elementElementTriP3
, cubic triangular Lagrange elementElementTriP4
, quartic triangular Lagrange elementElementTri15ParamPlate
, 15-parameter nonconforming triangular element for platesElementTriBDM1
, the lowest order Brezzi-Douglas-Marini elementMesh.draw().show()
will now visualize any mesh interactively (requires vedo)MeshTet1
MappingIsoparametric
is now about 2x faster for large meshes thanks
to additional cachingMeshHex2.save
did not work properlyMesh.load
ignores unparseable cell_sets
inserted by meshio
in MSH 4.1Mesh
string representation is now more informativeForm.assemble
no more allows keyword arguments with list
or
dict
type: from now on only DiscreteField
or 1d/2d ndarray
objects are
allowed and 1d ndarray
is passed automatically to Basis.interpolate
for
convenienceMeshLine
is now a function which initializes MeshLine1
and not an alias to MeshLine1
FacetBasis
is now a shorthand for BoundaryFacetBasis
and no
longer initializes InteriorFacetBasis
or MortarFacetBasis
if the keyword
argument side
is passed to the constructorMesh.define_boundary
methodElementTriCCR
and ElementTetCCR
, conforming Crouzeix-Raviart finite elementsMesh.mirrored
returned a wrong mesh when a point other than the origin was usedMeshLine
constructor accepted only numpy arrays and not plain Python listsMesh.element_finder
(and CellBasis.probes
, CellBasis.interpolator
) was not working properly for a small number of elements (<5) or a large number of input points (>1000)MeshTet
and MeshTri.element_finder
are now more robust against degenerate elementsMesh.element_finder
(and CellBasis.probes
, CellBasis.interpolator
) raises exception if the query point is outside of the domainBasis
, a shorthand for CellBasis
CellBasis
, a new preferred name for InteriorBasis
BoundaryFacetBasis
, a new preferred name for ExteriorFacetBasis
skfem.utils.penalize
, an alternative to condense
and enforce
for
essential boundary conditionsInteriorBasis.point_source
, with ex38
ElementTetDG
, similar to ElementTriDG
for tetrahedral meshesMeshLine1.element_finder
Mesh
base class which is "immutable" and uses
Element
classes to define the ordering of nodes; better support for
high-order and other more general mesh types in the futureMeshTri2
, MeshQuad2
, MeshTet2
and MeshHex2
InteriorBasis.probes
; like InteriorBasis.interpolator
but returns a matrix
that operates on solution vectors to interpolate them at the given pointsDiscreteField
, e.g., multiplication, summation
and subtraction are now explicitly supported inside the form definitionsMeshHex.to_meshtet
for splitting hexahedra into tetrahedraMeshHex.element_finder
for interpolating finite element solutions
on hexahedral meshes via InteriorBasis.interpolator
Mesh.with_boundaries
, a functional replacement to
Mesh.define_boundary
, i.e. defining boundaries via Boolean lambda functionMesh.with_subdomains
for defining subdomains via Boolean lambda functionskfem.utils.projection
, a replacement of skfem.utils.project
with a different, more intuitive order of argumentsskfem.utils.enforce
for setting essential boundary conditions by
changing matrix rows to zero and diagonals to one.skfem.utils.project
in favor of skfem.utils.projection
Mesh.define_boundary
in favor of Mesh.with_boundaries
Mesh.{refine,scale,translate}
; the replacements are Mesh.{refined,scaled,translated}
skfem.models.helpers
; available as skfem.helpers
DiscreteField.{f,df,ddf,hod}
; available as DiscreteField.{value,grad,hess,grad3,...}
skfem.utils.L2_projection
skfem.utils.derivative
Mesh.refined
no more attempts to fix the indexing of Mesh.boundaries
after refineskfem.utils.solve
now uses scipy.sparse.eigs
instead of scipy.sparse.eigsh
by default;
the old behavior can be retained by explicitly passing solver=solver_scipy_eigs_sym()
skfem.visuals.matplotlib
related to 1D plottingside
keyword argument to FacetBasis
in favor of the more
explicit InteriorFacetBasis
and MortarFacetBasis
.InteriorFacetBasis
for integrating over the interior facets, e.g.,
evaluating error estimators with jumps and implementing DG methods.MortarFacetBasis
for integrating over the mortar mesh.InteriorBasis.with_element
for reinitializing an equivalent basis
that uses a different element.Form.partial
for applying functools.partial
to the form function
wrapped by Form
.asm
.Mesh2D.mirror
in favor of the more general Mesh.mirrored
.Mesh.refine
, Mesh.scale
and Mesh.translate
in favor of
Mesh.refined
, Mesh.scaled
and Mesh.translated
.Mesh.refined
, Mesh.scaled
, and Mesh.translated
. The new methods
return a copy instead of modifying self
.Mesh.mirrored
for mirroring a mesh using a normal and a point.Functional
now supports forms that evaluate to vectors or other
tensors.ElementHex0
, piecewise constant element for hexahedral meshes.FacetBasis.trace
for restricting existing solutions to lower
dimensional meshes on boundaries or interfaces.MeshLine.refined
now correctly performs adaptive refinement of
one-dimensional meshes.ElementLineP0
, one-dimensional piecewise constant element.skfem.helpers.curl
now calculates the rotated gradient for
two-dimensional elements.MeshTet.init_ball
for meshing a ball.ElementQuad0
was not compatible with FacetBasis
.TestEx32
more robust.tests
from the PyPI distribution.L2_projection
will be replaced by project
.derivative
will be replaced by project
.MeshTet.element_finder
and MeshLine.element_finder
for using
InteriorBasis.interpolator
.ElementTriCR
, the nonconforming Crouzeix-Raviart element for Stokes flow.ElementTetCR
, tetrahedral nonconforming Crouzeix-Raviart element.ElementTriHermite
, an extension of ElementLineHermite
to triangular
meshes.Mesh.validate
for unsigned Mesh.t
.Mesh3D.boundary_edges
: tested to run on a laptop
with over 10 million elements.ElementHex2
, a triquadratic hexahedral element.MeshTri.init_circle
, constructor for a circle mesh.Mesh3D.boundary_edges
(and, consequently, Basis.find_dofs
) was slow
and used lots of memory due to an exhaustive search of all edges.project
will only support functions like lambda x: x[0]
instead of lambda x, y, z: x
in the future.BilinearForm
and LinearForm
now take
an optional argument dtype
which defaults to np.float64
but can be also np.complex64
.Dofs.__or__
and Dofs.__add__
, for merging degree-of-freedom sets
(i.e. Dofs
objects) using |
and +
operators.Dofs.drop
and Dofs.keep
, for further filtering the degree-of-freedom setsbilinear_form
, linear_form
, and
functional
(deprecated since 1.0.0).FacetBasis
did not initialize with ElementQuadP
.MeshQuad._splitquads
aliased as MeshQuad.to_meshtri
.Mesh.__add__
, for merging meshes using +
operator: duplicated nodes are
joined.ElementHexS2
, a 20-node quadratic hexahedral serendipity element.ElementLineMini
, MINI-element for one-dimensional mesh.Mesh3D.boundary_edges
was broken in case of hexahedral meshes.skfem.utils.project
did not work for ElementGlobal
.ElementTetMini
, MINI-element for tetrahedral mesh.Mesh3D.boundary_edges
incorrectly returned all edges where both nodes are on
the boundary.bilinear_form
, linear_form
, and functional
.Basis.interpolate
returns DiscreteField
objects instead of ndarray tuples.Basis.interpolate
works now properly for vectorial and high-order elements
by interpolating all components and higher order derivatives.Form.assemble
accepts now any keyword arguments (with type DiscreteField
)
that are passed over to the forms.skfem.importers
to skfem.io
.skfem.models.helpers
to skfem.helpers
.skfem.utils.solve
will now expand also the solutions of eigenvalue problems.BilinearForm
, LinearForm
, and Functional
.skfem.io.json
for serialization of meshes to/from json-files.ElementLinePp
, p-th order one-dimensional elements.ElementQuadP
, p-th order quadrilateral elements.ElementQuadDG
for transforming quadrilateral H1 elements to DG elements.ElementQuadBFS
, Bogner-Fox-Schmit element for biharmonic problems.ElementTriMini
, MINI-element for Stokes problems.ElementComposite
for using multiple elements in one bilinear form.ElementQuadS2
, quadratic Serendipity element.ElementLineHermite
, cubic Hermite element for Euler-Bernoulli beams.Mesh.define_boundary
for defining named boundaries.Basis.find_dofs
for finding degree-of-freedom indices.Mesh.from_basis
for defining high-order meshes.Basis.split
for splitting multicomponent solutions.MortarMapping
with basic support for mortar methods in 2D.Basis
constructors now accept quadrature
keyword argument for specifying
a custom quadrature rule.