scikit-fem

scikit-fem 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 one-dimensional, triangular, quadrilateral, tetrahedral and hexahedral finite elements
includes special elements such as Raviart-Thomas, Nédélec, MINI, Crouzeix-Raviart, Argyris, ...

Installation

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.

Examples

Solve the Poisson problem (see also ex01.py):

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

# create the mesh
m = MeshTri().refined(4)
# or, with your own points and elements:
# m = MeshTri(points, elements)

e = ElementTriP1()
basis = Basis(m, e)  # shorthand for CellBasis

@BilinearForm
def laplace(u, v, _):
    return dot(grad(u), grad(v))

@LinearForm
def rhs(v, _):
    return 1.0 * v

A = laplace.assemble(basis)
b = rhs.assemble(basis)

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

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

# plot the solution
from skfem.visuals.matplotlib import plot, savefig
plot(m, x, shading='gouraud', colorbar=True)
savefig('solution.png')

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.

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.

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

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 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.

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/scikit-fem-docker-action. The releases are tested in kinnala/scikit-fem-release-tests.

Licensing

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.

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:

Reporting a bug
Writing an example
Improving the tests
Finding typos in the documentation.

By contributing code to scikit-fem, you are agreeing to release it under BSD-3-Clause, see LICENSE.md.

Citing the library

We appreciate if you cite the following article in your scientific works:

@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}
}

Changelog

The format is based on Keep a Changelog, and this project adheres to Semantic Versioning with respect to documented and/or tested features.

Unreleased

Changed: 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)
Changed: Writing w['u'] and w.u inside the form definition is now equivalent (previously w.u == w['u'].value)
Changed: Mesh.draw now uses matplotlib by default, calling Mesh.draw("vedo") uses vedo
Changed: skfem.visuals.matplotlib now uses jet as the default colormap
Changed: BoundaryFacetBasis is now an alias of FacetBasis instead of other way around
Deprecated: DiscreteField.value remains for backwards-compatibility but is now deprecated and can be dropped
Added: Mesh.plot, a wrapper to skfem.visuals.*.plot
Added: Basis.plot, a wrapper to skfem.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 and load now preserve the orientation of boundaries and interfaces
Added: OrientedBoundary which is a subclass of ndarray for facet index arrays with the orientation information (0 or 1 per facet) available as OrientedBoundary.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 if boundaries=True is given
Added: Mesh.facets_satisfying allows specifying the keyword argument normal for orienting the resulting interface
Added: FacetBasis constructor now has the keyword argument side which allows changing the side of the facet used to calculate the basis function values and gradients
Fixed: Improvements to backwards compatibility in asm/assemble keyword arguments
Fixed: Save format issue with meshio 5.3.0+
Fixed: CellBasis did not properly support elements argument
Fixed: Basis.interpolate did not properly interpolate all components of ElementComposite

[5.2.0] - 2021-12-27

Added: Basis.project, a more general and easy to use alternative for projection
Added: Basis and FacetBasis kwargs elements and facets can now be a string refering to subdomain and boundary tags
Added: ElementQuadRT0, lowest-order quadrilateral Raviart-Thomas element
Fixed: Functional returned only the first component for forms with non-scalar output

[5.1.0] - 2021-11-30

Added: skfem.helpers.mul for matrix multiplication
Added: Basis.split will now also split ElementVector into its components
Fixed: ElementDG was not included in the wildcard import
Fixed: Automatic visualization of MeshTri2 and MeshQuad2 in Jupyter notebooks raised exception

[5.0.0] - 2021-11-21

Changed: meshio is now an optional dependency
Changed: ElementComposite uses DiscreteField() 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 of MeshLine1, MeshTri1 and MeshQuad1
Fixed: Refinement of quadratic meshes will now fall back to the refinement algorithm of first-order 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 of Basis.get_dofs
Deprecated: Merging DofsView objects via + and | because of surprising behavior in some edge cases

[4.0.1] - 2021-10-15

Fixed: MappingIsoparametric can now be pickled

[4.0.0] - 2021-09-27

Added: Mesh.save/Mesh.load now exports/imports Mesh.subdomains and Mesh.boundaries
Added: 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']
Added: 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 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 flag boundaries_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 for ElementTriDG(ElementTriP1()) etc.
Added: ElementTriSkeletonP0 and ElementTriSkeletonP1 for defining Lagrange multipliers on the skeleton mesh (see Example 40)
Added: TrilinearForm for assembling a sparse 3-tensor, 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 and ElementWedge1, the lowest order wedge mesh and element
Added: ElementTriP3, cubic triangular Lagrange element
Added: ElementTriP4, quartic triangular Lagrange element
Added: ElementTri15ParamPlate, 15-parameter nonconforming triangular element for plates
Added: ElementTriBDM1, the lowest order Brezzi-Douglas-Marini 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 unparseable cell_sets inserted by meshio in MSH 4.1
Changed: Mesh string representation is now more informative
Changed: Form.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 convenience
Changed: MeshLine is now a function which initializes MeshLine1 and not an alias to MeshLine1
Changed: FacetBasis is now a shorthand for BoundaryFacetBasis and no longer initializes InteriorFacetBasis or MortarFacetBasis if the keyword argument side is passed to the constructor
Removed: the deprecated Mesh.define_boundary method

[3.2.0] - 2021-08-02

Added: ElementTriCCR and ElementTetCCR, conforming Crouzeix-Raviart 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 (and CellBasis.probes, CellBasis.interpolator) was not working properly for a small number of elements (<5) or a large number of input points (>1000)
Fixed: MeshTet and MeshTri.element_finder are now more robust against degenerate elements
Fixed: Mesh.element_finder (and CellBasis.probes, CellBasis.interpolator) raises exception if the query point is outside of the domain

[3.1.0] - 2021-06-18

Added: Basis, a shorthand for CellBasis
Added: CellBasis, a new preferred name for InteriorBasis
Added: BoundaryFacetBasis, a new preferred name for ExteriorFacetBasis
Added: utils.penalize, an alternative to condense and enforce for essential boundary conditions
Added: InteriorBasis.point_source, with ex38
Added: ElementTetDG, similar to ElementTriDG for tetrahedral meshes
Fixed: MeshLine1.element_finder

[3.0.0] - 2021-04-19

Added: Completely rewritten 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 future
Added: New quadratic mesh types: MeshTri2, MeshQuad2, MeshTet2 and MeshHex2
Added: InteriorBasis.probes; like InteriorBasis.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 via InteriorBasis.interpolator
Added: Mesh.with_boundaries, a functional replacement to Mesh.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 of skfem.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 of skfem.utils.projection
Deprecated: Mesh.define_boundary in favor of Mesh.with_boundaries
Removed: Mesh.{refine,scale,translate}; the replacements are Mesh.{refined,scaled,translated}
Removed: skfem.models.helpers; available as skfem.helpers
Removed: DiscreteField.{f,df,ddf,hod}; available as DiscreteField.{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 of Mesh.boundaries after refine
Changed: skfem.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()
Fixed: High memory usage in skfem.visuals.matplotlib related to 1D plotting

[2.5.0] - 2021-02-13

Deprecated: side keyword argument to FacetBasis in favor of the more explicit InteriorFacetBasis and MortarFacetBasis.
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 applying functools.partial to the form function wrapped by Form.
Fixed: Include explicit Python 3.9 support.

[2.4.0] - 2021-01-20

Deprecated: List and tuple keyword argument types to asm.
Deprecated: Mesh2D.mirror in favor of the more general Mesh.mirrored.
Deprecated: Mesh.refine, Mesh.scale and Mesh.translate in favor of Mesh.refined, Mesh.scaled and Mesh.translated.
Added: Mesh.refined, Mesh.scaled, and Mesh.translated. The new methods return a copy instead of modifying self.
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 one-dimensional meshes.

[2.3.0] - 2020-11-24

Added: ElementLineP0, one-dimensional piecewise constant element.
Added: skfem.helpers.curl now calculates the rotated gradient for two-dimensional elements.
Added: MeshTet.init_ball for meshing a ball.
Fixed: ElementQuad0 was not compatible with FacetBasis.

[2.2.3] - 2020-10-16

Fixed: Remove an unnecessary dependency.

[2.2.2] - 2020-10-15

Fixed: Make the preconditioner in TestEx32 more robust.

[2.2.1] - 2020-10-15

Fixed: Remove tests from the PyPI distribution.

[2.2.0] - 2020-10-14

Deprecated: L2_projection will be replaced by project.
Deprecated: derivative will be replaced by project.
Added: MeshTet.element_finder and MeshLine.element_finder for using InteriorBasis.interpolator.
Added: ElementTriCR, the nonconforming Crouzeix-Raviart element for Stokes flow.
Added: ElementTetCR, tetrahedral nonconforming Crouzeix-Raviart element.
Added: ElementTriHermite, an extension of ElementLineHermite to triangular meshes.
Fixed: Fix Mesh.validate for unsigned Mesh.t.

[2.1.1] - 2020-10-01

Fixed: Further optimizations to Mesh3D.boundary_edges: tested to run on a laptop with over 10 million elements.

[2.1.0] - 2020-09-30

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] - 2020-08-21

Deprecated: project will only support functions like lambda x: x[0] instead of lambda x, y, z: x in the future.
Added: Support for complex-valued forms: BilinearForm and LinearForm now take an optional argument dtype which defaults to np.float64 but can be also np.complex64.
Added: Dofs.__or__ and Dofs.__add__, for merging degree-of-freedom sets (i.e. Dofs objects) using | and + operators.
Added: Dofs.drop and Dofs.keep, for further filtering the degree-of-freedom sets
Removed: Support for old-style decorators bilinear_form, linear_form, and functional (deprecated since 1.0.0).
Fixed: FacetBasis did not initialize with ElementQuadP.

[1.2.0] - 2020-07-07

Added: MeshQuad._splitquads aliased as MeshQuad.to_meshtri.
Added: Mesh.__add__, for merging meshes using + operator: duplicated nodes are joined.
Added: ElementHexS2, a 20-node quadratic hexahedral serendipity element.
Added: ElementLineMini, MINI-element for one-dimensional mesh.
Fixed: Mesh3D.boundary_edges was broken in case of hexahedral meshes.
Fixed: skfem.utils.project did not work for ElementGlobal.

[1.1.0] - 2020-05-18

Added: ElementTetMini, MINI-element for tetrahedral mesh.
Fixed: Mesh3D.boundary_edges incorrectly returned all edges where both nodes are on the boundary.

[1.0.0] - 2020-04-22

Deprecated: Old-style form constructors bilinear_form, linear_form, and functional.
Changed: Basis.interpolate returns DiscreteField objects instead of ndarray tuples.
Changed: Basis.interpolate works now properly for vectorial and high-order elements by interpolating all components and higher order derivatives.
Changed: Form.assemble accepts now any keyword arguments (with type DiscreteField) that are passed over to the forms.
Changed: Renamed skfem.importers to skfem.io.
Changed: Renamed skfem.models.helpers to skfem.helpers.
Changed: skfem.utils.solve will now expand also the solutions of eigenvalue problems.
Added: New-style form constructors BilinearForm, LinearForm, and Functional.
Added: skfem.io.json for serialization of meshes to/from json-files.
Added: ElementLinePp, p-th order one-dimensional elements.
Added: ElementQuadP, p-th order quadrilateral elements.
Added: ElementQuadDG for transforming quadrilateral H1 elements to DG elements.
Added: ElementQuadBFS, Bogner-Fox-Schmit element for biharmonic problems.
Added: ElementTriMini, MINI-element for Stokes problems.
Added: ElementComposite for using multiple elements in one bilinear form.
Added: ElementQuadS2, quadratic Serendipity element.
Added: ElementLineHermite, cubic Hermite element for Euler-Bernoulli beams.
Added: Mesh.define_boundary for defining named boundaries.
Added: Basis.find_dofs for finding degree-of-freedom indices.
Added: Mesh.from_basis for defining high-order meshes.
Added: Basis.split for splitting multicomponent solutions.
Added: MortarMapping with basic support for mortar methods in 2D.
Added: Basis constructors now accept quadrature keyword argument for specifying a custom quadrature rule.

Name		Name	Last commit message	Last commit date
Latest commit History 2,741 Commits
.github/workflows		.github/workflows
docs		docs
paper		paper
skfem		skfem
tests		tests
.flake8		.flake8
.gitignore		.gitignore
.zenodo.json		.zenodo.json
LICENSE		LICENSE
Makefile		Makefile
README.md		README.md
mypy.ini		mypy.ini
pyproject.toml		pyproject.toml
requirements.txt		requirements.txt
setup.cfg		setup.cfg
shell.nix		shell.nix

License

HadrienNU/scikit-fem

Folders and files

Latest commit

History

Repository files navigation