scikit-fem: finite element assembly in pure Python

One job: bilinear forms into sparse matrices

The README is precise about the scope. scikit-fem is a pure Python 3.10+ library for performing finite element assembly, and its main purpose is the transformation of bilinear forms into sparse matrices and linear forms into vectors. That is the entire contract. It is not a mesh generator with a solver bolted on, it is not a multiphysics framework, and it does not try to model your problem for you. What it does is remove the part of finite element work that is tedious and error prone, namely the index bookkeeping that turns element-local integrals into global sparse matrices, and it does that with an interface where a form is a decorated Python function. The audience is a Python-literate scientist or student who knows the weak form of their PDE and wants to assemble it, not a commercial CAE user. The library is BSD-3-Clause licensed, the default branch is master, and the last push to the repository was on 2026-10-20.

No compiled code, and only two required dependencies

The library advertises minimal dependencies and contains no compiled code, which is the property that decides whether it fits your environment. The project metadata lists the runtime dependencies as exactly numpy and scipy, with a requires-python of >=3.10 and classifiers covering Python 3.10 through 3.14. The build system is flit_core, and the module is distributed under the name skfem, so the import is skfem even though the distribution is scikit-fem. That matters for environments with strict policies about compiled extensions, since there is nothing to build. It also means the performance story rests entirely on NumPy and SciPy doing the heavy lifting, which is the subject of the benchmark further down. The package is also usable in a browser through Pyodide, and the repository ships a paper/ directory with the Journal of Open Source Software article for citation, with the DOI 10.21105/joss.02369.

Installing with the extras that make examples runnable

The most recent release installs with one command.

pip install scikit-fem[all]

Removing the [all] extra skips the optional dependencies, which the README names as meshio for mesh input and output and matplotlib for creating simple visualizations, leaving only numpy and scipy. A few submodules need more: skfem.autodiff requires jax and jaxlib, and skfem.supermeshing requires shapely. The dependency picture is therefore layered rather than all or nothing, and the choice that matters most in practice is whether to take meshio, because without it you cannot load a mesh file from disk and the examples that read .msh files will not run. The test requirements listed in requirements.txt go further still, adding ipython, pytest, pyamg, mypy, flake8, sphinx, jax, jaxlib, pep517 and shapely, which tells you the development environment is considerably larger than the install a user needs.

The Poisson example, and how a form becomes a matrix

The README's first example solves the Poisson problem, and it is worth reading line by line because it shows the whole model. A mesh is created and refined, a basis is built from the mesh and an element, and then the two forms are defined.

from skfem import *
from skfem.helpers import dot, grad

# create the mesh
mesh = MeshTri().refined(4)

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)

The decorators are the mechanism. A bilinear form is a function of two trial and test arguments and a set of quadrature points, written as plain Python with array operations, and the assembler works out the element contributions and scatters them into a sparse matrix. The linear form works the same way and produces a vector. The signature takes three arguments, with the third being the quadrature points, which the form here ignores.

Boundary conditions and the solve follow in the same style.

# Dirichlet boundary conditions
A, b = enforce(A, b, D=mesh.boundary_nodes())

# solve the linear system
x = solve(A, b)

The enforce call returns the modified matrix and right-hand side, and solve hands the sparse system to SciPy. The result is plotted with mesh.plot or written out with mesh.save for a VTK file, both of which need the optional extras.

Meshes from constructors, from files, or from your own points

Meshes come from three directions, and the README shows all of them. The special constructors build simple shapes, so a one-dimensional mesh from a coordinate array, a triangle from explicit points and connectivity, or a tensor-product tetrahedral mesh.

import numpy as np
from skfem import MeshLine, MeshTri, MeshTet

mesh = MeshLine(np.array([0., .5, 1.]))
mesh = MeshTri.load("docs/examples/meshes/square.msh")  # requires meshio
mesh = MeshTet.init_tensor(*((np.linspace(0, 1, 60),) * 3))

The load call reads a mesh file and needs meshio installed, which is the concrete reason to take the [all] extra even for a quick trial. A mesh can also be constructed from your own points and elements, with the comment in the Poisson example showing the alternative form of MeshTri(points, elements). The element itself is chosen separately from the mesh, and the element catalogue is wide: one-dimensional, triangular, quadrilateral, tetrahedral and hexahedral elements are supported, plus special elements such as Raviart-Thomas, Nedelec, MINI, Crouzeix-Raviart and Argyris. A second-order tetrahedral basis is a one-line change.

from skfem import Basis, ElementTetP2

basis = Basis(mesh, ElementTetP2())  # quadratic tetrahedron
A = laplace.assemble(basis)  # type: scipy.sparse.csr_matrix

The benchmark says assembly is usually not the bottleneck

The README runs a benchmark from docs/examples/performance.py that compares the time spent on assembly with the time spent on the linear solve, solving the Laplace equation with linear tetrahedral elements and the default direct sparse solver of scipy.sparse.linalg.spsolve, with the numbers calculated on a Macbook Pro M2. The README warns that the timings are only illustrative, since they depend on the type of element used, the number of quadrature points, the type of linear solver and the complexity of the forms. The table runs from 64 degrees of freedom at 0.00633 seconds of assembly and 0.00159 of solve, through 8000 degrees of freedom at 0.05288 and 0.33537, to 125000 at 1.10499 with a nan solve time, and 1030301 degrees of freedom at 12.58563 seconds of assembly with nan for the solve. The conclusion the README draws is that in many cases the time spent on the linear solve will significantly dominate, and a Python package such as scikit-fem can be fast enough for assembling the finite element matrices. Read that last part carefully. It is a claim about assembly, not about solving.

Where a Python assembly step stops making sense

The benchmark table is also the clearest statement of the limits. The crossover happens somewhere between a few thousand and a few tens of thousands of degrees of freedom, where the solve time overtakes assembly and then grows far faster, reaching 62.69 seconds of solve against 0.54 of assembly at 64000 degrees of freedom. The nan entries in the last rows are the solve timings running out of usable measurement, not a claim that the solve became free. So the case for a pure Python assembler weakens as the problem gets large, and it weakens on the solve side rather than the assembly side. Two other cases point elsewhere. The first is an environment where you need a compiled finite element stack with a maintained solver and preconditioners, since scikit-fem assembles the matrix and then hands it to SciPy, so the solver is whatever SciPy offers rather than a purpose-built one. The second is a codebase that already uses a mature finite element package, where the cost of a second framework is not worth the readability of a Python decorator. The library is also a research code in the sense that the project is a single authored package with a journal paper, an issue tracker and a Discussions board, not a foundation-governed project with a long support history.

Editorial conclusion

scikit-fem fits a Python researcher or engineer who needs a working finite element model quickly, in a language their colleagues already use, and who can live with a pure Python assembly step. It does not fit a production solver that needs the assembly to stay negligible next to the solve, because the README's own benchmark shows the solve overtaking assembly and then running away, with the largest degrees of freedom entries recording nan for the solve time. Before you adopt it, read the benchmark table in the README with your own element type and quadrature count in mind, since the numbers depend on both, and install it with the optional extras, because meshio and matplotlib are what turn the examples into something you can actually look at.

Frequently asked questions

Official sources

1. kinnala/scikit-fem on GitHub
2. License: BSD-3-Clause
3. Project website
4. README
5. Releases