API reference

Contents

API reference#

Mesh utilities#

scifem.create_entity_markers(domain: dolfinx.mesh.Mesh, dim: int, entities_list: list[TaggedEntities]) → dolfinx.mesh.MeshTags[source]#

Mark entities of specified dimension according to a geometrical marker function.

Parameters:
  • domain – A dolfinx.mesh.Mesh object

  • dim – Dimension of the entities to mark

  • entities_list –

    A list of tuples with the following elements:

    • index 0: The tag to assign to the entities

    • index 1: A function that takes a point and returns a boolean array indicating whether the point is inside the entity

    • index 2: Optional, if True, the entities will be marked on the boundary

Returns:

A dolfinx.mesh.MeshTags object with the corresponding entities marked. If an entity satisfies multiple input marker functions, it is not deterministic what value the entity gets.

Note:

scifem.reverse_mark_entities(entity_map: IndexMap, entities: ndarray[tuple[Any, ...], dtype[int32]]) → ndarray[tuple[Any, ...], dtype[int32]][source]#

Communicate entities marked on a single process to all processes that ghosts or owns this entity.

Parameters:
  • entity_map – Index-map describing entity ownership

  • entities – Local indices of entities to communicate

Returns:

Local indices marked on any process sharing this entity

scifem.extract_submesh(mesh: Mesh, entity_tag: MeshTags, tags: Sequence[int]) → SubmeshData[source]#

Generate a sub-mesh from a subset of tagged entities in a meshtag object.

Parameters:
  • mesh – The mesh to extract the submesh from.

  • entity_tag – MeshTags object containing marked entities.

  • tags – What tags the marked entities used in the submesh should have.

Returns:

A tuple (submesh, subcell_to_parent_entity, subvertex_to_parent_vertex, subnode_to_parent_node, entity_tag_on_submesh).

scifem.transfer_meshtags_to_submesh(entity_tag: MeshTags, submesh: Mesh, *, vertex_to_parent: _EntityMap | ndarray[tuple[Any, ...], dtype[int32]], cell_to_parent: _EntityMap | ndarray[tuple[Any, ...], dtype[int32]]) → tuple[MeshTags, ndarray[tuple[Any, ...], dtype[int32]]][source]#

Transfer a entity_tag from a parent mesh to a submesh.

Parameters:
  • entity_tag – Tag to transfer

  • submesh – Submesh to transfer tag to

  • vertex_to_parent – Mapping from submesh vertices to parent mesh vertices

  • cell_to_parent – Mapping from submesh cells to parent entities

Returns:

submesh_tag is the tag on the submesh and sub_to_parent_entity_map is a mapping from submesh entities in the tag to the corresponding entities in the parent.

Return type:

A tuple (submesh_tag, sub_to_parent_entity_map) where

scifem.find_interface(cell_tags: MeshTags, id_0: tuple[int, ...], id_1: tuple[int, ...]) → ndarray[tuple[Any, ...], dtype[int32]][source]#

Given to sets of cells, find the facets that are shared between them.

Parameters:
  • cell_tags – MeshTags object marking cells.

  • id_0 – Tags to extract for domain 0

  • id_1 – Tags to extract for domain 1

Returns:

The facets shared between the two domains.

scifem.compute_interface_data(cell_tags: MeshTags, facet_indices: ndarray[tuple[Any, ...], dtype[int32]], include_ghosts: bool = False) → ndarray[tuple[Any, ...], dtype[int32]][source]#

Compute interior facet integrals that are consistently ordered according to the cell_tags, such that the data (cell0, facet_idx0, cell1, facet_idx1) is ordered such that cell_tags[cell0]`<`cell_tags[cell1], i.e the cell with the lowest cell marker is considered the “+” restriction”.

Parameters:
  • cell_tags – MeshTags that must contain an integer marker for all cells adjacent to the facet_indices

  • facet_indices – List of facets (local index) that are on the interface.

  • include_ghosts – If True integration entities will include facets that are ghosts on

  • facets. (the current process. This is for instance useful for interpolation on interior)

Returns:

The integration data.

scifem.compute_subdomain_exterior_facets(mesh: Mesh, ct: MeshTags, markers: Sequence[int]) → ndarray[tuple[Any, ...], dtype[int32]][source]#

Find the the facets that are considered to be on the “exterior” boundary of a subdomain.

The subdomain is defined as the collection of cells in ct that is marked with any of the markers. The exterior boundary of the subdomain is defined as the collection of facets that are only connected to a single cell within the subdomain.

Note

Ghosted facets are included in the resulting array.

Parameters:
  • mesh – Mesh to extract subdomains from

  • ct – MeshTags object marking subdomains

  • markers – The tags making up the “new” mesh

Returns:

The exterior facets

Function spaces and degrees of freedom#

scifem.create_real_functionspace(mesh: Mesh, value_shape: tuple[int, ...] = ()) → FunctionSpace[source]#

Create a real function space.

Parameters:
  • mesh – The mesh the real space is defined on.

  • value_shape – The shape of the values in the real space.

Returns:

The real valued function space.

Note

For scalar elements value shape is ().

scifem.create_space_of_simple_functions(mesh: Mesh, cell_tag: MeshTags, tags: Sequence[int | int32 | int64], value_shape: tuple[()] | tuple[int | tuple[int]] | None = None) → FunctionSpace[source]#

Create a space of simple functions.

This is a space that represents piecewise constant functions of N patches, where N is the number of input tags. Each patch is defined by the cells in cell_tag which is marked with the corresponding tag value.

Note

All cells are expected to have a tag in tags.

Parameters:
  • mesh – The mesh the function space is defined on.

  • cell_tag – The mesh tags defining the different patches of cells.

  • tags – The set of unique values within cell_tag that defines the different patches.

  • value_shape – The shape of the values in the space.

Returns:

The space of simple functions.

scifem.vertex_to_dofmap(V: FunctionSpace) → ndarray[tuple[Any, ...], dtype[int32]][source]#

Create a map from the vertices (local to the process) to the correspondning degrees of freedom.

Parameters:

V – The function space

Returns:

An array mapping local vertex i to local degree of freedom

Note

If using a blocked space this map is not unrolled for the DofMap block size.

scifem.dof_to_vertexmap(V: FunctionSpace) → ndarray[tuple[Any, ...], dtype[int32]][source]#

Create a map from the degrees of freedom to the vertices of the mesh. As not every degree of freedom is associated with a vertex, every dof that is not associated with a vertex returns -1

Parameters:

V – The function space

Returns:

An array mapping local dof i to a local vertex

Assembly#

scifem.assemble_scalar(J: Form | Form, entity_maps: dict[Mesh, ndarray[tuple[Any, ...], dtype[int32]]] | None = None) → floating | complexfloating[source]#

Assemble a scalar form and gather result across processes

Parameters:
  • form – The form to assemble.

  • entity_maps – Maps of entities on related submeshes to the domain used in J.

Returns:

The accumulated value of the assembled form.

scifem.norm(expr: Expr, norm_type: Literal['L2', 'H1', 'H10'], entity_maps: dict[Mesh, ndarray[tuple[Any, ...], dtype[int32]]] | None = None) → Form[source]#

Compile the norm of an UFL expression into a DOLFINx form.

Parameters:
  • expr – UFL expression

  • norm_type – Type of norm

  • entity_maps – Mapping for Constants and Coefficients within the expression that lives on a submesh.

Boundary conditions and sources#

scifem.interpolate_function_onto_facet_dofs(Q: FunctionSpace, expr: Expr, facets: ndarray[tuple[Any, ...], dtype[int32]]) → Function[source]#

Create a function \(u_h\in Q\) such that \(u_h=\text{expr}\) for all dofs belonging to a subset of facets. All other dofs are set to zero.

Note

The resulting function is only correct in the “normal” direction, i.e. \(u_{bc}\cdot n = expr\), while the tangential component is uncontrolled. This makes it hard to visualize the function when outputting it to file, either through interpolation to an appropriate DG space, or to a point-cloud.

Parameters:
  • Q – The function space to create the function $u_h$ in.

  • expr – The expression to evaluate.

  • facets – The facets on which to evaluate the expression.

class scifem.PointSource(V: FunctionSpace, points: ndarray[tuple[Any, ...], dtype[float32]] | ndarray[tuple[Any, ...], dtype[float64]], magnitude: floating | complexfloating = np.float64(1.0), tol: floating | None = None)[source]#

Class for defining a point source in a given function space.

apply_to_vector(b: Function | Vector | Vec, recompute: bool = False)[source]#

Apply the point sources to a vector.

Parameters:
  • b – The vector to apply the point sources to.

  • recompute – If the point sources should be recomputed before applying. Recomputation should only be done if the mesh geometry has been modified.

Note

The user is responsible for forward scattering of the vector after applying the point sources.

Note

If a PETSc vector is passed in, one has to call b.assemble() prior to solving the linear system (post scattering).

compute_cell_contributions()[source]#

Compute the basis function values at the point sources.

recompute_sources(tol: float | None = None)[source]#

Recompute the what cells the point sources collide with.

This function should be called if the mesh geometry has been modified.

Parameters:

tol – Tolerance for point location. If None, the tolerance from initialization is used.

Interpolation#

scifem.interpolation_matrix(expr: Expr, Q: FunctionSpace) → MatrixCSR[source]#

Create the interpolation matrix \(\Lambda\) of a UFL-expression such that

\[\begin{split}\begin{align*} \Lambda: V &\rightarrow Q \\ \Lambda u &= \sum_{i=0}^{N_Q-1}\sum_{j=0}^{N_V-1} \phi_i l_i(expr(\psi_j))u_j \end{align*}\end{split}\]

where \(l_j\) is the dual basis of the space \(Q\) with basis functions \(\phi_j\), and \(\psi_j\) are the basis functions of the space \(V\).

Parameters:
  • expr – The UFL expression

  • Q – Output interpolation space

Returns:

Interpolation matrix as a MatrixCSR.

scifem.prepare_interpolation_data(expr: Expr, Q: FunctionSpace, interpolation_entities: ndarray[tuple[Any, ...], dtype[int32]] | None = None) → ndarray[tuple[Any, ...], dtype[inexact]][source]#

Convenience function for preparing data required for assembling the interpolation matrix

\[\begin{split}\begin{align*} \Lambda: V &\rightarrow Q \\ \Lambda u &= \sum_{i=0}^{N_Q-1}\sum_{j=0}^{N_V-1} \phi_i l_i(expr(\psi_j))u_j \end{align*}\end{split}\]

where \(l_j\) is the dual basis of the space \(Q\) with basis functions \(\phi_j\), and \(\psi_j\) are the basis functions of the space \(V\).

Parameters:
  • expr – The UFL expression containing a trial function from space V

  • Q – Output interpolation space

  • interpolation_entities – Entities of the domain of the input space V that one should evaluate the expr at. If not provided, it is assumed that we are integrating over all cells in V and that Q is defined on the same grid.

Returns:

Interpolation data per cell, as an numpy array.

scifem.petsc_interpolation_matrix(expr: Expr, Q: FunctionSpace, use_petsc: bool = False) → Mat[source]#

Create the interpolation matrix \(\Lambda\) of a UFL-expression such that

\[\begin{split}\begin{align*} \Lambda: V &\rightarrow Q \\ \Lambda u &= \sum_{i=0}^{N_Q-1}\sum_{j=0}^{N_V-1} \phi_i l_i(expr(\psi_j))u_j \end{align*}\end{split}\]

where \(l_j\) is the dual basis of the space \(Q\) with basis functions \(\phi_j\), and \(\psi_j\) are the basis functions of the space \(V\).

Parameters:
  • expr – The UFL expression

  • Q – Output interpolation space

Returns:

Interpolation matrix as a PETSc.Mat.

Facet submeshes#

Moving data between a mesh and a submesh of its facets, in both directions. The classes prepare the Expression, the connectivities and the facet orientations once, for repeated use; the functions are one-off calls to them.

scifem.interpolation.interpolate_to_surface_submesh(u_volume: Function, u_surface: Function, submesh_facets: ndarray[tuple[Any, ...], dtype[int32]], integration_entities: ndarray[tuple[Any, ...], dtype[int32]], entity_maps: list[_EntityMap] | None = None)[source]#

Interpolate a function u_volume into the function u_surface.

See SurfaceSubmeshInterpolation, which prepares the interpolation once for repeated use.

Note

Does not work for DG as no dofs are associated with the facets in versions of DOLFINx prior to FEniCS/dolfinx#4140, which is included in version 0.11.0 and later.

Parameters:
  • u_volume – Function to interpolate data from

  • u_surface – Function to interpolate data to

  • submesh_facets – Cells in facet mesh

  • integration_entities – Integration entities on the parent mesh corresponding to the facets in submesh_facets

  • entity_maps – Entity maps for an expression with coefficients on other meshes

class scifem.interpolation.SurfaceSubmeshInterpolation(expr: Expr, V_surface: FunctionSpace, submesh_facets: ndarray[tuple[Any, ...], dtype[int32]], integration_entities: ndarray[tuple[Any, ...], dtype[int32]], entity_maps: list[_EntityMap] | None = None)[source]#

Interpolation of a volume expression into a space on a facet submesh, prepared once.

The expression is compiled into a dolfinx.fem.Expression at the surface element’s interpolation points, and the connectivities and permutation info it needs are created, once. apply() then evaluates the expression on the parent facets, at the current values of its coefficients, and interpolates the result into each submesh cell. For DOLFINx < 0.11 the values also have to be reordered from each facet’s orientation in its cell to its own (compute_entity_closure_permutations()); that permutation is computed here too.

Note

Does not work for DG as no dofs are associated with the facets in versions of DOLFINx prior to FEniCS/dolfinx#4140, which is included in version 0.11.0 and later.

Parameters:
  • expr – Function or expression on the parent mesh to interpolate from

  • V_surface – Space on the facet submesh to interpolate into

  • submesh_facets – Cells in facet mesh

  • integration_entities – Integration entities on the parent mesh corresponding to the facets in submesh_facets

  • entity_maps – Entity maps for an expression with coefficients on other meshes

property V_surface: FunctionSpace#

The space on the facet submesh interpolated into.

apply(u_surface: Function)[source]#

Interpolate the expression, at its coefficients’ current values, into u_surface.

Parameters:

u_surface – Function in V_surface, overwritten on the submesh cells, and scattered forward.

property expression: Expression#

The compiled expression, at the surface element’s interpolation points.

scifem.interpolation.interpolate_from_surface_submesh(u_surface: Function, u_volume: Function, submesh_facets: ndarray[tuple[Any, ...], dtype[int32]], integration_entities: ndarray[tuple[Any, ...], dtype[int32]], entity_maps: list[_EntityMap] | None = None)[source]#

Extend a function u_surface on a facet submesh by zero into the function u_volume.

The reverse of interpolate_to_surface_submesh(): u_volume takes the values of u_surface at its nodes on the facets and is zero at every other node. See SurfaceSubmeshExtension, which can be reused across calls.

Parameters:
  • u_surface – Function on the facet submesh to extend

  • u_volume – Function on the parent mesh to extend into, continuous Lagrange or Piola-mapped

  • submesh_facets – Cells in facet mesh

  • integration_entities – Integration entities on the parent mesh corresponding to the facets in submesh_facets

  • entity_maps – The submesh’s entity map, needed for a Piola-mapped u_volume

class scifem.interpolation.SurfaceSubmeshExtension(V_surface: FunctionSpace, V_volume: FunctionSpace, submesh_facets: ndarray[tuple[Any, ...], dtype[int32]], integration_entities: ndarray[tuple[Any, ...], dtype[int32]], entity_maps: list[_EntityMap] | None = None)[source]#

Extension by zero of a function on a facet submesh into a space on its parent mesh.

The reverse of SurfaceSubmeshInterpolation: the dofs on the closure of each submesh facet are set from the surface function and all other volume dofs are zero. Dofs shared by several facets get the mean of the facets’ values (see apply()).

Supported volume spaces are continuous Lagrange and, given entity_maps, Piola-mapped spaces such as RT and N1curl, of which only the facet trace (the normal or tangential component) is set.

For a surface space that is the trace of the volume space, this is a right inverse of interpolate_to_surface_submesh().

Parameters:
  • V_surface – A space on the facet submesh: with as many components as V_volume for Lagrange, and vector-valued of the geometric dimension otherwise.

  • V_volume – A continuous space on the parent mesh.

  • submesh_facets – Cells of the submesh to extend from.

  • integration_entities – (cell, local facet) on the parent mesh of each of submesh_facets, as for interpolate_to_surface_submesh().

  • entity_maps – The submesh’s entity map, needed for a volume space that is not Lagrange.

Raises:
  • ValueError – If V_volume is discontinuous, the value sizes do not fit, or entity_maps is missing where it is needed.

  • NotImplementedError – If V_surface needs dof transformations, the volume element’s pullback cannot be inverted, or its interpolation points differ between the facets of a cell.

property V_surface: FunctionSpace#

The surface space.

property V_volume: FunctionSpace#

The volume space.

apply(u_surface: Function, u_volume: Function)[source]#

Set u_volume to the extension of u_surface.

Unlike interpolate_to_surface_submesh() and scifem.interpolate_function_onto_facet_dofs(), this does not call dolfinx.fem.Function.interpolate(), which would set every dof of each parent cell and zero boundary nodes that lie on none of that cell’s facets. The values are instead added straight into the volume dofs of each facet’s closure, scaled by weights, then summed onto their owners with a reverse scatter, i.e. u_volume[volume_dofs] += weights * (basis @ u_surface[surface_dofs]). Each boundary dof thus gets the mean of its writes: its value for a compatible (continuous) surface space, and an average of the facets’ values where they disagree.

Parameters:
  • u_surface – Function in V_surface, with consistent ghosts.

  • u_volume – Function in V_volume, overwritten, ghosts included.

property basis: ndarray[tuple[Any, ...], dtype[floating]]#

The map from each facet’s surface_dofs to its volume_dofs, (facets, volume dofs, surface dofs).

property surface_dofs: ndarray[tuple[Any, ...], dtype[int32]]#

The surface dofs of each facet’s submesh cell, unrolled by block, (facets, surface dofs).

property volume_dofs: ndarray[tuple[Any, ...], dtype[int32]]#

The volume dofs of each facet’s closure, local to the process and unrolled by block, (facets, volume dofs).

property weights: ndarray[tuple[Any, ...], dtype[floating]]#

One over the number of writes to each of volume_dofs, over all facets and processes, with the same shape.

Both directions line up the dofs of a sub-entity as seen from its cell with those of the entity taken as a cell of its own, which these compute:

scifem.interpolation.compute_entity_closure_permutations(V: FunctionSpace, dim: int, entities: ndarray[tuple[Any, ...], dtype[int32]], size: int | None = None) → ndarray[tuple[Any, ...], dtype[int32]][source]#

Permutations taking each entity’s closure dofs from its cell’s orientation to its own.

Row i reorders data given at the closure dofs of the local entity entities[i, 1] of cell entities[i, 0], in the order of that entity in the cell, into the order of the entity as a cell of its own (as in a submesh), with its vertices ordered by global index: data_entity = data_cell[perm[i]]. Every cell sharing an entity therefore agrees on the order. It is basix’s permute_subentity_closure_inv with the cell’s permutation info, computed once per distinct (cell permutation info, local entity) pair.

Parameters:
  • V – Space on the parent mesh whose element defines the permutations.

  • dim – Topological dimension of the entities, below that of the cells.

  • entities – (cell, local entity) pairs, shape (num_entities, 2), such as facet integration entities for dim = tdim - 1.

  • size – Length of each permutation. Defaults to the number of closure dofs of an entity.

Returns:

The permutations, shape (num_entities, size).

Raises:

ValueError – If size is not given and the entities’ closures differ in size, as for the triangular and quadrilateral facets of a prism.

scifem.interpolation.compute_entity_closure_dofs(V: FunctionSpace, dim: int, entities: ndarray[tuple[Any, ...], dtype[int32]]) → ndarray[tuple[Any, ...], dtype[int32]][source]#

The dofs of each entity’s closure, from its cell’s dofmap, in the entity’s own orientation.

Entry [i, j] is the dof of V (a node, for a blocked space) at the j-th closure dof of the local entity entities[i, 1] of cell entities[i, 0], taken as a cell of its own. Cells sharing an entity give the same row, and for facets it lines up with the dofs of the same element’s trace on a facet submesh. See compute_entity_closure_permutations().

Parameters:
  • V – Space on the parent mesh.

  • dim – Topological dimension of the entities, below that of the cells.

  • entities – (cell, local entity) pairs, shape (num_entities, 2).

Returns:

The dofs, shape (num_entities, num closure dofs).

Evaluation and geometry#

scifem.evaluate_function(u: Function, points: Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str], broadcast=True) → ndarray[tuple[Any, ...], dtype[float64]][source]#

Evaluate a function at a set of points.

Parameters:
  • u – The function to evaluate.

  • points – The points to evaluate the function at.

  • broadcast –

    If True, the values will be broadcasted to all processes.

    Note

    Uses a global MPI call to broadcast values, thus this has to be called on all active processes synchronously.

    Note

    If the function is discontinuous, different processes may return different values for the same point. In this case, the value returned is the maximum value across all processes.

Returns:

The values of the function evaluated at the points.

scifem.find_cell_extrema(u: Expr, cell: int, kind: Callable[[Sequence[T]], T], x0: ndarray[tuple[Any, ...], dtype[_ScalarT]] | None = None, jit_options: dict | None = None, method: str | None = None, options: dict | None = None, tol: float | None = None) → tuple[ndarray[tuple[Any, ...], dtype[floating]], floating][source]#

Find the extrema of a ufl.core.expr.Expr within a cell.

Parameters:
  • u – The expression to find the extrema of.

  • cell – The local index of the cell to search in

  • kind – If we search for minima or maxima

  • x_0 – The point to start the initial search at

  • method – Optimization algorithm to use for local problem

  • options – Options for optimization method

  • tol – Tolerance for scipy minimize

Returns:

The point (in physical space) of the extrema and the value at the extrema.

scifem.compute_extrema(u: Expr, kind=typing.Callable[[typing.Sequence[~T]], ~T], p: int = 1, num_candidates: int = 3, x0: ndarray[tuple[Any, ...], dtype[floating]] | None = None, method: str | None = None, options: dict | None = None, tol: float | None = None, jit_options: dict[str, Any] | None = None) → tuple[floating, ndarray[tuple[Any, ...], dtype[floating]]][source]#

Find the extrema of an expression across its integration domain.

Parameters:
  • u – The expression

  • p – Integer to compute the cell-wise L^p norm to speed up search.

  • kind – min or max.

  • num_candidates – The number of candidates with the largert average $L^p$ norm over a cell.

  • x0 – Initial point in reference cell to start search at.

  • method – Optimization algorithm to use for local problem

  • options – Options for optimization method

  • tol – Tolerance for scipy minimize

Returns:

The value at the extrema and the physical point

scifem.closest_point_projection(mesh: Mesh, cells: ndarray[tuple[Any, ...], dtype[int32]], target_points: ndarray[tuple[Any, ...], dtype[float64 | float32]], tol_x: float | None = None, tol_dist: float = 1e-10, tol_grad: float = 1e-10, max_iter: int = 2000, max_ls_iter: int = 250, num_threads: int = 1)[source]#

Projects a 3D point onto a cell in a potentially higher order mesh.

Uses the Goldstein-Levitin-Polyak Gradient projection method, where potential simplex constraints are handled by an exact projection using a primal-dual root finding method. See: - Held, M., Wolfe, P., Crowder, H.: Validation of subgradient optimization (1974) - Laurent Condat. Fast Projection onto the Simplex and the l1 Ball. (2016) - Dimitri P. Bertsekas, “On the Goldstein-Levitin-Polyak gradient projection method,” (1976)

Parameters:
  • mesh – The mesh containing the cells.

  • cells – The local indices of the cells to project onto.

  • target_points – (n, 3) numpy array, the 3D points to project.

  • tol_x – Tolerance for changes between iterates in the reference coordinates. If None, uses the square root of machine precision.

  • tol_grad – Tolerance for determination based on the gradient norm. The gradient is scaled by the Jacobian to account for stretching.

  • tol_dist – Tolerance used to determine if the projected point is close enough to the target point to stop optimization.

  • max_iter – int, the maximum number of iterations for the projected gradient method.

  • max_ls_iter – int, the maximum number of line search iterations.

  • num_threads – int, the number of threads to use for parallel projection.

Returns:

A tuple of arrays containing the closest points (in physical space) and reference coordinates for each cell to each target point.

Solvers#

class scifem.NewtonSolver(**kwargs)[source]#
A: Mat#
property F#

The list of residuals where each entry is a dolfinx.fem.Form.

property J#

The Jacobian blocks represented as lists of lists where each entry is a dolfinx.fem.Form.

b: Vec#
max_iterations: int#
set_post_solve_callback(callback: Callable[[NewtonSolver], None])[source]#

Set a callback function that is called after each Newton iteration.

set_pre_solve_callback(callback: Callable[[NewtonSolver], None])[source]#

Set a callback function that is called before each Newton iteration.

solve(atol=1e-06, rtol=1e-08, beta=1.0) → int[source]#

Solve the nonlinear problem using Newton’s method.

Parameters:
  • atol – Absolute tolerance for the update.

  • rtol – Relative tolerance for the update.

  • beta – Damping parameter for the update.

Returns:

The number of Newton iterations used to converge.

Note

The tolerance is on the 0-norm of the update.

Periodic meshes#

Building a periodic mesh, moving data on and off it, and the vertex correspondence the rebuild runs on – found either from the coordinates or from the $Periodic section of a gmsh model.

class scifem.periodic.PeriodicNodes(replaced: ~numpy.ndarray[tuple[~typing.Any, ...], ~numpy.dtype[~numpy.int64]] = <factory>, partner: ~numpy.ndarray[tuple[~typing.Any, ...], ~numpy.dtype[~numpy.int64]] = <factory>, num_nodes_global: int = 0)[source]#

Node pairs to identify, resolved to roots.

The pairs are given in the mesh’s input global numbering, so they say nothing about how the mesh is distributed and can come from anywhere that knows it – scifem.periodic.extract_gmsh_periodic_nodes() reads them out of a $Periodic section, but nothing here depends on that.

Only the process that has the pairs holds them; every other one passes an empty set, which is what the defaults are for – PeriodicNodes() says “nothing here” without the caller having to build two empty arrays to say it.

Parameters:
  • replaced – 0-based node indices that are to be replaced, ascending and without repeats. These are values of input_global_indices.

  • partner – For each entry of replaced, the node it is identified with. Never itself replaced, so no further resolution is needed.

  • num_nodes_global – The size of the input global numbering, i.e. one past its largest index. Not mesh.geometry.index_map().size_global, which is smaller whenever the mesh was built from a node set with entries no cell references. Taken from root and broadcast, so the default stands on every other process; on root it has to be set, and periodic_correspondence_from_nodes() checks that it was.

class scifem.periodic.VertexCorrespondence(indicator_vertices: ndarray[tuple[Any, ...], dtype[int32]], indicator_facets: ndarray[tuple[Any, ...], dtype[int32]], src_owner: ndarray[tuple[Any, ...], dtype[int32]], dest_owner: ndarray[tuple[Any, ...], dtype[int32]], partner_vertex: ndarray[tuple[Any, ...], dtype[int32]])[source]#

Which vertices of mesh are identified with which, and which ranks hold each end.

This is what scifem.periodic rebuilds from: it consumes nothing else and never evaluates a coordinate, which is what lets the pairs be found either geometrically or topologically.

Stores the data of dolfinx.geometry.PointOwnershipData for the partner_vertex, over query points that are the images of indicator_vertices – the vertices given up to the partner side – plus one extra array, indicator_facets, the facets given up with them.

The two halves are keyed independently, so indicator_vertices[i] is not the vertex replaced by partner_vertex[i]: indicator_vertices and src_owner are keyed on what this process gives up, dest_owner and partner_vertex on what other processes gave up to it, and the two sides of a pair rarely live on the same process. Splitting a 6x6 unit square over three ranks gives one rank 7 and 0, and another 0 and 8. In serial the two lengths coincide, which makes the assumption easy to form and wrong to act on. They line up only after the exchange, which is what compute_insert_position reorders.

Parameters:
  • indicator_vertices – Local vertices, owned and ghost, that are to be replaced by their partner vertex. Broadened across processes: a vertex marked on its owner is marked on every process that ghosts it.

  • indicator_facets – Local facets, owned and ghost, lying on the seam, i.e. the exterior facets all of whose vertices are in indicator_vertices. Broadened the same way.

  • src_owner – For each entry of indicator_vertices, the rank owning one of the cells its partner vertex belongs to. A vertex is shared by several cells, so which one this names is arbitrary – _build_periodic_mesh recovers the rest of the ranks holding a cell at that vertex, which all need the seam cells too. Note also that this is a cell owner, not the owner of the partner vertex, which the rank in question may merely ghost.

  • dest_owner – For each vertex this process is the far side of, the rank that asked. Must be sorted ascending: the packing groups by destination and relies on it.

  • partner_vertex – For each entry of dest_owner, the local vertex that replaces the vertex that rank gave up. Same length as dest_owner.

scifem.periodic.check_cells_stayed_distinct(mesh)[source]#

Raise if two cells of mesh carry the same vertices.

This is what going too coarse across a periodic direction does. With two cells between the two sides of a seam, a cell’s opposite facets are identified with each other and the two cells end up on the same vertices – they are one topological cell drawn twice. Nothing else notices: no vertex repeats within a cell, and the vertex count, cell count and volume are all what a correct mesh would have. Measured on a 2x2 doubly periodic square, all four cells collapse onto a single vertex set.

Not a manifold test, which would be easier and wrong. A facet with more than two cells is also what a legitimately non-manifold mesh has – three sheets meeting along an edge, the stem of a T – and refusing that would refuse a geometry nobody asked us to object to. What separates the two is whether any of those cells are the same cell, so that is what is asked, directly.

Ghost cells are included, which is what lets a duplicate pair split across two processes be seen: the two share every facet, so each is ghosted onto the other’s owner.

Collective, and a postcondition – call it on the rebuilt mesh, not the input.

Parameters:

mesh – The mesh to check. One cell type, as everywhere else here: the dofmap is read as a rectangular array.

Raises:

RuntimeError – If two distinct cells share a vertex set.

scifem.periodic.check_facet_ghosting(mesh)[source]#

Raise unless every facet between two processes carries both of its cells.

The rebuild assumes it: a cell is shipped across the seam so that the facet it will be glued along ends up with two incident cells, and if the mesh already fails that away from the seam, the result is a mesh whose interior facet integrals cannot be assembled and whose diagnosis points at the seam rather than at the input.

interprocess_facets is what makes this checkable. It names the facets on the interprocess boundary from the facet index map, so it is the same set whatever the ghost mode – unlike a local “one incident cell” test, which cannot tell an exterior facet from an unghosted interprocess one. Under shared_facet every such facet has two incident cells; under none every one of them has one.

Collective, and a no-op in serial, where there are no interprocess facets.

Parameters:

mesh – The mesh to check.

Raises:

RuntimeError – If any interprocess facet has other than two incident cells.

scifem.periodic.create_periodic_mesh(mesh, indicator, mapping_function, tag_base: int = 1101) → tuple[Mesh, ndarray[tuple[Any, ...], dtype[int32]], ndarray[tuple[Any, ...], dtype[int32]]][source]#

Create a periodic mesh that takes all facets that satisfy the indicator function, and map the vertices of these facets to the vertices that satisfies the mapping function.

Note

The cell ownership does not change, only additional ghosts are added to a given process

Note

The vertex ownership does not change, only additional ghosts are added to a given process

Note

This is match_vertices_geometric() followed by _build_periodic_mesh. Only the first half evaluates indicator and mapping_function; a reader that knows the vertex pairs already, such as one for the $Periodic section of a gmsh file, builds a VertexCorrespondence and calls the second half directly.

Parameters:
  • mesh – The mesh to make periodic. It has to carry a layer of ghost cells across every interprocess facet; see check_facet_ghosting().

  • indicator – Marks the entities to be replaced, given coordinates as (3, n).

  • mapping_function – Maps a marked vertex to the one it is identified with, given coordinates as (3, n).

  • tag_base – The first of scifem.periodic.mesh.NUM_CONSENSUS_TAGS consecutive MPI tags for the consensus exchanges inside the rebuild. Only worth setting when another such exchange can be in flight on an overlapping communicator at the same time – a nested call, or two sub-communicators that share ranks – since the tags would then have to be spaced apart.

Returns:

A tuple (new_mesh, replaced_vertices, replacement_map) where new_mesh is the new mesh with periodicity, replaced_vertices is a list of vertices of the input mesh that has been replaced (local to process). replacement_map is a map from the old vertices (local to process) to the new vertices (local to process).

Note

This map does not contain additional ghost vertices added to the process that has taken over the facet or given away a facet.

Example

mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 7, 19)

def indicator(x):
    return numpy.isclose(x[1], 1)

def map(x):
    values = x.copy()
    values[1] -= 1
    return values

periodic_mesh, _, _ = create_periodic_mesh(mesh, indicator, map)
scifem.periodic.create_periodic_mesh_from_igi(mesh, replaced_igi, partner_igi, num_nodes_global, root: int = 0, tag_base: int = 1101) → tuple[Mesh, ndarray[tuple[Any, ...], dtype[int32]], ndarray[tuple[Any, ...], dtype[int32]]][source]#

Make mesh periodic from the node pairs .

The point of VertexCorrespondence is that it is the seam between finding the periodic pairs and rebuilding the mesh from them. Everything geometric – the indicator, the mapping function, the tolerance, the point searches – lives on the match_vertices_geometric() side of it, and _build_periodic_mesh sees only the struct. So a reader that already knows the pairing, as gmsh does, fills the same fields and reuses the rebuild unchanged: no indicator, no mapping_function, and therefore no tolerance to tune and no risk of a snap onto the wrong vertex. It also handles rotational and reflective periodicity, which the coordinate mapping can only express if the caller writes the transform by hand.

Collective.

Parameters:
  • mesh – The mesh read from the same gmsh model, so that mesh.geometry.input_global_indices is the node numbering the pairs use.

  • replaced_igi – Corresponding node pairs, as 0-based gmsh node tags. Held on root only; ignored elsewhere. Every partner must be a root – a node that is not itself a replaced – so chains through a corner have to be resolved first, which scifem.periodic.extract_gmsh_periodic_nodes() does for a gmsh model.

  • partner_igi – Corresponding node pairs, as 0-based gmsh node tags. Held on root only; ignored elsewhere. Every partner must be a root – a node that is not itself a replaced – so chains through a corner have to be resolved first, which scifem.periodic.extract_gmsh_periodic_nodes() does for a gmsh model.

  • num_nodes_global – The number of nodes in the gmsh model. Not mesh.geometry.index_map().size_global, which is smaller when create_mesh drops nodes no cell references.

  • root – The rank holding the pairs.

  • tag_base – The first of scifem.periodic.mesh.NUM_CONSENSUS_TAGS consecutive MPI tags for the consensus exchanges inside the rebuild. Only worth setting when another such exchange can be in flight on an overlapping communicator at the same time – a nested call, or two sub-communicators that share ranks – since the tags would then have to be spaced apart.

Returns:

As create_periodic_mesh().

Note

To go straight from a .msh file, use scifem.periodic.read_periodic_mesh_from_msh(), which reads the pairs out of the model before the reader finalizes it.

scifem.periodic.extract_gmsh_periodic_nodes(model, include_high_order: bool = False, tol: float | None = 1e-08) → PeriodicNodes[source]#

Collect the $Periodic node pairs of model, resolved to roots.

Runs where the gmsh model lives, so serially on the reading rank.

Parameters:
  • model – An initialised gmsh.model carrying a meshed, periodic geometry.

  • include_high_order – Keep the nodes that are not cell vertices. The correspondence this feeds is between vertices, and on a higher-order mesh entities_to_geometry(mesh, 0, vertices) returns each vertex’s corner node, so the default is what is wanted. The resolution below treats the extra nodes no differently, so the flag costs nothing either way.

  • tol – Absolute tolerance for checking each pair against the affine transform gmsh recorded with it; None skips the check. Pairs whose entity stored no transform are never checked. The check earns its place on a file from elsewhere, where nothing has yet compared $Periodic against the coordinates, and is close to tautological for a model this process just built.

Returns:

The pairs, as PeriodicNodes.

Raises:

RuntimeError – If the pairs cycle, disagree on a root, or contradict the affine transform recorded with them.

scifem.periodic.match_vertices_geometric(mesh, indicator, mapping_function, max_chain_length: int | None = None) → VertexCorrespondence[source]#

Pair up the vertices of the seam by evaluating mapping_function on them.

Selects the seam with indicator, moves each selected vertex with mapping_function, and snaps the image onto the nearest vertex of the mesh, checking that it actually landed there.

Parameters:
  • mesh – The mesh to make periodic.

  • indicator – Marks the vertices to be replaced, given coordinates as (3, n).

  • mapping_function – Maps a marked vertex to the one it is identified with, given coordinates as (3, n).

  • max_chain_length – How many times mapping_function may be re-applied to reach a vertex outside indicator, for a mapping that applies one offset per call and so needs several passes to carry a corner to its root. Defaults to mesh.topology.dim, the number of directions such a mesh can be periodic in. Exceeding it raises, which is how a cyclic mapping is caught.

Returns:

The correspondence scifem.periodic rebuilds from.

scifem.periodic.periodic_correspondence_from_nodes(mesh, pairs: PeriodicNodes, root: int = 0) → VertexCorrespondence[source]#

Turn node pairs held on one process into a distributed vertex correspondence.

The topological half of finding the pairs: the identification is given, as indices into the mesh’s input global numbering, and this resolves it against the distribution. No coordinate is read and no tolerance is involved, which is what separates it from scifem.periodic.match_vertices_geometric().

The pairs arrive on one process while the vertices they name are spread over every one, and neither side knows where the other is. A post office resolves that: input global index i is looked after by a fixed rank, scifem.mpi_utils.index_owner, which every process can compute without asking anyone.

  1. every process registers the boundary vertices it holds with the post offices for their indices, saying whether it owns each one;

  2. root sends each pair to the post office of its partner index, which knows who owns that vertex. It tells that owner which pair it answers, and forwards the pair to the post office of the replaced index, which passes it to every process holding a copy;

  3. those processes then ask the partner’s owner directly, which is what tells it who needs the cells at that vertex.

Only the boundary vertices are registered, since an identification pairs nothing else, so the post office stays proportional to the surface rather than the volume. Nothing is gathered: no process holds more than its own block of indices, except root, which holds the pairs it was given.

The rank named for a partner is its vertex owner, which is unique – keeping the join single-valued – and always owns a cell incident to the vertex, which is what src_owner of scifem.periodic.VertexCorrespondence requires.

Collective.

Parameters:
  • mesh – The mesh the pairs refer to, so that mesh.geometry.input_global_indices is the numbering pairs is written in.

  • pairs – The node pairs, meaningful on root only.

  • root – The rank holding pairs.

Returns:

The correspondence scifem.periodic rebuilds from.

Raises:

RuntimeError – If a pair names a node that is not a vertex of the mesh.

scifem.periodic.read_periodic_mesh_from_msh(filename, comm, rank: int = 0, gdim: int = 3, partitioner=None, tag_base: int = 1101, **kwargs)[source]#

Read a .msh file and make the mesh periodic from its $Periodic section.

Owns the gmsh session, because the pairs have to be read out of the model before it is finalized and the usual readers finalize it on the way out.

Collective.

Parameters:
  • filename – The .msh file. Read on rank only.

  • comm – The communicator to distribute the mesh over.

  • rank – The rank that reads the file.

  • gdim – Geometric dimension of the mesh.

  • partitioner – Cell partitioner, passed through to model_to_mesh.

  • tag_base – Passed through to create_periodic_mesh_from_igi().

  • kwargs – Further arguments for model_to_mesh, such as ghost_mode where the installed DOLFINx takes it there.

Returns:

(periodic_mesh, replaced_vertices, replacement_map), as scifem.periodic.create_periodic_mesh().

scifem.periodic.resolve_to_roots(replaced, partner)[source]#

Follow every pair to a node that is not itself replaced.

Parameters:
  • replaced – 0-based node tags, with repeats and possibly several partners each.

  • partner – The node paired with each entry of replaced.

Returns:

each distinct replaced once, and the node it ultimately resolves to.

Return type:

(unique_replaced, root)

Raises:

RuntimeError – If the pairs cycle, or if two routes out of one node disagree on where it ends up.

scifem.periodic.transfer_function_to_parent_mesh(u: Function, parent_mesh: Mesh) → Function[source]#

Transfer a function from a periodic mesh to the parent_mesh it was created from.

Use this to visualize a solution. dolfinx.io.VTXWriter and dolfinx.io.VTKFile.write_function() place one output point per degree of freedom, which a periodic mesh cannot supply a coordinate for: a degree of freedom on the seam belongs to cells on opposite sides of the domain. On the parent mesh the two sides are distinct nodes again, so only the seam is duplicated.

Parameters:
Returns:

A function on parent_mesh, in the same space as u

Note

The transfer is cell-wise, as create_periodic_mesh leaves the cells and the geometry untouched. Merging the seam can change a cell’s orientation, so the two meshes need not agree on the dof transformations of a cell; the cell-wise dolfinx.fem.Function.interpolate() used here accounts for that, so elements that apply their transformations at assembly time, such as RT, N1curl and BDM, transfer correctly too. The writers still require Lagrange or discontinuous Lagrange, so interpolate before writing.

Raises:

ValueError – If parent_mesh does not have the same cells as the mesh of u

scifem.periodic.transfer_meshtags_to_periodic_mesh(mesh: Mesh, periodic_mesh: Mesh, replaced_vertices: ndarray[tuple[Any, ...], dtype[int32]], meshtags: MeshTags) → MeshTags[source]#

Transfer a mesh tag from a mesh to the periodic mesh.

Note

Entities that have been replaced (vertices, edges, faces) are removed from the mesh tag

Parameters:
  • mesh – The original mesh

  • periodic_mesh – The periodic mesh

  • replaced_vertices – The vertices that have been replaced (local to process)

  • meshtags – The mesh tag to transfer

The two MPI tag constants are documented from the module that defines them, since that is where their values are written down.

scifem.periodic.mesh.DEFAULT_TAG_BASE = 1101#

First of the MPI tags the rebuild uses for its consensus exchanges, when the caller names none. The value is arbitrary; what matters is that a caller can move it.

scifem.periodic.mesh.NUM_CONSENSUS_TAGS = 5#

How many consecutive tags from tag_base the rebuild consumes, so that a caller making several overlapping calls knows how far apart to space them.

PETSc utilities#

scifem.petsc.apply_lifting_and_set_bc(b: Vec, a: Iterable[Form] | Iterable[Iterable[Form]], bcs: Iterable[DirichletBC] | Iterable[Iterable[DirichletBC]], x: Vec | None = None, alpha: float = 1.0)[source]#

Apply lifting to a vector and set boundary conditions.

Convenience function to apply lifting and set boundary conditions for multiple matrix types. This modifies the vector b such that

\[\begin{split}b = \begin{pmatrix} b_{free} \\ b_{bc} \end{pmatrix} = \begin{pmatrix} b_{free} - \alpha \sum_{i=0}^n a[i] (u_{bc}[i] - x[i])\\ u_{bc}[0]\\ \vdots\\ u_{bc}[n] \end{pmatrix}.\end{split}\]

where \(b_{free}\) is the free part of the vector, \(b_{bc}\) is the part that has boundary conditions applied, \(u_{bc}[i]\) is the value of the ith boundary condition.

Parameters:
  • b – The vector to apply lifting to.

  • a – Sequence of forms to apply lifting from. If the system is blocked or nested, this is a nested list of forms.

  • bcs – The boundary conditions to apply. If the form is blocked or nested, this is a list, while if it is a single form, this is a nested list.

  • x – Vector to subtract from the boundary conditions. Usually used in a Newton iteration.

  • alpha – The scaling factor for the boundary conditions.

scifem.petsc.ghost_update(x: Vec, insert_mode: InsertMode, scatter_mode: ScatterMode) → None[source]#

Ghost update a vector

scifem.petsc.zero_petsc_vector(b: Vec) → None[source]#

Zero a PETSc vector, including ghosts

XDMF output#

class scifem.xdmf.BaseXDMFFile[source]#
close() → None[source]#

Close the XDMF file.

abstract property data_arrays: list[ndarray[tuple[Any, ...], dtype[floating]]]#

The data arrays.

abstract property data_names: list[str]#

The names of the data.

write(time: float) → None[source]#

Write the point cloud at a given time.

Parameters:

time – The time value.

class scifem.xdmf.FunctionSpaceData(points: ForwardRef('npt.NDArray[np.float64]'), bs: ForwardRef('int'), num_dofs_global: ForwardRef('int'), num_dofs_local: ForwardRef('int'), local_range: ForwardRef('npt.NDArray[np.int64]'), comm: ForwardRef('MPI.Intracomm'))[source]#

Data class for function space information.

bs: int#

Alias for field number 1

comm: Intracomm#

Alias for field number 5

local_range: ndarray[tuple[Any, ...], dtype[int64]]#

Alias for field number 4

num_dofs_global: int#

Alias for field number 2

num_dofs_local: int#

Alias for field number 3

points: ndarray[tuple[Any, ...], dtype[float64]]#

Alias for field number 0

property points_out: ndarray[tuple[Any, ...], dtype[float64]]#

Pad points to be 3D

class scifem.xdmf.NumpyXDMFFile(filename: PathLike, arrays: list[ndarray[tuple[Any, ...], dtype[floating]]], function_space_data: FunctionSpaceData, filemode: Literal['r', 'a', 'w'] = 'w', backend: Literal['h5py', 'adios2'] = 'adios2', array_names: list[str] | None = None)[source]#
property data_arrays: list[ndarray[tuple[Any, ...], dtype[floating]]]#

The data arrays.

property data_names: list[str]#

The names of the data.

class scifem.xdmf.XDMFData(*args, **kwargs)[source]#
class scifem.xdmf.XDMFFile(filename: PathLike, functions: Sequence[Function], filemode: Literal['r', 'a', 'w'] = 'w', backend: Literal['h5py', 'adios2'] = 'adios2')[source]#
property data_arrays: list[ndarray[tuple[Any, ...], dtype[floating]]]#

The data arrays.

property data_names: list[str]#

The names of the data.

scifem.xdmf.check_function_space(functions: Sequence[Function]) → FunctionSpaceData | None[source]#

Check that all functions are in the same function space, and return the function space data.

Parameters:

functions – The functions to check.

Returns:

The function space data if all functions are in the same function space, otherwise None.

scifem.xdmf.create_function_space_data(V: FunctionSpace) → FunctionSpaceData[source]#

Create function space data from a function space.

Parameters:

V – The function space.

Returns:

The function space data.

scifem.xdmf.create_pointcloud(filename: PathLike, functions: Sequence[Function]) → None[source]#

Create a point cloud from a list of functions to be visualized in Paraview. The point cloud is written to a file in XDMF format.

Parameters:
  • filename – The name of the file to write the point cloud to.

  • functions – The functions to write to the point cloud.

Note

This is useful for visualizing functions in quadrature spaces.

Note

Any function space that can call tabulate_dof_coordinates can be used.

Note

ADIOS2 is the preferred backend for writing HDF5 files, and will be used if available. If ADIOS2 is not available, h5py will be used.

scifem.xdmf.h5pyfile(h5name, filemode='r', force_serial: bool = False, comm=None)[source]#

Context manager for opening an HDF5 file with h5py.

Parameters:
  • h5name – The name of the HDF5 file.

  • filemode – The file mode.

  • force_serial – Force serial access to the file.

  • comm – The MPI communicator

scifem.xdmf.write_hdf5_adios(functions: Sequence[Function], h5name: Path, data: FunctionSpaceData | None = None) → None[source]#

Write the point cloud to an HDF5 file using ADIOS2.

Parameters:
  • functions – The functions to write to the point cloud.

  • h5name – The name of the file to write the point cloud to.

  • data – The function space data.

Note

All input functions has to share the same function space.

scifem.xdmf.write_hdf5_h5py(functions: Sequence[Function], h5name: Path, data: FunctionSpaceData | None = None) → None[source]#

Write the point cloud to an HDF5 file using h5py.

Parameters:
  • functions – The functions to write to the point cloud.

  • h5name – The name of the file to write the point cloud to.

  • data – The function space data.

Note

All input functions has to share the same function space.

scifem.xdmf.write_xdmf(functions: Sequence[Function], filename: PathLike, h5name: Path, xdmfdata: XDMFData) → None[source]#

Write the XDMF file for the point cloud.

Parameters:
  • functions – List of functions to write to the point cloud.

  • filename – The name of the file to write the XDMF to.

  • h5name – The name of the HDF5 file.

  • xdmfdata – The XDMF data.

Note

This function does not check the validity of the input data, i.e. all functions in functions are assumed to share the same function space, including block_size.

Compat functions#

Layer for small backward compatibility wrappers for DOLFINx

scifem.compat.compute_integration_domains(integral_type: IntegralType, topology: Topology, entities: ndarray[tuple[Any, ...], dtype[int32]]) → ndarray[tuple[Any, ...], dtype[int32]][source]#

dolfinx.fem.compute_integration_domains() across DOLFINx versions.

Older versions also take the dimension of entities, which follows from integral_type: the cells for a cell integral, the facets otherwise.

Parameters:
  • integral_type – The type of integral the entities are for.

  • topology – The topology of the mesh the entities belong to.

  • entities – The entities, local to the process.

Returns:

The integration entities, flattened, as returned by DOLFINx.

scifem.compat.create_cell_permutations(topology: Topology)[source]#

Compute the packed per-cell permutation info, across DOLFINx versions.

Parameters:

topology – The topology to compute the permutation info of.

scifem.compat.create_cpp_finite_element(constructor, element, gdim: int, block_shape: tuple[int, ...])[source]#

Create a blocked C++ finite element across the supported DOLFINx versions.

FEniCS/dolfinx#4511 added the geometric dimension to the constructor. Versions from before block_shape existed raise TypeError, for the caller to fall back on block_size.

Parameters:
  • constructor – dolfinx.cpp.fem.FiniteElement_float32 or _float64.

  • element – The C++ Basix element.

  • gdim – Geometric dimension of the mesh.

  • block_shape – Block shape of the element.

scifem.compat.create_partitioner(ghost_mode: GhostMode = GhostMode.shared_facet, max_facet_to_cell_links: int = 2)[source]#

Create a partitioner across the supported DOLFINx versions.

The constructor has changed shape more than once and none of the forms is introspectable, so they are told apart by the TypeError the call itself raises:

  1. communicator and ghost mode only, everything else through setters;

  2. the same with max_facet_to_cell_links passed positionally;

  3. as (2) without the communicator.

Parameters:
  • comm – The communicator of the new partitioner.

  • ghost_mode – The ghosting mode to use.

  • max_facet_to_cell_links – Maximum number of cells that can share a facet. Used by forms (2) and (3) only, which take it directly; form (1) does not accept it.

Returns:

The default partitioner

scifem.compat.get_facet_permutations(topology: Topology) → ndarray[tuple[Any, ...], dtype[uint8]][source]#

The permutation of every facet of every cell, across DOLFINx versions.

Each value encodes how a facet is oriented as seen from a cell, relative to a low-to-high ordering of its global vertex indices, as FFCx uses for quadrature_permutation.

Parameters:

topology – The topology to compute the facet permutations of.

Returns:

The permutations, shape (num_cells, num_facets_per_cell), ghost cells included.

Layer for small backward compatibility wrappers for UFL

scifem.ufl_compat.apply_pullback_inverse(pullback: AbstractPullback, expr: Expr, domain: Mesh) → Expr[source]#

Map expr from the physical cell to the reference cell, across UFL versions.

Uses pullback.apply_inverse, added in FEniCS/ufl#511, when the installed UFL has it, and otherwise a closed form for the identity, contravariant Piola and covariant Piola pullbacks.

Parameters:
  • pullback – The element’s pullback.

  • expr – A physical-cell expression.

  • domain – The domain whose Jacobian relates the two cells.

Returns:

expr pulled back to the reference cell.

Raises:

NotImplementedError – If this UFL has no apply_inverse and there is no closed form for pullback here.