Discretizations

Abstract interface

Error handling

exception meshmode.discretization.ElementGroupTypeError

A TypeError specific for handling element groups. This exception may be raised to indicate whenever an improper operation or function is applied to a particular subclass of ElementGroupBase.

exception meshmode.discretization.NoninterpolatoryElementGroupError

A specialized ElementGroupTypeError that may be raised whenever non-interpolatory element groups are being used for interpolation.

Base classes

class meshmode.discretization.ElementGroupBase(mesh_el_group: MeshElementGroup, order: int, index: int | None = None)

Defines a discrete function space on a homogeneous (in terms of element type and order) subset of a Discretization. These correspond one-to-one with meshmode.mesh.MeshElementGroup. Responsible for all bulk data handling in Discretization.

mesh_el_group

order

is_affine

A bool flag that is True if the local-to-global parametrization of all the elements in the group is affine. Based on meshmode.mesh.MeshElementGroup.is_affine.

nelements

The total number of polygonal elements in the meshmode.mesh.MeshElementGroup.

nunit_dofs

The number of degrees of freedom (“DOFs”) associated with a single element.

ndofs

The total number of degrees of freedom (“DOFs”) associated with the entire element group.

dim

The number of spatial dimensions in which the functions in space operate.

shape

Returns a subclass of modepy.Shape representing the reference element defining the element group.

space

Returns a modepy.FunctionSpace representing the underlying polynomial space defined on the element group’s reference element.

__init__(mesh_el_group: MeshElementGroup, order: int, index: int | None = None) → None

discretization_key() → Hashable

Return a hashable, equality-comparable object that fully describes the per-element discretization used by this element group. (This should cover all parts of the Ciarlet Triple: reference element, shape functions, and the linear functionals defining the degrees of freedom.) The object should be independent, however, of the (global) elements that make up the group.

The structure of the element is not specified, but it must be globally unique to this element group.

class meshmode.discretization.ElementGroupFactory(*args, **kwargs)

A typing.Protocol specifying the interface for group factories.

__call__(mesh_el_group: MeshElementGroup, index: int | None = None) → ElementGroupBase

Create a new ElementGroupBase for the given mesh_el_group.

class meshmode.discretization.NodalElementGroupBase(mesh_el_group: MeshElementGroup, order: int, index: int | None = None)

Base class for nodal element groups, defined as finite elements equipped with nodes. Nodes are specific locations defined on the reference element (shape) defining a degree of freedom by point evaluation at that location. Such element groups have an associated quadrature rule to perform numerical integration, but are not necessarily usable (unisolvent) for interpolation.

Inherits from ElementGroupBase.

unit_nodes

Returns a numpy.ndarray of shape (dim, nunit_dofs) of reference coordinates of interpolation nodes.

Note: this method dispatches to the nodes of the underlying quadrature rule. This means, for interpolatory element groups, interpolation nodes are collocated with quadrature nodes.

abstract quadrature_rule()

Returns a modepy.Quadrature object for the element group.

class meshmode.discretization.ElementGroupWithBasis(mesh_el_group: MeshElementGroup, order: int, index: int | None = None)

Base class for element groups which possess an explicit basis for the underlying function space space.

Inherits from ElementGroupBase.

abstract basis_obj()

Returns the modepy.Basis which spans the underlying space.

is_orthonormal_basis()

Returns a bool flag that is True if the basis corresponding to the element group is orthonormal with respect to the \(L^2\) inner-product.

class meshmode.discretization.InterpolatoryElementGroupBase(mesh_el_group: MeshElementGroup, order: int, index: int | None = None)

An element group equipped with both an explicit basis for the underlying space, and a set of nodal locations on the shape. These element groups are unisolvent in the Ciarlet sense, meaning the dimension of space matches the number of interpolatory nodal locations. These element groups are therefore suitable for interpolation and differentiation.

Inherits from NodalElementGroupBase and ElementGroupWithBasis.

abstract mass_matrix()

Return a numpy.ndarray of shape (nunit_nodes, nunit_nodes), which is defined as the operator \(M\), with

\[M_{ij} = \int_{K} \phi_i \cdot \phi_j \mathrm{d}x,\]

where \(K\) denotes a cell and \(\phi_i\) is the basis spanning the underlying space.

abstract diff_matrices()

Return a dim-long tuple of numpy.ndarray of shape (nunit_nodes, nunit_nodes), each of which, when applied to an array of nodal values, take derivatives in the reference \((r, s, t)\) directions.

class meshmode.discretization.ModalElementGroupBase(mesh_el_group: MeshElementGroup, order: int, index: int | None = None)

An element group equipped with a function space and a hierarchical basis that is orthonormal with respect to the \(L^2\) inner product.

Inherits from ElementGroupWithBasis.

Discretization class

meshmode.discretization.num_reference_derivative(discr: Discretization, ref_axes: Iterable[int], vec: DOFArray) → DOFArray

Parameters:

ref_axes – an Iterable of indices that define the sequence of derivatives to vec. For example, (0, 1, 1) would take a third partial derivative, one in the first axis and two in the second axis.

class meshmode.discretization.Discretization(actx: ArrayContext, mesh: Mesh, group_factory: ElementGroupFactory, real_dtype: dtype | None = None, _force_actx_clone: bool = True)

An unstructured composite discretization.

real_dtype

complex_dtype

mesh

dim

ambient_dim

ndofs

groups

is_nodal

A bool indicating whether the Discretization is defined over element groups subclasses of NodalElementGroupBase.

is_modal

A bool indicating whether the Discretization is defined over element groups subclasses of ModalElementGroupBase.

__init__(actx: ArrayContext, mesh: Mesh, group_factory: ElementGroupFactory, real_dtype: dtype | None = None, _force_actx_clone: bool = True) → None

Parameters:

• actx – an arraycontext.ArrayContext used to perform computation needed during initial set-up of the discretization.

• mesh – a Mesh over which the discretization is built.

• group_factory – an ElementGroupFactory.

• real_dtype – The numpy data type used for representing real data, either numpy.float32 or numpy.float64.

copy(actx: ArrayContext | None = None, mesh: Mesh | None = None, group_factory: ElementGroupFactory | None = None, real_dtype: dtype | None = None) → Discretization

Creates a new object of the same type with all arguments that are not None replaced. The copy is not recursive (e.g. it does not call meshmode.mesh.Mesh.copy()).

empty(actx: ArrayContext, dtype: dtype | None = None) → DOFArray

Return an empty DOFArray.

Parameters:

dtype – type special value ‘c’ will result in a vector of dtype complex_dtype. If None (the default), a real vector will be returned.

zeros(actx: ArrayContext, dtype: dtype | None = None) → DOFArray

Return a zero-initialized DOFArray.

Parameters:

dtype – type special value ‘c’ will result in a vector of dtype complex_dtype. If None (the default), a real vector will be returned.

empty_like(array: DOFArray) → DOFArray

zeros_like(array: DOFArray) → DOFArray

nodes(cached: bool = True) → ndarray

Parameters:

cached – A bool indicating whether the computed nodes should be stored for future use.

Returns:

object array of shape (ambient_dim,) containing DOFArrays of (global) nodal locations on the mesh.

num_reference_derivative(ref_axes: Iterable[int], vec: DOFArray) → DOFArray

quad_weights() → DOFArray

Returns:

a DOFArray with quadrature weights.

Element Groups for Composite Polynomial Discretization

Group types

meshmode.discretization.poly_element.mass_matrix(grp: InterpolatoryElementGroupBase) → ndarray

meshmode.discretization.poly_element.diff_matrices(grp: InterpolatoryElementGroupBase) → Tuple[ndarray]

Simplicial group types

class meshmode.discretization.poly_element.ModalSimplexElementGroup(mesh_el_group: MeshElementGroup, order: int, index: int | None = None)

class meshmode.discretization.poly_element.InterpolatoryQuadratureSimplexElementGroup(mesh_el_group: MeshElementGroup, order: int, index: int | None = None)

Elemental discretization supplying a high-order quadrature rule with a number of nodes matching the number of polynomials in \(P^k\), hence usable for differentiation and interpolation.

No interpolation nodes are present on the boundary of the simplex.

class meshmode.discretization.poly_element.QuadratureSimplexElementGroup(mesh_el_group: MeshElementGroup, order: int, index: int | None = None)

Elemental discretization supplying a high-order quadrature rule with a number of nodes which does not necessarily match the number of polynomials in \(P^k\). This discretization therefore excels at quadarature, but is not necessarily usable for interpolation.

No interpolation nodes are present on the boundary of the simplex.

class meshmode.discretization.poly_element.PolynomialWarpAndBlend2DRestrictingElementGroup(mesh_el_group: MeshElementGroup, order: int, index: int | None = None)

Elemental discretization with a number of nodes matching the number of polynomials in \(P^k\), hence usable for differentiation and interpolation. Interpolation nodes edge-clustered for avoidance of Runge phenomena. Nodes are present on the boundary of the simplex. Provides nodes in two and fewer dimensions, based on the 2D warp-and-blend nodes and their facial restrictions.

Uses modepy.warp_and_blend_nodes().

class meshmode.discretization.poly_element.PolynomialWarpAndBlend3DRestrictingElementGroup(mesh_el_group: MeshElementGroup, order: int, index: int | None = None)

Elemental discretization with a number of nodes matching the number of polynomials in \(P^k\), hence usable for differentiation and interpolation. Interpolation nodes edge-clustered for avoidance of Runge phenomena. Nodes are present on the boundary of the simplex. Provides nodes in two and fewer dimensions, based on the 3D warp-and-blend nodes and their facial restrictions.

Uses modepy.warp_and_blend_nodes().

class meshmode.discretization.poly_element.PolynomialRecursiveNodesElementGroup(mesh_el_group, order, family, index=None)

Elemental discretization with a number of nodes matching the number of polynomials in \(P^k\), hence usable for differentiation and interpolation. Interpolation nodes edge-clustered for avoidance of Runge phenomena. Depending on the family argument, nodes may be present on the boundary of the simplex. See [Isaac20] for details.

Supports a choice of the base family of 1D nodes, see the documentation of the family argument to recursivenodes.recursive_nodes().

Requires recursivenodes to be installed.

[Isaac20]

Tobin Isaac. Recursive, parameter-free, explicitly defined interpolation nodes for simplices. Arxiv preprint.

New in version 2020.2.

class meshmode.discretization.poly_element.PolynomialEquidistantSimplexElementGroup(mesh_el_group: MeshElementGroup, order: int, index: int | None = None)

Elemental discretization with a number of nodes matching the number of polynomials in \(P^k\), hence usable for differentiation and interpolation. Interpolation nodes are present on the boundary of the simplex.

New in version 2016.1.

class meshmode.discretization.poly_element.PolynomialGivenNodesElementGroup(mesh_el_group, order, unit_nodes, index=None)

Elemental discretization with a number of nodes matching the number of polynomials in \(P^k\), hence usable for differentiation and interpolation. Uses nodes given by the user.

Tensor product group types

class meshmode.discretization.poly_element.ModalTensorProductElementGroup(mesh_el_group: MeshElementGroup, order: int, index: int | None = None)

class meshmode.discretization.poly_element.GaussLegendreTensorProductElementGroup(mesh_el_group, order, index=None)

Elemental discretization supplying a high-order quadrature rule with a number of nodes matching the number of polynomials in the tensor product basis, hence usable for differentiation and interpolation.

No interpolation nodes are present on the boundary of the hypercube.

class meshmode.discretization.poly_element.LegendreGaussLobattoTensorProductElementGroup(mesh_el_group, order, index=None)

Elemental discretization supplying a high-order quadrature rule with a number of nodes matching the number of polynomials in the tensor product basis, hence usable for differentiation and interpolation. Nodes sufficient for unisolvency are present on the boundary of the hypercube.

Uses legendre_gauss_lobatto_nodes().

class meshmode.discretization.poly_element.EquidistantTensorProductElementGroup(mesh_el_group, order, index=None)

Elemental discretization supplying a high-order quadrature rule with a number of nodes matching the number of polynomials in the tensor product basis, hence usable for differentiation and interpolation. Nodes sufficient for unisolvency are present on the boundary of the hypercube.

Uses equidistant_nodes().

Group factories

class meshmode.discretization.poly_element.HomogeneousOrderBasedGroupFactory(order: int)

Element group factory for a single type of meshmode.mesh.MeshElementGroup and fixed order.

mesh_group_class

group_class

order

__init__(order: int) → None

Parameters:

order – integer denoting the order of the ElementGroupBase. The exact interpretation of the order is left to each individual class, as given by group_class.

__call__(mesh_el_group, index=None)

Returns:

an element group of type group_class and order order.

class meshmode.discretization.poly_element.TypeMappingGroupFactory(order, mesh_group_class_to_factory)

Element group factory that supports multiple types of MeshElementGroups, defined through the mapping mesh_group_class_to_factory.

order

mesh_group_class_to_factory

A dict from MeshElementGroups to factory callables that return a corresponding ElementGroupBase of order order.

__init__(order, mesh_group_class_to_factory)

Parameters:

mesh_group_class_to_factory – a dict from MeshElementGroup subclasses to ElementGroupBase subclasses or ElementGroupFactory instances.

__call__(mesh_el_group, index=None)

Create a new ElementGroupBase for the given mesh_el_group.

Simplicial group factories

class meshmode.discretization.poly_element.ModalSimplexGroupFactory(order: int)

class meshmode.discretization.poly_element.InterpolatoryQuadratureSimplexGroupFactory(order: int)

class meshmode.discretization.poly_element.QuadratureSimplexGroupFactory(order: int)

class meshmode.discretization.poly_element.PolynomialWarpAndBlend2DRestrictingGroupFactory(order: int)

class meshmode.discretization.poly_element.PolynomialWarpAndBlend3DRestrictingGroupFactory(order: int)

class meshmode.discretization.poly_element.PolynomialRecursiveNodesGroupFactory(order, family)

class meshmode.discretization.poly_element.PolynomialEquidistantSimplexGroupFactory(order: int)

New in version 2016.1.

class meshmode.discretization.poly_element.PolynomialGivenNodesGroupFactory(order, unit_nodes)

Tensor product group factories

class meshmode.discretization.poly_element.ModalTensorProductGroupFactory(order: int)

class meshmode.discretization.poly_element.GaussLegendreTensorProductGroupFactory(order: int)

class meshmode.discretization.poly_element.LegendreGaussLobattoTensorProductGroupFactory(order: int)

class meshmode.discretization.poly_element.EquidistantTensorProductGroupFactory(order: int)

Type-based group factories

class meshmode.discretization.poly_element.InterpolatoryEdgeClusteredGroupFactory(order)

Element group factory for all supported MeshElementGroups that constructs (recommended) element groups with edge-clustered nodes that can be used for interpolation.

class meshmode.discretization.poly_element.InterpolatoryQuadratureGroupFactory(order)

Element group factory for all supported MeshElementGroups that constructs (recommended) element groups with nodes that can be used for interpolation and high-order quadrature.

class meshmode.discretization.poly_element.InterpolatoryEquidistantGroupFactory(order)

Element group factory for all supported MeshElementGroups that constructs (recommended) element groups with equidistant nodes that can be used for interpolation.

class meshmode.discretization.poly_element.QuadratureGroupFactory(order)

Element group factory for all supported MeshElementGroups that constructs (recommended) element groups with nodes that can be used for high-order quadrature, but (not necessarily) for interpolation.

class meshmode.discretization.poly_element.ModalGroupFactory(order)

Element group factory for all supported MeshElementGroups that constructs (recommended) element groups with modal degrees of freedom.

Visualization

meshmode.discretization.visualization.make_visualizer(actx, discr, vis_order=None, element_shrink_factor=None, force_equidistant=False)

Parameters:

• vis_order – order of the visualization DOFs.

• element_shrink_factor – number in \((0, 1]\).

• force_equidistant – if True, the visualization is done on equidistant nodes. If plotting high-order Lagrange VTK elements, this needs to be set to True.

class meshmode.discretization.visualization.Visualizer(connection, element_shrink_factor=None, is_equidistant=False, _vtk_linear_connectivity=None, _vtk_lagrange_connectivity=None)

show_scalar_in_mayavi(field, **kwargs)

show_scalar_in_matplotlib_3d(field, **kwargs)

write_vtk_file(file_name, names_and_fields, compressor=None, real_only=False, overwrite=False, use_high_order=None, par_manifest_filename=None, par_file_names=None)

Write a Vtk XML file (typical extension .vtu) containing the visualization data in names_and_fields. Can optionally also write manifests for distributed memory simulation (typical extension .pvtu). See also write_parallel_vtk_file() for a convenience wrapper.

Parameters:

• names_and_fields – A list of tuples (name, value), where name is a string and value is a DOFArray or a constant, or an object array of those. value may also be a data class (see dataclasses), whose attributes will be inserted into the visualization with their names prefixed by name. If value is None, then there is no data to write and the corresponding name will not appear in the data file. If value is None, it should be None collectively across all ranks for parallel writes; otherwise the behavior of this routine is undefined.

• overwrite – If True, silently overwrite existing files.

• use_high_order – Writes arbitrary order Lagrange VTK elements. These elements are described in this blog post and are available in VTK 8.1 and newer.

• par_manifest_filename – If not None write a distributed-memory manifest with this file name if file_name matches the first entry in par_file_names.

• par_file_names – A list of file names of visualization files to include in the distributed-memory manifest.

Changed in version 2020.2:

• Added par_manifest_filename and par_file_names.

• Added use_high_order.

write_parallel_vtk_file(mpi_comm, file_name_pattern, names_and_fields, compressor=None, real_only=False, overwrite=False, use_high_order=None, par_manifest_filename=None)

A convenience wrapper around write_vtk_file() for distributed-memory visualization.

Parameters:

• mpi_comm – An object that supports mpi_comm.Get_rank() and mpi_comm.Get_size() method calls, typically (but not necessarily) an instance of mpi4py.Comm. This is used to determine the current rank as well as the total number of files being written. May also be None in which case a unit-size communicator is assumed.

• file_name_pattern – A file name pattern (required to end in .vtu) that will be used with str.format() with an (integer) argument of rank to obtain the per-rank file name. Relative path names are also supported.

• par_manifest_filename – as in write_vtk_file(). If not given, par_manifest_filename is synthesized by substituting rank 0 into file_name_pattern and replacing the file extension with .pvtu.

See write_vtk_file() for the meaning of the remainder of the arguments.

New in version 2020.2.

write_vtkhdf_file(file_name: str, names_and_fields: List[Tuple[str, Any]], *, comm=None, use_high_order: bool = False, real_only: bool = False, overwrite: bool = False, h5_file_options: Dict[str, Any] | None = None, dset_options: Dict[str, Any] | None = None) → None

Write a VTK HDF5 file (typical extension '.hdf') containing the visualization fields in names_and_fields.

This function requires h5py and has support for parallel writes through mpi4py.

Parameters:

• comm – an mpi4py.Comm-like interface that supports Get_rank, Get_size, scan and reduce. The last two are required to gather global information about the points and cells in the discretizations.

• h5_file_options – a dict passed directly to h5py.File that allows controlling chunking, compatibility, etc.

• dataset_options – a dict passed directly to h5py.Group.create_dataset().

write_xdmf_file(file_name, names_and_fields, attrs=None, h5_file_options=None, dataset_options=None, real_only=False, overwrite=False)

Write an XDMF file (with an .xmf extension) containing the arrays in names_and_fields. The heavy data is written to binary HDF5 files, which requires installing h5py. Distributed memory visualization is not yet supported.

Parameters:

• names_and_fields – a list of (name, array), where array is an array-like object (see Visualizer.write_vtk_file()).

• attrs – a dict of scalar attributes that will be saved in the root HDF5 group.

• h5_file_options – a dict passed directly to h5py.File that allows controlling chunking, compatibility, etc.

• dataset_options – a dict passed directly to h5py.Group.create_dataset().

copy_with_same_connectivity(actx, discr, skip_tests=False)

Makes a copy of the visualizer for a Discretization with the same group structure as the original discretization. This can be useful when the geometry is mapped (e.g. using affine_map()) and the connectivity can be reused.

The “same group structure” here means that the two discretizations should have the same group types, number of elements, degrees of freedom, etc.

Parameters:

skip_tests – If True, no checks in the group structure of the discretizations are performed.

meshmode.discretization.visualization.write_nodal_adjacency_vtk_file(file_name, mesh, compressor=None, overwrite=False)