Simple finite element assemblers
Project description
scikitfem
scikitfem
is a pure Python 3.7+ 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
 has minimal dependencies
 contains no compiled code
 supports onedimensional, triangular, quadrilateral, tetrahedral and hexahedral finite elements
 includes special elements such as RaviartThomas, Nédélec, MINI, CrouzeixRaviart, Argyris, ...
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 = {scikitfem: A {P}ython package for finite element assembly},
journal = {Journal of Open Source Software}
}
Installation
The most recent release can be installed simply by
pip install scikitfem[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.
Examples
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 secondorder 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.
Benchmark
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
.
Degreesoffreedom  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 
Documentation
The project is documented using Sphinx under docs/
.
Built version can be found from Read the Docs.
Here are direct links to additional resources:
Getting help
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__)
Dependencies
The minimal dependencies for installing scikitfem
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.
Testing
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/scikitfemdockeraction.
The releases are tested in
kinnala/scikitfemreleasetests.
Licensing
The contents of skfem/
and the PyPI package scikitfem
are licensed under
the 3clause BSD license. Some examples under docs/examples/
or snippets
in the documentation may have a different license.
Acknowledgements
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.
Contributing
We are happy to welcome any contributions to the library. Reasonable projects for first timers include:
By contributing code to scikitfem, you are agreeing to release it under BSD3Clause, see LICENSE.md.
Changelog
The format is based on Keep a Changelog, and this project adheres to Semantic Versioning with respect to documented and/or tested features.
[7.0.1]  20220803
 Fixed: Updated changelog was missing.
[7.0.0]  20220803
 Changed: Removed the optimization of using
DiscreteField.is_zero
in the helpers to skip the evaluation of zero components inElementComposite
to improve type stability with respect to the size of the underlying numpy arrays; this is technically a backwards incompatible change and might affect selfcreated helper functions  Deprecated:
FacetBasis.trace
in favor ofBasis.interpolator
andBasis.project
 Added: Output of
Basis.interpolator
supports trailing axes; can be now passed toBasis.project
for (inexact) interpolation between meshes  Added: Renamed
ElementTriRT0
toElementTriRT1
and added alias for backwards compatibility  Added: Renamed
ElementTetRT0
toElementTetRT1
and added alias for backwards compatibility  Added: Renamed
ElementQuadRT0
toElementQuadRT1
and added alias for backwards compatibility  Added:
ElementTriRT2
, the second order RaviartThomas element  Added:
ElementHexRT1
, the first order RaviartThomas element for hexahedral meshes  Added:
Basis.project
now better supportsElementComposite
[6.0.0]  20220315  Added:
solver_iter_cg
, a simple pure Python conjugate gradient solver for environments that do not have sparse solver libraries (e.g., Pyodide)  Added:
ElementTriP2B
andElementTriP1B
, new aliases forElementTriMini
andElementTriCCR
 Added:
ElementTriP1G
andElementTriP2G
, variants ofElementTriP1
andElementTriP2
usingElementGlobal
so that second derivatives are available (useful, e.g., for stabilized methods and the Stokes problem)  Added:
Basis.plot3
, a wrapper toskfem.visuals.*.plot3
 Fixed: Calculation of size in
Basis.__repr__
was slow and incorrect  Fixed: Subclasses of
ElementHdiv
did not work together withFacetBasis
[6.0.0]  20220315
 Changed:
DiscreteField
is now a subclass ofndarray
instead ofNamedTuple
and, consequently, the components ofDiscreteField
cannot no more be indexed inside forms likeu[1]
(useu.grad
instead)  Changed: Writing
w['u']
andw.u
inside the form definition is now equivalent (previouslyw.u == w['u'].value
)  Changed:
Mesh.draw
now usesmatplotlib
by default, callingMesh.draw("vedo")
usesvedo
 Changed:
skfem.visuals.matplotlib
now usesjet
as the default colormap  Changed:
BoundaryFacetBasis
is now an alias ofFacetBasis
instead of other way around  Deprecated:
DiscreteField.value
remains for backwardscompatibility but is now deprecated and can be dropped  Added:
Mesh.plot
, a wrapper toskfem.visuals.*.plot
 Added:
Basis.plot
, a wrapper toskfem.visuals.*.plot
 Added:
Basis.refinterp
now supports vectorial fields  Added:
skfem.visuals.matplotlib.plot
now has a basic quiver plot for vector fields  Added:
Mesh.facets_around
which constructs a set of facets around a subdomain  Added:
Mesh.save
andload
now preserve the orientation of boundaries and interfaces  Added:
OrientedBoundary
which is a subclass ofndarray
for facet index arrays with the orientation information (0 or 1 per facet) available asOrientedBoundary.ori
 Added:
FacetBasis
will use the facet orientations (if present) to calculate traces and normal vectors  Added:
skfem.visuals.matplotlib.draw
will visualize the orientations ifboundaries=True
is given  Added:
Mesh.facets_satisfying
allows specifying the keyword argumentnormal
for orienting the resulting interface  Added:
FacetBasis
constructor now has the keyword argumentside
which allows changing the side of the facet used to calculate the basis function values and gradients  Added:
Basis.boundary
method to quickly initialize the correspondingFacetBasis
 Fixed: Improvements to backwards compatibility in
asm
/assemble
keyword arguments  Fixed: Save format issue with meshio 5.3.0+
 Fixed:
CellBasis
did not properly supportelements
argument  Fixed:
Basis.interpolate
did not properly interpolate all components ofElementComposite
[5.2.0]  20211227
 Added:
Basis.project
, a more general and easy to use alternative forprojection
 Added:
Basis
andFacetBasis
kwargselements
andfacets
can now be a string refering to subdomain and boundary tags  Added:
ElementQuadRT0
, lowestorder quadrilateral RaviartThomas element  Fixed:
Functional
returned only the first component for forms with nonscalar output
[5.1.0]  20211130
 Added:
skfem.helpers.mul
for matrix multiplication  Added:
Basis.split
will now also splitElementVector
into its components  Fixed:
ElementDG
was not included in the wildcard import  Fixed: Automatic visualization of
MeshTri2
andMeshQuad2
in Jupyter notebooks raised exception
[5.0.0]  20211121
 Changed:
meshio
is now an optional dependency  Changed:
ElementComposite
usesDiscreteField()
to represent zero  Added: Support more argument types in
Basis.get_dofs
 Added: Version information in
skfem.__version__
 Added: Preserve
Mesh.boundaries
during uniform refinement ofMeshLine1
,MeshTri1
andMeshQuad1
 Fixed: Refinement of quadratic meshes will now fall back to the refinement algorithm of firstorder meshes instead of crashing
 Fixed: Edge cases in the adaptive refine of
MeshTet1
that failed to produce a valid mesh  Deprecated:
Basis.find_dofs
in favor ofBasis.get_dofs
 Deprecated: Merging
DofsView
objects via+
and
because of surprising behavior in some edge cases
[4.0.1]  20211015
 Fixed:
MappingIsoparametric
can now be pickled
[4.0.0]  20210927
 Added:
Mesh.save
/Mesh.load
now exports/importsMesh.subdomains
andMesh.boundaries
 Added:
Mesh.load
now optionally writes any mesh data to a list passed via the keyword argumentout
, e.g.,out=data
wheredata = ['point_data']
 Added:
Mesh.load
(andskfem.io.meshio.from_file
) now supports the additional keyword argumentforce_meshio_type
for loading mesh files that have multiple element types written in the same file, one element type at a time  Added:
asm
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 Example 41)  Added:
Mesh.with_boundaries
now allows the definition of internal boundaries/interfaces via the flagboundaries_only=False
 Added:
MeshTri1DG
,MeshQuad1DG
,MeshHex1DG
,MeshLine1DG
; new mesh types for describing meshes with a discontinuous topology, e.g., periodic meshes (see Example 42)  Added:
ElementHexDG
for transforming hexahedral H1 elements to DG/L2 elements.  Added:
ElementTriP1DG
,ElementQuad1DG
,ElementHex1DG
,ElementLineP1DG
; shorthands forElementTriDG(ElementTriP1())
etc.  Added:
ElementTriSkeletonP0
andElementTriSkeletonP1
for defining Lagrange multipliers on the skeleton mesh (see Example 40)  Added:
TrilinearForm
for assembling a sparse 3tensor, e.g., when dealing with unknown material data  Added:
MeshTri.oriented
for CCW oriented triangular meshes which can be useful for debugging or interfacing to external tools  Added: partial support for
MeshWedge1
andElementWedge1
, the lowest order wedge mesh and element  Added:
ElementTriP3
, cubic triangular Lagrange element  Added:
ElementTriP4
, quartic triangular Lagrange element  Added:
ElementTri15ParamPlate
, 15parameter nonconforming triangular element for plates  Added:
ElementTriBDM1
, the lowest order BrezziDouglasMarini element  Added:
Mesh.draw().show()
will now visualize any mesh interactively (requires vedo)  Added: Adaptive refinement for
MeshTet1
 Fixed:
MappingIsoparametric
is now about 2x faster for large meshes thanks to additional caching  Fixed:
MeshHex2.save
did not work properly  Fixed:
Mesh.load
ignores unparseablecell_sets
inserted bymeshio
in MSH 4.1  Changed:
Mesh
string representation is now more informative  Changed:
Form.assemble
no more allows keyword arguments withlist
ordict
type: from now on onlyDiscreteField
or 1d/2dndarray
objects are allowed and 1dndarray
is passed automatically toBasis.interpolate
for convenience  Changed:
MeshLine
is now a function which initializesMeshLine1
and not an alias toMeshLine1
 Changed:
FacetBasis
is now a shorthand forBoundaryFacetBasis
and no longer initializesInteriorFacetBasis
orMortarFacetBasis
if the keyword argumentside
is passed to the constructor  Removed: the deprecated
Mesh.define_boundary
method
[3.2.0]  20210802
 Added:
ElementTriCCR
andElementTetCCR
, conforming CrouzeixRaviart finite elements  Fixed:
Mesh.mirrored
returned a wrong mesh when a point other than the origin was used  Fixed:
MeshLine
constructor accepted only numpy arrays and not plain Python lists  Fixed:
Mesh.element_finder
(andCellBasis.probes
,CellBasis.interpolator
) was not working properly for a small number of elements (<5) or a large number of input points (>1000)  Fixed:
MeshTet
andMeshTri.element_finder
are now more robust against degenerate elements  Fixed:
Mesh.element_finder
(andCellBasis.probes
,CellBasis.interpolator
) raises exception if the query point is outside of the domain
[3.1.0]  20210618
 Added:
Basis
, a shorthand forCellBasis
 Added:
CellBasis
, a new preferred name forInteriorBasis
 Added:
BoundaryFacetBasis
, a new preferred name forExteriorFacetBasis
 Added:
utils.penalize
, an alternative tocondense
andenforce
for essential boundary conditions  Added:
InteriorBasis.point_source
, withex38
 Added:
ElementTetDG
, similar toElementTriDG
for tetrahedral meshes  Fixed:
MeshLine1.element_finder
[3.0.0]  20210419
 Added: Completely rewritten
Mesh
base class which is "immutable" and usesElement
classes to define the ordering of nodes; better support for highorder and other more general mesh types in the future  Added: New quadratic mesh types:
MeshTri2
,MeshQuad2
,MeshTet2
andMeshHex2
 Added:
InteriorBasis.probes
; likeInteriorBasis.interpolator
but returns a matrix that operates on solution vectors to interpolate them at the given points  Added: More overloads for
DiscreteField
, e.g., multiplication, summation and subtraction are now explicitly supported inside the form definitions  Added:
MeshHex.to_meshtet
for splitting hexahedra into tetrahedra  Added:
MeshHex.element_finder
for interpolating finite element solutions on hexahedral meshes viaInteriorBasis.interpolator
 Added:
Mesh.with_boundaries
, a functional replacement toMesh.define_boundary
, i.e. defining boundaries via Boolean lambda function  Added:
Mesh.with_subdomains
for defining subdomains via Boolean lambda function  Added:
skfem.utils.projection
, a replacement ofskfem.utils.project
with a different, more intuitive order of arguments  Added:
skfem.utils.enforce
for setting essential boundary conditions by changing matrix rows to zero and diagonals to one.  Deprecated:
skfem.utils.project
in favor ofskfem.utils.projection
 Deprecated:
Mesh.define_boundary
in favor ofMesh.with_boundaries
 Removed:
Mesh.{refine,scale,translate}
; the replacements areMesh.{refined,scaled,translated}
 Removed:
skfem.models.helpers
; available asskfem.helpers
 Removed:
DiscreteField.{f,df,ddf,hod}
; available asDiscreteField.{value,grad,hess,grad3,...}
 Removed: Python 3.6 support
 Removed:
skfem.utils.L2_projection
 Removed:
skfem.utils.derivative
 Changed:
Mesh.refined
no more attempts to fix the indexing ofMesh.boundaries
after refine  Changed:
skfem.utils.solve
now usesscipy.sparse.eigs
instead ofscipy.sparse.eigsh
by default; the old behavior can be retained by explicitly passingsolver=solver_scipy_eigs_sym()
 Fixed: High memory usage in
skfem.visuals.matplotlib
related to 1D plotting
[2.5.0]  20210213
 Deprecated:
side
keyword argument toFacetBasis
in favor of the more explicitInteriorFacetBasis
andMortarFacetBasis
.  Added:
InteriorFacetBasis
for integrating over the interior facets, e.g., evaluating error estimators with jumps and implementing DG methods.  Added:
MortarFacetBasis
for integrating over the mortar mesh.  Added:
InteriorBasis.with_element
for reinitializing an equivalent basis that uses a different element.  Added:
Form.partial
for applyingfunctools.partial
to the form function wrapped byForm
.  Fixed: Include explicit Python 3.9 support.
[2.4.0]  20210120
 Deprecated: List and tuple keyword argument types to
asm
.  Deprecated:
Mesh2D.mirror
in favor of the more generalMesh.mirrored
.  Deprecated:
Mesh.refine
,Mesh.scale
andMesh.translate
in favor ofMesh.refined
,Mesh.scaled
andMesh.translated
.  Added:
Mesh.refined
,Mesh.scaled
, andMesh.translated
. The new methods return a copy instead of modifyingself
.  Added:
Mesh.mirrored
for mirroring a mesh using a normal and a point.  Added:
Functional
now supports forms that evaluate to vectors or other tensors.  Added:
ElementHex0
, piecewise constant element for hexahedral meshes.  Added:
FacetBasis.trace
for restricting existing solutions to lower dimensional meshes on boundaries or interfaces.  Fixed:
MeshLine.refined
now correctly performs adaptive refinement of onedimensional meshes.
[2.3.0]  20201124
 Added:
ElementLineP0
, onedimensional piecewise constant element.  Added:
skfem.helpers.curl
now calculates the rotated gradient for twodimensional elements.  Added:
MeshTet.init_ball
for meshing a ball.  Fixed:
ElementQuad0
was not compatible withFacetBasis
.
[2.2.3]  20201016
 Fixed: Remove an unnecessary dependency.
[2.2.2]  20201015
 Fixed: Make the preconditioner in
TestEx32
more robust.
[2.2.1]  20201015
 Fixed: Remove
tests
from the PyPI distribution.
[2.2.0]  20201014
 Deprecated:
L2_projection
will be replaced byproject
.  Deprecated:
derivative
will be replaced byproject
.  Added:
MeshTet.element_finder
andMeshLine.element_finder
for usingInteriorBasis.interpolator
.  Added:
ElementTriCR
, the nonconforming CrouzeixRaviart element for Stokes flow.  Added:
ElementTetCR
, tetrahedral nonconforming CrouzeixRaviart element.  Added:
ElementTriHermite
, an extension ofElementLineHermite
to triangular meshes.  Fixed: Fix
Mesh.validate
for unsignedMesh.t
.
[2.1.1]  20201001
 Fixed: Further optimizations to
Mesh3D.boundary_edges
: tested to run on a laptop with over 10 million elements.
[2.1.0]  20200930
 Added:
ElementHex2
, a triquadratic hexahedral element.  Added:
MeshTri.init_circle
, constructor for a circle mesh.  Fixed:
Mesh3D.boundary_edges
(and, consequently,Basis.find_dofs
) was slow and used lots of memory due to an exhaustive search of all edges.
[2.0.0]  20200821
 Deprecated:
project
will only support functions likelambda x: x[0]
instead oflambda x, y, z: x
in the future.  Added: Support for complexvalued forms:
BilinearForm
andLinearForm
now take an optional argumentdtype
which defaults tonp.float64
but can be alsonp.complex64
.  Added:
Dofs.__or__
andDofs.__add__
, for merging degreeoffreedom sets (i.e.Dofs
objects) using
and+
operators.  Added:
Dofs.drop
andDofs.keep
, for further filtering the degreeoffreedom sets  Removed: Support for oldstyle decorators
bilinear_form
,linear_form
, andfunctional
(deprecated since 1.0.0).  Fixed:
FacetBasis
did not initialize withElementQuadP
.
[1.2.0]  20200707
 Added:
MeshQuad._splitquads
aliased asMeshQuad.to_meshtri
.  Added:
Mesh.__add__
, for merging meshes using+
operator: duplicated nodes are joined.  Added:
ElementHexS2
, a 20node quadratic hexahedral serendipity element.  Added:
ElementLineMini
, MINIelement for onedimensional mesh.  Fixed:
Mesh3D.boundary_edges
was broken in case of hexahedral meshes.  Fixed:
skfem.utils.project
did not work forElementGlobal
.
[1.1.0]  20200518
 Added:
ElementTetMini
, MINIelement for tetrahedral mesh.  Fixed:
Mesh3D.boundary_edges
incorrectly returned all edges where both nodes are on the boundary.
[1.0.0]  20200422
 Deprecated: Oldstyle form constructors
bilinear_form
,linear_form
, andfunctional
.  Changed:
Basis.interpolate
returnsDiscreteField
objects instead of ndarray tuples.  Changed:
Basis.interpolate
works now properly for vectorial and highorder elements by interpolating all components and higher order derivatives.  Changed:
Form.assemble
accepts now any keyword arguments (with typeDiscreteField
) that are passed over to the forms.  Changed: Renamed
skfem.importers
toskfem.io
.  Changed: Renamed
skfem.models.helpers
toskfem.helpers
.  Changed:
skfem.utils.solve
will now expand also the solutions of eigenvalue problems.  Added: Newstyle form constructors
BilinearForm
,LinearForm
, andFunctional
.  Added:
skfem.io.json
for serialization of meshes to/from jsonfiles.  Added:
ElementLinePp
, pth order onedimensional elements.  Added:
ElementQuadP
, pth order quadrilateral elements.  Added:
ElementQuadDG
for transforming quadrilateral H1 elements to DG elements.  Added:
ElementQuadBFS
, BognerFoxSchmit element for biharmonic problems.  Added:
ElementTriMini
, MINIelement for Stokes problems.  Added:
ElementComposite
for using multiple elements in one bilinear form.  Added:
ElementQuadS2
, quadratic Serendipity element.  Added:
ElementLineHermite
, cubic Hermite element for EulerBernoulli beams.  Added:
Mesh.define_boundary
for defining named boundaries.  Added:
Basis.find_dofs
for finding degreeoffreedom indices.  Added:
Mesh.from_basis
for defining highorder meshes.  Added:
Basis.split
for splitting multicomponent solutions.  Added:
MortarMapping
with basic support for mortar methods in 2D.  Added:
Basis
constructors now acceptquadrature
keyword argument for specifying a custom quadrature rule.
Project details
Release history Release notifications  RSS feed
Download files
Download the file for your platform. If you're not sure which to choose, learn more about installing packages.
Source Distribution
Built Distribution
Hashes for scikit_fem7.0.1py3noneany.whl
Algorithm  Hash digest  

SHA256  6a5c708f6bdf9e466c49a310970c2dd5c462b5b8131f53eb0361b5b53ea05215 

MD5  8233204c0a90b7eb0ae2428c9ffc748e 

BLAKE2256  1cd2b28f4256579062d27da5745649e44761dae259545cb3a56b2da3daca6d6d 