Detailed API description

This section contains API documentation for the most commonly used interfaces of the library.

Module: skfem.mesh

This module defines finite element meshes. Meshes can be created using built-in constructors or loaded from external formats using meshio. The supported types are

Default constructor creates a mesh for the unit square:

>>> from skfem.mesh import MeshTri
>>> MeshTri()
<skfem MeshTri1 object>
  Number of elements: 2
  Number of vertices: 4
  Number of nodes: 4

Each mesh type has several constructors; see the docstring, e.g., help(MeshTri) or click MeshTri in the online documentation. Importing from external formats can be done with the constructor load().

Abstract class: Mesh

class skfem.mesh.Mesh(doflocs: ~numpy.ndarray, t: ~numpy.ndarray, _boundaries: ~typing.Dict[str, ~numpy.ndarray] | None = None, _subdomains: ~typing.Dict[str, ~numpy.ndarray] | None = None, elem: ~typing.Type[~skfem.element.element.Element] = <class 'skfem.element.element.Element'>, affine: bool = False, sort_t: bool = False, validate: bool = True)

A mesh defined by Element class.

Mesh is defined as a combination of elements/cells by specifying the spatial locations of the finite element nodes.

elements_satisfying(test: Callable[[ndarray], ndarray]) → ndarray

Return elements whose midpoints satisfy some condition.

Parameters:

test – A function which returns True for the element midpoints that are to be included in the return set.

facets_satisfying(test: Callable[[ndarray], ndarray], boundaries_only: bool = False, normal: ndarray | None = None) → ndarray

Return facets whose midpoints satisfy some condition.

Parameters:

  • test – A function which returns True for the facet midpoints that are to be included in the return set.

  • boundaries_only – If True, include only boundary facets.

  • normal – If given, used to orient the set of facets.

classmethod load(filename: str, out: List[str] | None = None, **kwargs)

Load a mesh using meshio.

Parameters:

  • filename – The filename of the mesh file.

  • out – Optional list of meshio.Mesh attribute names, overwrite with the corresponding data. E.g., ['point_data', 'cell_data'].

nodes_satisfying(test: Callable[[ndarray], ndarray], boundaries_only: bool = False) → ndarray

Return nodes that satisfy some condition.

Parameters:

  • test – A function which returns True for the set of nodes that are to be included in the return set.

  • boundaries_only – If True, include only boundary facets.

refined(times_or_ix: int | ndarray = 1)

Return a refined mesh.

Parameters:

times_or_ix – Either an integer giving the number of uniform refinements or an array of element indices for adaptive refinement.

save(filename: str | PathLike, point_data: Dict[str, ndarray] | None = None, cell_data: Dict[str, ndarray] | None = None, **kwargs) → None

Export the mesh and fields using meshio.

Parameters:

  • filename – The output filename, with suffix determining format; e.g. .msh, .vtk, .xdmf

  • point_data – Data related to the vertices of the mesh.

  • cell_data – Data related to the elements of the mesh.

with_boundaries(boundaries: Dict[str, Callable[[ndarray], ndarray] | ndarray], boundaries_only: bool = True)

Return a copy of the mesh with named boundaries.

Parameters:

  • boundaries – A dictionary of index arrays or lambda functions with the names of the boundaries as keys. The midpoint of the facet should return True for the corresponding lambda function if the facet belongs to the boundary.

  • boundaries_only – If True, consider only facets on the boundary of the domain.

with_subdomains(subdomains: Dict[str, Callable[[ndarray], ndarray] | ndarray])

Return a copy of the mesh with named subdomains.

Parameters:

subdomains – A dictionary of lambda functions with the names of the subdomains as keys. The midpoint of the element should return True for the corresponding lambda function if the element belongs to the subdomain.

Class: MeshTri

skfem.mesh.MeshTri

alias of MeshTri1

class skfem.mesh.MeshTri1(doflocs: ~numpy.ndarray = <factory>, t: ~numpy.ndarray = <factory>, _boundaries: ~typing.Dict[str, ~numpy.ndarray] | None = None, _subdomains: ~typing.Dict[str, ~numpy.ndarray] | None = None, elem: ~typing.Type[~skfem.element.element.Element] = <class 'skfem.element.element_tri.element_tri_p1.ElementTriP1'>, affine: bool = True, sort_t: bool = True, validate: bool = True)

A standard first-order triangular mesh.

__init__(doflocs: ~numpy.ndarray = <factory>, t: ~numpy.ndarray = <factory>, _boundaries: ~typing.Dict[str, ~numpy.ndarray] | None = None, _subdomains: ~typing.Dict[str, ~numpy.ndarray] | None = None, elem: ~typing.Type[~skfem.element.element.Element] = <class 'skfem.element.element_tri.element_tri_p1.ElementTriP1'>, affine: bool = True, sort_t: bool = True, validate: bool = True) → None
classmethod init_circle(nrefs: int = 3, smoothed: bool = False) → Mesh2D

Initialize a circle mesh.

Works by repeatedly refining the following mesh and moving new nodes to the boundary:

       *
     / | \
   /   |   \
 /     |     \
*------O------*
 \     |     /
   \   |   /
     \ | /
       *
Parameters:

  • nrefs – Number of refinements, by default 3.

  • smoothed – If True, apply smoothing after each refine.

classmethod init_lshaped() → Mesh2D

Initialize a mesh for the L-shaped domain.

The mesh topology is as follows:

*-------*
| \     |
|   \   |
|     \ |
|-------O-------*
|     / | \     |
|   /   |   \   |
| /     |     \ |
*---------------*
classmethod init_refdom()

Initialize a mesh corresponding to the reference domain.

classmethod init_sqsymmetric() → Mesh2D

Initialize a symmetric mesh of the unit square.

The mesh topology is as follows:

*------*------*
|\     |     /|
|  \   |   /  |
|    \ | /    |
*------*------*
|    / | \    |
|  /   |   \  |
|/     |     \|
O------*------*
classmethod init_symmetric() → Mesh2D

Initialize a symmetric mesh of the unit square.

The mesh topology is as follows:

*------------*
|\          /|
|  \      /  |
|    \  /    |
|     *      |
|    /  \    |
|  /      \  |
|/          \|
O------------*
classmethod init_tensor(x: ndarray, y: ndarray)

Initialize a tensor product mesh.

The mesh topology is as follows:

*---------------*
|'-.|'-.|`'---._|
|---+---+-------|
|\  |\  |'.     |
| \ | \ |  '-.  |
|  \|  \|     '.|
*---------------*
Parameters:

  • x – The nodal coordinates in dimension x.

  • y – The nodal coordinates in dimension y.

classmethod load(filename: str, out: List[str] | None = None, **kwargs)

Load a mesh using meshio.

Parameters:

  • filename – The filename of the mesh file.

  • out – Optional list of meshio.Mesh attribute names, overwrite with the corresponding data. E.g., ['point_data', 'cell_data'].

class skfem.mesh.MeshTri2(doflocs: ~numpy.ndarray = <factory>, t: ~numpy.ndarray = <factory>, _boundaries: ~typing.Dict[str, ~numpy.ndarray] | None = None, _subdomains: ~typing.Dict[str, ~numpy.ndarray] | None = None, elem: ~typing.Type[~skfem.element.element.Element] = <class 'skfem.element.element_tri.element_tri_p2.ElementTriP2'>, affine: bool = False, sort_t: bool = False, validate: bool = True)

A quadratic triangular mesh.

classmethod init_circle(nrefs: int = 3, smoothed: bool = False) → MeshTri2

Initialize a circle mesh.

Works by repeatedly refining the following mesh and moving new nodes to the boundary:

       *
     / | \
   /   |   \
 /     |     \
*------O------*
 \     |     /
   \   |   /
     \ | /
       *
Parameters:

  • nrefs – Number of refinements, by default 3.

  • smoothed – If True, apply smoothing after each refine.

classmethod load(filename: str, out: List[str] | None = None, **kwargs)

Load a mesh using meshio.

Parameters:

  • filename – The filename of the mesh file.

  • out – Optional list of meshio.Mesh attribute names, overwrite with the corresponding data. E.g., ['point_data', 'cell_data'].

Class: MeshQuad

skfem.mesh.MeshQuad

alias of MeshQuad1

class skfem.mesh.MeshQuad1(doflocs: ~numpy.ndarray = <factory>, t: ~numpy.ndarray = <factory>, _boundaries: ~typing.Dict[str, ~numpy.ndarray] | None = None, _subdomains: ~typing.Dict[str, ~numpy.ndarray] | None = None, elem: ~typing.Type[~skfem.element.element.Element] = <class 'skfem.element.element_quad.element_quad1.ElementQuad1'>, affine: bool = False, sort_t: bool = False, validate: bool = True)

A standard first-order quadrilateral mesh.

If t is provided, order of vertices in each element should match the numbering:

3---2
|   |
0---1
__init__(doflocs: ~numpy.ndarray = <factory>, t: ~numpy.ndarray = <factory>, _boundaries: ~typing.Dict[str, ~numpy.ndarray] | None = None, _subdomains: ~typing.Dict[str, ~numpy.ndarray] | None = None, elem: ~typing.Type[~skfem.element.element.Element] = <class 'skfem.element.element_quad.element_quad1.ElementQuad1'>, affine: bool = False, sort_t: bool = False, validate: bool = True) → None
classmethod init_refdom()

Initialize a mesh corresponding to the reference domain.

classmethod init_tensor(x: ndarray, y: ndarray)

Initialize a tensor product mesh.

The mesh topology is as follows:

*-------------*
|   |  |      |
|---+--+------|
|   |  |      |
|   |  |      |
|   |  |      |
*-------------*
Parameters:

  • x – The nodal coordinates in dimension x.

  • y – The nodal coordinates in dimension y.

classmethod load(filename: str, out: List[str] | None = None, **kwargs)

Load a mesh using meshio.

Parameters:

  • filename – The filename of the mesh file.

  • out – Optional list of meshio.Mesh attribute names, overwrite with the corresponding data. E.g., ['point_data', 'cell_data'].

to_meshtri(x: ndarray | None = None, style: str | None = None)

Split each quadrilateral into two triangles.

Parameters:

  • x – Elementwise constant function to preserve. If given, returned as the second additional output parameter.

  • style – Optionally, specify a splitting style ‘x’ for crisscross splitting.

class skfem.mesh.MeshQuad2(doflocs: ~numpy.ndarray = <factory>, t: ~numpy.ndarray = <factory>, _boundaries: ~typing.Dict[str, ~numpy.ndarray] | None = None, _subdomains: ~typing.Dict[str, ~numpy.ndarray] | None = None, elem: ~typing.Type[~skfem.element.element.Element] = <class 'skfem.element.element_quad.element_quad2.ElementQuad2'>, affine: bool = False, sort_t: bool = False, validate: bool = True)

A quadratic quadrilateral mesh.

classmethod load(filename: str, out: List[str] | None = None, **kwargs)

Load a mesh using meshio.

Parameters:

  • filename – The filename of the mesh file.

  • out – Optional list of meshio.Mesh attribute names, overwrite with the corresponding data. E.g., ['point_data', 'cell_data'].

Class: MeshTet

skfem.mesh.MeshTet

alias of MeshTet1

class skfem.mesh.MeshTet1(doflocs: ~numpy.ndarray = <factory>, t: ~numpy.ndarray = <factory>, _boundaries: ~typing.Dict[str, ~numpy.ndarray] | None = None, _subdomains: ~typing.Dict[str, ~numpy.ndarray] | None = None, elem: ~typing.Type[~skfem.element.element.Element] = <class 'skfem.element.element_tet.element_tet_p1.ElementTetP1'>, affine: bool = True, sort_t: bool = False, validate: bool = True)

A standard first-order tetrahedral mesh.

__init__(doflocs: ~numpy.ndarray = <factory>, t: ~numpy.ndarray = <factory>, _boundaries: ~typing.Dict[str, ~numpy.ndarray] | None = None, _subdomains: ~typing.Dict[str, ~numpy.ndarray] | None = None, elem: ~typing.Type[~skfem.element.element.Element] = <class 'skfem.element.element_tet.element_tet_p1.ElementTetP1'>, affine: bool = True, sort_t: bool = False, validate: bool = True) → None
classmethod init_ball(nrefs: int = 3)

Initialize a ball mesh.

Parameters:

nrefs – Number of refinements, by default 3.

classmethod init_refdom()

Initialize a mesh corresponding to the reference domain.

classmethod init_tensor(x: ndarray, y: ndarray, z: ndarray)

Initialize a tensor product mesh.

Parameters:

  • x – The nodal coordinates in dimension x.

  • y – The nodal coordinates in dimension y.

  • z – The nodal coordinates in dimension z.

classmethod load(filename: str, out: List[str] | None = None, **kwargs)

Load a mesh using meshio.

Parameters:

  • filename – The filename of the mesh file.

  • out – Optional list of meshio.Mesh attribute names, overwrite with the corresponding data. E.g., ['point_data', 'cell_data'].

class skfem.mesh.MeshTet2(doflocs: ~numpy.ndarray = <factory>, t: ~numpy.ndarray = <factory>, _boundaries: ~typing.Dict[str, ~numpy.ndarray] | None = None, _subdomains: ~typing.Dict[str, ~numpy.ndarray] | None = None, elem: ~typing.Type[~skfem.element.element.Element] = <class 'skfem.element.element_tet.element_tet_p2.ElementTetP2'>, affine: bool = False, sort_t: bool = False, validate: bool = True)

A quadratic tetrahedral mesh.

classmethod init_ball(nrefs: int = 3) → MeshTet2

Initialize a ball mesh.

Parameters:

nrefs – Number of refinements, by default 3.

classmethod load(filename: str, out: List[str] | None = None, **kwargs)

Load a mesh using meshio.

Parameters:

  • filename – The filename of the mesh file.

  • out – Optional list of meshio.Mesh attribute names, overwrite with the corresponding data. E.g., ['point_data', 'cell_data'].

Class: MeshHex

skfem.mesh.MeshHex

alias of MeshHex1

class skfem.mesh.MeshHex1(doflocs: ~numpy.ndarray = <factory>, t: ~numpy.ndarray = <factory>, _boundaries: ~typing.Dict[str, ~numpy.ndarray] | None = None, _subdomains: ~typing.Dict[str, ~numpy.ndarray] | None = None, elem: ~typing.Type[~skfem.element.element.Element] = <class 'skfem.element.element_hex.element_hex1.ElementHex1'>, affine: bool = False, sort_t: bool = False, validate: bool = True)

A standard first-order hexahedral mesh.

If t is provided, order of vertices in each element should match the numbering:

  2---6
 /   /|
4---7 3
|   |/
1---5
__init__(doflocs: ~numpy.ndarray = <factory>, t: ~numpy.ndarray = <factory>, _boundaries: ~typing.Dict[str, ~numpy.ndarray] | None = None, _subdomains: ~typing.Dict[str, ~numpy.ndarray] | None = None, elem: ~typing.Type[~skfem.element.element.Element] = <class 'skfem.element.element_hex.element_hex1.ElementHex1'>, affine: bool = False, sort_t: bool = False, validate: bool = True) → None
classmethod init_tensor(x: ndarray, y: ndarray, z: ndarray)

Initialize a tensor product mesh.

Parameters:

  • x – The nodal coordinates in dimension x.

  • y – The nodal coordinates in dimension y.

  • z – The nodal coordinates in dimension z.

classmethod load(filename: str, out: List[str] | None = None, **kwargs)

Load a mesh using meshio.

Parameters:

  • filename – The filename of the mesh file.

  • out – Optional list of meshio.Mesh attribute names, overwrite with the corresponding data. E.g., ['point_data', 'cell_data'].

to_meshtet()

Split each hexahedron into six tetrahedra.

Class: MeshLine

skfem.mesh.MeshLine

alias of MeshLine1

class skfem.mesh.MeshLine1(doflocs: ~numpy.ndarray = <factory>, t: ~numpy.ndarray = <factory>, _boundaries: ~typing.Dict[str, ~numpy.ndarray] | None = None, _subdomains: ~typing.Dict[str, ~numpy.ndarray] | None = None, elem: ~typing.Type[~skfem.element.element.Element] = <class 'skfem.element.element_line.element_line_p1.ElementLineP1'>, affine: bool = True, sort_t: bool = False, validate: bool = True)

A one-dimensional mesh.

__init__(doflocs: ~numpy.ndarray = <factory>, t: ~numpy.ndarray = <factory>, _boundaries: ~typing.Dict[str, ~numpy.ndarray] | None = None, _subdomains: ~typing.Dict[str, ~numpy.ndarray] | None = None, elem: ~typing.Type[~skfem.element.element.Element] = <class 'skfem.element.element_line.element_line_p1.ElementLineP1'>, affine: bool = True, sort_t: bool = False, validate: bool = True) → None

Module: skfem.assembly

This module performs the finite element assembly. The basic workflow is the following:

  1. Initialize Mesh and Element.

>>> import skfem as fem
>>> m = fem.MeshTri()
>>> e = fem.ElementTriP1()
>>> m
<skfem MeshTri1 object>
  Number of elements: 2
  Number of vertices: 4
  Number of nodes: 4

  1. Create CellBasis or FacetBasis objects.

>>> basis = fem.CellBasis(m, e)

  1. Define the forms using BilinearForm, LinearForm, or Functional.

>>> form_a = fem.BilinearForm(lambda u, v, w: u * v)
>>> form_l = fem.LinearForm(lambda v, w: w.x[0] ** 2 * v)

Mathematically the above forms are

\[a(u,v) = \int_\Omega u v \,\mathrm{d}x \quad \mathrm{and} \quad l(v) = \int_\Omega x^2v \,\mathrm{d}x.\]

  1. Create the matrices/vectors using assemble().

>>> A = form_a.assemble(basis)
>>> b = form_l.assemble(basis)
>>> A.toarray()
array([[0.08333333, 0.04166667, 0.04166667, 0.        ],
       [0.04166667, 0.16666667, 0.08333333, 0.04166667],
       [0.04166667, 0.08333333, 0.16666667, 0.04166667],
       [0.        , 0.04166667, 0.04166667, 0.08333333]])
>>> b
array([0.0162037 , 0.15046296, 0.06712963, 0.09953704])
skfem.assembly.asm(form: ~skfem.assembly.form.form.Form, *args, to=<function _sum>, **kwargs) → Any

Perform finite element assembly.

A shorthand for skfem.assembly.Form.assemble() which, in addition, supports assembling multiple bases at once and summing the result.

Abstract class: AbstractBasis

Subclasses of AbstractBasis represent a global finite element basis evaluated at quadrature points.

class skfem.assembly.basis.AbstractBasis(mesh: ~skfem.mesh.mesh.Mesh, elem: ~skfem.element.element.Element, mapping: ~skfem.mapping.mapping.Mapping | None = None, intorder: int | None = None, quadrature: ~typing.Tuple[~numpy.ndarray, ~numpy.ndarray] | None = None, refdom: ~typing.Type[~skfem.refdom.Refdom] = <class 'skfem.refdom.Refdom'>, dofs: ~skfem.assembly.dofs.Dofs | None = None, disable_doflocs: bool = False)

Finite element basis at global quadrature points.

Please see the following implementations:

get_dofs(facets: Any | None = None, elements: Any | None = None, nodes: Any | None = None, skip: List[str] | None = None) → Any

Find global DOF numbers.

Accepts an array of facet/element indices. However, various argument types can be turned into an array of facet/element indices.

Get all boundary DOFs:

>>> import numpy as np
>>> from skfem import MeshTri, Basis, ElementTriP1
>>> m = MeshTri().refined()
>>> basis = Basis(m, ElementTriP1())
>>> basis.get_dofs().flatten()
array([0, 1, 2, 3, 4, 5, 7, 8], dtype=int32)

Get DOFs via a function query:

>>> import numpy as np
>>> from skfem import MeshTri, Basis, ElementTriP1
>>> m = MeshTri().refined()
>>> basis = Basis(m, ElementTriP1())
>>> basis.get_dofs(lambda x: np.isclose(x[0], 0)).flatten()
array([0, 2, 5], dtype=int32)

Get DOFs via named boundaries:

>>> import numpy as np
>>> from skfem import MeshTri, Basis, ElementTriP1
>>> m = (MeshTri()
...      .refined()
...      .with_boundaries({'left': lambda x: np.isclose(x[0], 0)}))
>>> basis = Basis(m, ElementTriP1())
>>> basis.get_dofs('left').flatten()
array([0, 2, 5], dtype=int32)

Get DOFs via named subdomains:

>>> from skfem import MeshTri, Basis, ElementTriP1
>>> m = (MeshTri()
...      .refined()
...      .with_subdomains({'left': lambda x: x[0] < .5}))
>>> basis = Basis(m, ElementTriP1())
>>> basis.get_dofs(elements='left').flatten()
array([0, 2, 4, 5, 6, 8], dtype=int32)

Further reduce the set of DOFs:

>>> import numpy as np
>>> from skfem import MeshTri, Basis, ElementTriArgyris
>>> m = MeshTri().refined()
>>> basis = Basis(m, ElementTriArgyris())
>>> basis.get_dofs(lambda x: np.isclose(x[0], 0)).nodal.keys()
dict_keys(['u', 'u_x', 'u_y', 'u_xx', 'u_xy', 'u_yy'])
>>> basis.get_dofs(lambda x: np.isclose(x[0], 0)).all(['u', 'u_x'])
array([ 0,  1, 12, 13, 30, 31], dtype=int32)

Skip some DOF names altogether:

>>> import numpy as np
>>> from skfem import MeshTri, Basis, ElementTriArgyris
>>> m = MeshTri().refined()
>>> basis = Basis(m, ElementTriArgyris())
>>> basis.get_dofs(lambda x: np.isclose(x[0], 0),
...                skip=['u_x', 'u_y']).nodal.keys()
dict_keys(['u', 'u_xx', 'u_xy', 'u_yy'])

Combine several boundaries into one call:

>>> import numpy as np
>>> from skfem import MeshTri, Basis, ElementTriP1
>>> m = (MeshTri()
...      .with_boundaries({'left': lambda x: np.isclose(x[0], 0),
...                        'right': lambda x: np.isclose(x[0], 1)}))
>>> basis = Basis(m, ElementTriP1())
>>> basis.get_dofs({'left', 'right'}).flatten()
array([0, 1, 2, 3], dtype=int32)
Parameters:

  • facets – An array of facet indices. If None, find facets by self.mesh.boundary_facets(). If callable, call self.mesh.facets_satisfying(facets) to get the facets. If string, use self.mesh.boundaries[facets] to get the facets. If list, tuple or set, use the combined facet indices.

  • elements – An array of element indices. See above.

  • nodes – An array of node indices.

  • skip – List of dofnames to skip.

interpolate(w: ndarray) → DiscreteField | Tuple[DiscreteField, ...]

Interpolate a solution vector to quadrature points.

Useful when a solution vector is needed in the forms, e.g., when evaluating functionals or when solving nonlinear problems.

Parameters:

w – A solution vector.

project(interp, **kwargs)

Perform \(L^2\) projection onto the basis.

See Performing projections for more information.

Class: CellBasis

skfem.assembly.Basis

alias of CellBasis

class skfem.assembly.CellBasis(mesh: Mesh, elem: Element, mapping: Mapping | None = None, intorder: int | None = None, elements: Any | None = None, quadrature: Tuple[ndarray, ndarray] | None = None, dofs: Dofs | None = None, disable_doflocs: bool = False)

For fields defined inside the domain.

CellBasis object is a combination of Mesh and Element.

>>> from skfem import *
>>> m = MeshTri.init_symmetric()
>>> e = ElementTriP1()
>>> basis = CellBasis(m, e)

The resulting objects are used in the assembly.

>>> from skfem.models.poisson import laplace
>>> K = asm(laplace, basis)
>>> K.shape
(5, 5)
__init__(mesh: Mesh, elem: Element, mapping: Mapping | None = None, intorder: int | None = None, elements: Any | None = None, quadrature: Tuple[ndarray, ndarray] | None = None, dofs: Dofs | None = None, disable_doflocs: bool = False)

Combine Mesh and Element into a set of precomputed global basis functions.

Parameters:

  • mesh – An object of type Mesh.

  • elem – An object of type Element.

  • mapping – An object of type skfem.mapping.Mapping. If None, uses mesh.mapping.

  • intorder – Optional integration order, i.e. the degree of polynomials that are integrated exactly by the used quadrature. Not used if quadrature is specified.

  • elements – Optional subset of element indices.

  • quadrature – Optional tuple of quadrature points and weights.

  • dofs – Optional Dofs object.

  • disable_doflocs – If True, the computation of global DOF locations is disabled. This may save memory on large meshes if DOF locations are not required.

interpolate(w: ndarray) → DiscreteField | Tuple[DiscreteField, ...]

Interpolate a solution vector to quadrature points.

Useful when a solution vector is needed in the forms, e.g., when evaluating functionals or when solving nonlinear problems.

Parameters:

w – A solution vector.

project(interp, elements=None, dtype=None)

Perform \(L^2\) projection onto the basis.

See Performing projections for more information.

Parameters:

  • interp – An object of type DiscreteField which is a function (to be projected) evaluated at global quadrature points. If a function is given, then DiscreteField is created by passing an array of global quadrature point locations to the function.

  • elements – Optionally perform the projection on a subset of elements. The values of the remaining DOFs are zero.

  • dtype – Set to np.complex64 or similar to use complex numbers.

Class: FacetBasis

skfem.assembly.BoundaryFacetBasis

alias of FacetBasis

class skfem.assembly.FacetBasis(mesh: Mesh, elem: Element, mapping: Mapping | None = None, intorder: int | None = None, quadrature: Tuple[ndarray, ndarray] | None = None, facets: Any | None = None, dofs: Dofs | None = None, side: int = 0, disable_doflocs: bool = False)

For integrating over facets of the mesh. Usually over the boundary.

__init__(mesh: Mesh, elem: Element, mapping: Mapping | None = None, intorder: int | None = None, quadrature: Tuple[ndarray, ndarray] | None = None, facets: Any | None = None, dofs: Dofs | None = None, side: int = 0, disable_doflocs: bool = False)

Precomputed global basis on boundary facets.

Parameters:

  • mesh – An object of type Mesh.

  • elem – An object of type Element.

  • mapping – An object of type skfem.mapping.Mapping. If None, uses mesh.mapping.

  • intorder – Optional integration order, i.e. the degree of polynomials that are integrated exactly by the used quadrature. Not used if quadrature is specified.

  • quadrature – Optional tuple of quadrature points and weights.

  • facets – Optional subset of facet indices.

  • dofs – Optional Dofs object.

  • disable_doflocs – If True, the computation of global DOF locations is disabled. This may save memory on large meshes if DOF locations are not required.

Class: InteriorFacetBasis

class skfem.assembly.InteriorFacetBasis(mesh: Mesh, elem: Element, mapping: Mapping | None = None, intorder: int | None = None, quadrature: Tuple[ndarray, ndarray] | None = None, facets: Any | None = None, dofs: Dofs | None = None, side: int = 0, disable_doflocs: bool = False)

For evaluating integrals over interior facets.

Useful for, e.g., a posteriori error estimators or implementing interior penalty/discontinuous Galerkin methods.

__init__(mesh: Mesh, elem: Element, mapping: Mapping | None = None, intorder: int | None = None, quadrature: Tuple[ndarray, ndarray] | None = None, facets: Any | None = None, dofs: Dofs | None = None, side: int = 0, disable_doflocs: bool = False)

Precomputed global basis on interior facets.

Abstract class: Form

Class: BilinearForm

class skfem.assembly.BilinearForm(form: ~typing.Callable | ~skfem.assembly.form.form.Form | None = None, dtype: ~typing.Type[~numpy.float64] | ~typing.Type[~numpy.complex64] = <class 'numpy.float64'>, nthreads: int = 0, **params)

A bilinear form for finite element assembly.

Bilinear forms are defined using functions that takes three arguments: trial function u, test function v, and a dictionary of additional parameters w.

>>> from skfem import BilinearForm, Basis, MeshTri, ElementTriP1
>>> form = BilinearForm(lambda u, v, _: u * v)
>>> form.assemble(Basis(MeshTri(), ElementTriP1())).toarray()
array([[0.08333333, 0.04166667, 0.04166667, 0.        ],
       [0.04166667, 0.16666667, 0.08333333, 0.04166667],
       [0.04166667, 0.08333333, 0.16666667, 0.04166667],
       [0.        , 0.04166667, 0.04166667, 0.08333333]])

Alternatively, you can use BilinearForm as a decorator:

>>> @BilinearForm
... def form(u, v, _):
...     return u * v

Inside the form definition, u and v are arrays containing the basis function values at quadrature points. They also contain the values of the derivatives as attributes:

>>> @BilinearForm
... def form(u, v, _):
...     # u.grad[0] is first derivative with respect to x, and so on
...     return u.grad[0] * v.grad[0] + u.grad[1] * v.grad[1]  # laplacian

We suggest you to use helper functions from skfem.helpers to make the forms more readable:

>>> from skfem.helpers import dot, grad
>>> @BilinearForm
... def form(u, v, w):
...     return dot(grad(u), grad(v))
assemble(*args, **kwargs)

Assemble the bilinear form into a sparse matrix.

Parameters:

  • ubasis – The Basis for u.

  • vbasis – Optionally, specify a different Basis for v.

  • **kwargs – Any additional keyword arguments are appended to w.

Class: LinearForm

class skfem.assembly.LinearForm(form: ~typing.Callable | ~skfem.assembly.form.form.Form | None = None, dtype: ~typing.Type[~numpy.float64] | ~typing.Type[~numpy.complex64] = <class 'numpy.float64'>, nthreads: int = 0, **params)

A linear form for finite element assembly.

Used similarly as BilinearForm with the expection that forms take two parameters v and w.

assemble(*args, **kwargs) → Any

Assemble system matrix as SciPy CSR sparse matrix.

Class: Functional

class skfem.assembly.Functional(form: ~typing.Callable | ~skfem.assembly.form.form.Form | None = None, dtype: ~typing.Type[~numpy.float64] | ~typing.Type[~numpy.complex64] = <class 'numpy.float64'>, nthreads: int = 0, **params)

A functional for postprocessing finite element solution.

Used similarly as BilinearForm with the expection that forms take one parameter w.

assemble(*args, **kwargs) → Any

Assemble system matrix as SciPy CSR sparse matrix.

elemental(v: AbstractBasis, **kwargs) → ndarray

Evaluate the functional elementwise.

Module: skfem.element

The module skfem.element defines finite elements in a very generic sense.

In order to use an element, you simply initialize the respective object and pass it to the constructor of CellBasis or FacetBasis. See below for a list of supported elements.

skfem.element.ElementH1()

\(H^1\)-conforming basis defined through a reference element.

skfem.element.ElementVector(elem[, dim])

Use the same element for each dimension.

skfem.element.ElementHdiv()

\(H(div)\)-conforming basis defined through a reference element.

skfem.element.ElementHcurl()

\(H(curl)\)-conforming basis defined through a reference element.

skfem.element.ElementGlobal()

Elements defined implicitly through global degrees-of-freedom.

skfem.element.ElementDG(elem)

Turn a finite element discontinuous by cutting connectivity.

skfem.element.ElementComposite(*elems)

Combine multiple elements.

skfem.element.ElementTriP1()

Piecewise linear element.

skfem.element.ElementTriP2()

Piecewise quadratic element.

skfem.element.ElementTriP3()

Piecewise cubic element.

skfem.element.ElementTriP4()

Piecewise quartic element.

skfem.element.ElementTriP0()

Piecewise constant element.

skfem.element.ElementTriCR()

The nonconforming Crouzeix-Raviart element.

skfem.element.ElementTriCCR

alias of ElementTriP2B

skfem.element.ElementTriRT0

alias of ElementTriRT1

skfem.element.ElementTriMorley()

The Morley element.

skfem.element.ElementTri15ParamPlate()

15-parameter nonconforming plate element.

skfem.element.ElementTriArgyris()

The Argyris element.

skfem.element.ElementTriMini

alias of ElementTriP1B

skfem.element.ElementTriHermite()

The Hermite element with 3 DOFs per vertex and one interior DOF.

skfem.element.ElementTriSkeletonP0()

Constant element for the mesh skeleton.

skfem.element.ElementTriSkeletonP1()

Linear element for the mesh skeleton.

skfem.element.ElementTriBDM1()

The lowest order Brezzi-Douglas-Marini element.

skfem.element.ElementQuad0()

Piecewise constant element.

skfem.element.ElementQuad1()

Bilinear element.

skfem.element.ElementQuad2()

Biquadratic element.

skfem.element.ElementQuadS2()

Serendipity element with 8 DOFs.

skfem.element.ElementQuadP(p)

Piecewise \(p\)'th order element.

skfem.element.ElementQuadBFS()

The Bogner-Fox-Schmit element.

skfem.element.ElementQuadRT0

alias of ElementQuadRT1

skfem.element.ElementTetP0()

Piecewise constant element.

skfem.element.ElementTetP1()

Piecewise linear element.

skfem.element.ElementTetP2()

Piecewise quadratic element.

skfem.element.ElementTetRT0

alias of ElementTetRT1

skfem.element.ElementTetN0

alias of ElementTetN1

skfem.element.ElementTetMini()

The MINI element, i.e. linear element with additional bubble DOF.

skfem.element.ElementTetCR()

Nonconforming Crouzeix-Raviart element.

skfem.element.ElementTetCCR()

Conforming Crouzeix-Raviart element.

skfem.element.ElementHex0()

Piecewise constant element.

skfem.element.ElementHex1()

Trilinear element.

skfem.element.ElementHex2()

Triquadratic element.

skfem.element.ElementHexS2()

Serendipity element with 20 DOFs.

skfem.element.ElementLineP0()

Piecewise constant element.

skfem.element.ElementLineP1()

Piecewise linear element.

skfem.element.ElementLineP2()

Piecewise quadratic element.

skfem.element.ElementLinePp(p)

Piecewise \(p\)'th order element.

skfem.element.ElementLineHermite()

\(H^2\)-conforming element with 4 DOFs.

skfem.element.ElementLineMini()

The MINI element, i.e. linear element with additional bubble DOF.

Note

The element global basis is calculated at quadrature points and stored inside DiscreteField objects. The different element types precalculate different fields of DiscreteField. E.g., for \(H(div)\) finite elements it is natural to precalculate DiscreteField.div. The high order derivatives are created only when using subclasses of ElementGlobal.

class skfem.element.DiscreteField(value=None, grad=None, div=None, curl=None, hess=None, grad3=None, grad4=None, grad5=None, grad6=None)

Module: skfem.utils

Function: solve

skfem.utils.solve(A: spmatrix, b: ndarray | spmatrix, x: ndarray | None = None, I: ndarray | None = None, solver: Callable[[...], ndarray] | Callable[[...], Tuple[ndarray, ndarray]] | None = None, **kwargs) → ndarray | Tuple[ndarray, ndarray]

Solve a linear system or a generalized eigenvalue problem.

The remaining keyword arguments are passed to the solver.

Parameters:

  • A – The system matrix

  • b – The right hand side vector or the mass matrix of a generalized eigenvalue problem.

  • solver – Choose one of the following solvers: skfem.utils.solver_direct_scipy() (default), skfem.utils.solver_eigen_scipy() (default), skfem.utils.solver_iter_pcg(), skfem.utils.solver_iter_krylov().

Function: condense

skfem.utils.condense(A: spmatrix, b: ndarray | spmatrix | None = None, x: ndarray | None = None, I: ndarray | DofsView | Dict[str, DofsView] | None = None, D: ndarray | DofsView | Dict[str, DofsView] | None = None, expand: bool = True) → spmatrix | Tuple[spmatrix, ndarray] | Tuple[spmatrix, spmatrix] | Tuple[spmatrix, ndarray, ndarray] | Tuple[spmatrix, ndarray, ndarray, ndarray] | Tuple[spmatrix, spmatrix, ndarray, ndarray]

Eliminate degrees-of-freedom from a linear system.

The user should provide the linear system A and b and either the set of DOFs to eliminate (D) or the set of DOFs to keep (I). Optionally, nonzero values for the eliminated DOFs can be supplied via x.

Note

Supports also generalized eigenvalue problems where b is a matrix.

Example

Suppose that the solution vector \(x\) can be split as

\[\begin{split}x = \begin{bmatrix} x_I\\ x_D \end{bmatrix}\end{split}\]

where \(x_D\) are known and \(x_I\) are unknown. This allows splitting the linear system as

\[\begin{split}\begin{bmatrix} A_{II} & A_{ID}\\ A_{DI} & A_{DD} \end{bmatrix} \begin{bmatrix} x_I\\ x_D \end{bmatrix} = \begin{bmatrix} b_I\\ b_D \end{bmatrix}\end{split}\]

which leads to the condensed system

\[A_{II} x_I = b_I - A_{ID} x_D.\]

As an example, let us assemble the matrix \(A\) and the vector \(b\) corresponding to the Poisson equation \(-\Delta u = 1\).

>>> import skfem as fem
>>> from skfem.models.poisson import laplace, unit_load
>>> m = fem.MeshTri().refined(2)
>>> basis = fem.CellBasis(m, fem.ElementTriP1())
>>> A = laplace.assemble(basis)
>>> b = unit_load.assemble(basis)

The condensed system is obtained with skfem.utils.condense(). Below we provide the DOFs to eliminate via the keyword argument D.

>>> AII, bI, xI, I = fem.condense(A, b, D=m.boundary_nodes())
>>> AII.toarray()
array([[ 4.,  0.,  0.,  0., -1., -1., -1., -1.,  0.],
       [ 0.,  4.,  0.,  0., -1.,  0., -1.,  0.,  0.],
       [ 0.,  0.,  4.,  0.,  0., -1.,  0., -1.,  0.],
       [ 0.,  0.,  0.,  4., -1., -1.,  0.,  0.,  0.],
       [-1., -1.,  0., -1.,  4.,  0.,  0.,  0.,  0.],
       [-1.,  0., -1., -1.,  0.,  4.,  0.,  0.,  0.],
       [-1., -1.,  0.,  0.,  0.,  0.,  4.,  0., -1.],
       [-1.,  0., -1.,  0.,  0.,  0.,  0.,  4., -1.],
       [ 0.,  0.,  0.,  0.,  0.,  0., -1., -1.,  4.]])
 >>> bI
 array([0.0625, 0.0625, 0.0625, 0.0625, 0.0625, 0.0625, 0.0625, 0.0625,
        0.0625])

By default, the eliminated DOFs are set to zero. Different values can be provided through the keyword argument x; see Example 14: Laplace with inhomogeneous boundary conditions.

Parameters:

  • A – The system matrix

  • b – The right hand side vector, or zero if x is given, or the mass matrix for generalized eigenvalue problems.

  • x – The values of the condensed degrees-of-freedom. If not given, assumed to be zero.

  • I – The set of degree-of-freedom indices to include.

  • D – The set of degree-of-freedom indices to dismiss.

  • expand – If True (default), returns also x and I. As a consequence, skfem.utils.solve() will expand the solution vector automatically.

Returns:

The condensed linear system and (optionally) information about the boundary values.

Return type:

CondensedSystem

Function: enforce

skfem.utils.enforce(A: spmatrix, b: ndarray | spmatrix | None = None, x: ndarray | None = None, I: ndarray | DofsView | Dict[str, DofsView] | None = None, D: ndarray | DofsView | Dict[str, DofsView] | None = None, diag: float = 1.0, overwrite: bool = False) → spmatrix | Tuple[spmatrix, ndarray] | Tuple[spmatrix, spmatrix]

Enforce degrees-of-freedom of a linear system.

An alternative to condense() which sets the matrix diagonals to one and right-hand side vector to the enforced degree-of-freedom value.

Note

The original system is both returned (for compatibility with skfem.utils.solve()) and optionally (if overwrite) modified (for performance).

Parameters:

  • A – The system matrix

  • b – Optionally, the right hand side vector.

  • x – The values of the enforced degrees-of-freedom. If not given, assumed to be zero.

  • I – Specify either this or D: The set of degree-of-freedom indices to solve for.

  • D – Specify either this or I: The set of degree-of-freedom indices to enforce (rows/diagonal set to zero/one).

  • overwrite – Optionally, the original system is both modified (for performance) and returned (for compatibility with skfem.utils.solve()). By default, False.

Returns:

A linear system with the enforced rows/diagonals set to zero/one.

Return type:

LinearSystem

Module: skfem.helpers

Helper functions for defining forms.

skfem.helpers.grad(u: DiscreteField)

Gradient.

skfem.helpers.div(u: DiscreteField)

Divergence.

skfem.helpers.curl(u: DiscreteField)

Curl.

skfem.helpers.d(u: DiscreteField)

Gradient, divergence or curl.

skfem.helpers.dd(u: DiscreteField)

Hessian (for ElementGlobal).

skfem.helpers.ddd(u: DiscreteField)

Third derivative (for ElementGlobal).

skfem.helpers.dddd(u: DiscreteField)

Fourth derivative (for ElementGlobal).

skfem.helpers.sym_grad(u: DiscreteField)

Symmetric gradient.

skfem.helpers.dot(u: DiscreteField | ndarray, v: DiscreteField | ndarray)

Dot product.

skfem.helpers.ddot(u: DiscreteField | ndarray, v: DiscreteField | ndarray)

Double dot product.

skfem.helpers.dddot(u: DiscreteField | ndarray, v: DiscreteField | ndarray)

Triple dot product.

skfem.helpers.mul(A: DiscreteField | ndarray, x: DiscreteField | ndarray)

Matrix multiplication.

skfem.helpers.trace(T)

Trace of matrix.

skfem.helpers.transpose(T)

Transpose of matrix.

skfem.helpers.prod(u: DiscreteField | ndarray, v: DiscreteField | ndarray, w: DiscreteField | ndarray | None = None)

Tensor product.

skfem.helpers.inv(A)

Inverse of an array A over trailing axis (if any).

Module: skfem.visuals

skfem.visuals.matplotlib.plot(m, u, **kwargs) → Axes
skfem.visuals.matplotlib.plot(m: MeshLine1, z: ndarray, **kwargs)
skfem.visuals.matplotlib.plot(m: MeshTri1, z: ndarray, **kwargs) → Axes
skfem.visuals.matplotlib.plot(m: MeshQuad1, z, **kwargs)
skfem.visuals.matplotlib.plot(basis: CellBasis, z: ndarray, **kwargs) → Axes

Plot functions defined on nodes of the mesh.