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.Meshobjectdim – Dimension of the entities to mark
entities_list –
A list of tuples with the following elements:
index 0: The tag to assign to the entitiesindex 1: A function that takes a point and returns a boolean array indicating whether the point is inside the entityindex 2: Optional, if True, the entities will be marked on the boundary
- Returns:
A
dolfinx.mesh.MeshTagsobject 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_tagfrom a parent mesh to asubmesh.- 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_tagis the tag on the submesh andsub_to_parent_entity_mapis 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
ctthat is marked with any of themarkers. 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).
Interpolation#
- scifem.interpolation_matrix(expr: Expr, Q: FunctionSpace) MatrixCSR[source]#
Create the interpolation matrix \(\Lambda\) of a
UFL-expressionsuch 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-expressionsuch 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.Expressionat 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_surfaceon a facet submesh by zero into the functionu_volume.The reverse of
interpolate_to_surface_submesh():u_volumetakes the values ofu_surfaceat its nodes on the facets and is zero at every other node. SeeSurfaceSubmeshExtension, 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 (seeapply()).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_volumefor 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 ofsubmesh_facets, as forinterpolate_to_surface_submesh().entity_maps – The submesh’s entity map, needed for a volume space that is not Lagrange.
- Raises:
ValueError – If
V_volumeis discontinuous, the value sizes do not fit, orentity_mapsis missing where it is needed.NotImplementedError – If
V_surfaceneeds 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_volumeto the extension ofu_surface.Unlike
interpolate_to_surface_submesh()andscifem.interpolate_function_onto_facet_dofs(), this does not calldolfinx.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 byweights, 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.
- property basis: ndarray[tuple[Any, ...], dtype[floating]]#
The map from each facet’s
surface_dofsto itsvolume_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).
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
ireorders data given at the closure dofs of the local entityentities[i, 1]of cellentities[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’spermute_subentity_closure_invwith 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 fordim = 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
sizeis 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 ofV(a node, for a blocked space) at thej-th closure dof of the local entityentities[i, 1]of cellentities[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. Seecompute_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.Exprwithin 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]#
-
- 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.
- 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$Periodicsection, 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, andperiodic_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
meshare identified with which, and which ranks hold each end.This is what
scifem.periodicrebuilds 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.PointOwnershipDatafor 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$Periodicsection of a gmsh file, builds aVertexCorrespondenceand 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_TAGSconsecutive 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)wherenew_meshis the new mesh with periodicity,replaced_verticesis a list of vertices of the input mesh that has been replaced (local to process).replacement_mapis 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
VertexCorrespondenceis 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 thematch_vertices_geometric()side of it, and_build_periodic_meshsees 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_indicesis 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 whencreate_meshdrops nodes no cell references.root – The rank holding the pairs.
tag_base – The first of
scifem.periodic.mesh.NUM_CONSENSUS_TAGSconsecutive 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:
Note
To go straight from a
.mshfile, usescifem.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
$Periodicnode pairs of model, resolved to roots.Runs where the gmsh model lives, so serially on the reading rank.
- Parameters:
model – An initialised
gmsh.modelcarrying 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.periodicrebuilds 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
iis looked after by a fixed rank,scifem.mpi_utils.index_owner, which every process can compute without asking anyone.every process registers the boundary vertices it holds with the post offices for their indices, saying whether it owns each one;
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;
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_ownerofscifem.periodic.VertexCorrespondencerequires.Collective.
- Parameters:
mesh – The mesh the pairs refer to, so that
mesh.geometry.input_global_indicesis the numbering pairs is written in.pairs – The node pairs, meaningful on root only.
root – The rank holding pairs.
- Returns:
The correspondence
scifem.periodicrebuilds 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
.mshfile and make the mesh periodic from its$Periodicsection.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
.mshfile. 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 asghost_modewhere the installed DOLFINx takes it there.
- Returns:
(periodic_mesh, replaced_vertices, replacement_map), asscifem.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_meshit was created from.Use this to visualize a solution.
dolfinx.io.VTXWriteranddolfinx.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:
u – The function on the periodic mesh
parent_mesh – The mesh that was passed to
scifem.periodic.create_periodic_mesh()
- Returns:
A function on
parent_mesh, in the same space asu
Note
The transfer is cell-wise, as
create_periodic_meshleaves 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-wisedolfinx.fem.Function.interpolate()used here accounts for that, so elements that apply their transformations at assembly time, such asRT,N1curlandBDM, transfer correctly too. The writers still require Lagrange or discontinuous Lagrange, so interpolate before writing.- Raises:
ValueError – If
parent_meshdoes not have the same cells as the mesh ofu
- 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_basethe 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
XDMF output#
- class scifem.xdmf.BaseXDMFFile[source]#
-
- abstract property data_names: list[str]#
The names of the data.
- 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
- num_dofs_global: int#
Alias for field number 2
- num_dofs_local: int#
Alias for field number 3
- 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_names: list[str]#
The names of the data.
- class scifem.xdmf.XDMFFile(filename: PathLike, functions: Sequence[Function], filemode: Literal['r', 'a', 'w'] = 'w', backend: Literal['h5py', 'adios2'] = 'adios2')[source]#
-
- 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
functionshas 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
functionshas 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
functionsare assumed to share the same function space, includingblock_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 fromintegral_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_shapeexisted raiseTypeError, for the caller to fall back onblock_size.- Parameters:
constructor –
dolfinx.cpp.fem.FiniteElement_float32or_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:
communicator and ghost mode only, everything else through setters;
the same with max_facet_to_cell_links passed positionally;
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
exprfrom 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:
exprpulled back to the reference cell.- Raises:
NotImplementedError – If this UFL has no
apply_inverseand there is no closed form forpullbackhere.