repository·master·Indexed 20 days ago
https://github.com/kinnala/scikit-femA pure Python 3.10+ library for finite element assembly. It specializes in transforming bilinear and linear forms into sparse matrices and vectors, supporting various element types across 1D, 2D, and 3D geometries. The library includes modules for mesh management (skfem.mesh), assembly (skfem.assembly), finite element definitions (skfem.element), and utilities for solving linear systems. It supports optional integrations with meshio, matplotlib, JAX for automatic differentiation, and PETSc for high-performance parallel assembly.
scikit-fem is a pure Python 3.8+ library designed for finite element assembly. Its primary function is to transform bilinear forms into sparse matrices and linear forms into vectors.
Supported mesh types include:
scikit-fem is a lightweight Python library designed for the assembly of finite element matrices ($A$) and vectors ($b$) used to solve partial differential equations (PDEs). It transforms a PDE into a system of linear equations $Ax=b$.
To solve a problem, a user typically follows these steps:
meshio, pacopy, and pyamg.scikit-fem focuses specifically on the assembly process. Users are responsible for implementing higher-level logic such as nonlinear iterations (e.g., Newton's method), adaptive mesh refinement, and specific boundary condition implementations (though helper routines are provided).The number of degrees of freedom (dofs) differs between basis types on the same mesh:
You can check the number of dofs using the .zeros().shape[0] attribute on a basis object.
# Example output for a mesh with 9 nodes and 8 elements
print(f'{basis_p1.zeros().shape[0]} dofs in P1 == {mesh.p.shape[1]} nodes in the mesh')
print(f'{basis_p0.zeros().shape[0]} dofs in P0 == {mesh.t.shape[1]} elements in the mesh')Finite element integration relies on sampling functions at specific locations called quadrature points.
basis.quadrature returns the points and weights for the reference element (e.g., the unit triangle).basis.mapping.F(points). This maps the local reference coordinates to global mesh coordinates.The resulting global_points array is organized as (coordinate, element_index, quadrature_index).
points, weights = basis_p1.quadrature
# Map local quadrature points to global mesh coordinates
global_points = basis_p1.mapping.F(points)
# global_points.shape is (2, num_elements, num_points_per_element)In skfem, a Mesh (such as MeshTri) is represented by two primary components:
mesh.p): A 2D array where the first row contains the x-coordinates and the second row contains the y-coordinates of the vertices.mesh.t): An array of connections (indices) that define how points are grouped to form elements (e.g., triangles). Each row/column in the topology specifies the indices of the points that form a single element.For a 2D triangular mesh, mesh.t contains the indices of the three vertices for each triangle.
import skfem
import numpy as np
mesh = skfem.MeshTri()
print(mesh.p) # Vertex coordinates
print(mesh.t) # Element topologyTo use an existing finite element function (like a previous solution in an iterative solver) within a BilinearForm or LinearForm, you must interpolate it from the nodes to the quadrature points using basis.interpolate(x).
Inside the form definition, the third argument w (for BilinearForm) or the single argument w (for LinearForm) is a dictionary containing user-provided arguments. You can access these via w.key or w['key']. By default, w.x provides global coordinates and w.h provides the local mesh parameter.
When calling .assemble(), pass the interpolated function as a keyword argument.
import skfem as fem
from skfem.helpers import grad, dot
# 1. Define the form using a keyword argument 'u_k' from the 'w' dictionary
@fem.BilinearForm
def bilinf(u, v, w):
return (w.u_k + .1) * dot(grad(u), grad(v))
# 2. Prepare the mesh and basis
m = fem.MeshTri().refined(3)
basis = fem.Basis(m, fem.ElementTriP1())
# 3. Assume 'x' is your current solution vector
x = 0. * basis.x
# 4. Interpolate 'x' and pass it to assemble
A = bilinf.assemble(basis, u_k=basis.interpolate(x))skfem.autodiff module allows you to linearize nonlinear variational forms to solve them using the Newton method. This requires the JAX optional dependency. You can use the NonlinearForm decorator to define these forms.Beyond standard elements, scikit-fem supports:
To ensure that different basis sets share a common set of quadrature points (which allows for exact integration of the highest-order set and provides a common space for computations), it is recommended to construct the highest-order basis first and then derive lower-order bases from it using the with_element() method.
For example, you can derive a P0 (constant) basis from a P1 (linear) basis.
import skfem
# Assume mesh and basis_p1 (P1) are already defined
basis_p0 = basis_p1.with_element(skfem.ElementTriP0())skfem.utils.enforce or skfem.utils.penalize. enforce works by modifying the matrix rows and diagonals (setting the row to zero and diagonal to one), while penalize is an alternative method for setting boundary conditions.When debugging or implementing custom forms, it is critical to know that the integrand function must return a 2D numpy array. The shape is always:
[number of elements] x [number of quadrature points per element]
You can verify this by using a debugger (like pdb) inside your form function and checking the shape of the returned value or the shape of component derivatives like u.grad[0].
After assembly, the unknowns in the linear system $Ax = b$ are ordered in the vector $x$ based on the mesh and element type. You can inspect how DOFs map to topological entities using the following attributes on a Basis object:
basis.nodal_dofs: DOFs corresponding to mesh nodes (vertices).basis.facet_dofs: DOFs corresponding to facets (shared between two elements; in 2D, these are edges).basis.edge_dofs: DOFs corresponding to edges (used in 3D meshes).basis.interior_dofs: DOFs that are internal to the element and not shared with neighbors.Note on nomenclature: In scikit-fem, facets are the entities shared between two elements in any dimension. For 2D meshes, facets are what are traditionally called edges.
from skfem import *
m = MeshHex()
basis = Basis(m, ElementHex2())
# Accessing different DOF types
nodal = basis.nodal_dofs
facets = basis.facet_dofs
edges = basis.edge_dofs
interior = basis.interior_dofs