Modernize 2024

.ci/install-for-firedrake.sh

@@ -28,7 +28,4 @@ python -m pip install --force-reinstall git+https://github.com/inducer/pytools.g
 pip uninstall -y loopy
 pip install "git+https://github.com/inducer/loopy.git#egg=loopy"
 
-# https://github.com/OP2/PyOP2/pull/627
-(cd /home/firedrake/firedrake/src/PyOP2 && git pull https://github.com/OP2/PyOP2.git e56d26f219e962cf9423fc84406a8a0656eb364f)
-
 pip install .

examples/parallel-vtkhdf.py

@@ -69,7 +69,7 @@ def main(*, ambient_dim: int) -> None:
 
     vector_field = actx.thaw(discr.nodes())
     scalar_field = actx.np.sin(vector_field[0])
-    part_id = 1.0 + comm.rank + discr.zeros(actx)     # type: ignore[operator]
+    part_id = 1.0 + comm.rank + discr.zeros(actx)
     logger.info("[%4d] fields: finished", comm.rank)
 
     from meshmode.discretization.visualization import make_visualizer

meshmode/discretization/connection/face.py

@@ -226,7 +226,7 @@ def make_face_restriction(
 
     # }}}
 
-    from meshmode.mesh import _ModepyElementGroup, make_mesh
+    from meshmode.mesh import ModepyElementGroup, make_mesh
     bdry_mesh_groups = []
     connection_data = {}
 
@@ -235,7 +235,7 @@ def make_face_restriction(
 
         mgrp = grp.mesh_el_group
 
-        if not isinstance(mgrp, _ModepyElementGroup):
+        if not isinstance(mgrp, ModepyElementGroup):
             raise NotImplementedError("can only take boundary of "
                     "meshes based on SimplexElementGroup and "
                     "TensorProductElementGroup")
@@ -301,7 +301,7 @@ def make_face_restriction(
                 bdry_unit_nodes = mp.edge_clustered_nodes_for_space(space, face)
 
                 vol_basis = mp.basis_for_space(
-                        mgrp._modepy_space, mgrp._modepy_shape).functions
+                        mgrp.space, mgrp.shape).functions
 
                 vertex_indices = np.empty(
                         (ngroup_bdry_elements, face.nvertices),

meshmode/discretization/poly_element.py

@@ -118,10 +118,10 @@
 def from_mesh_interp_matrix(grp: InterpolatoryElementGroupBase) -> np.ndarray:
     meg = grp.mesh_el_group
 
-    from meshmode.mesh import _ModepyElementGroup
-    assert isinstance(meg, _ModepyElementGroup)
+    from meshmode.mesh import ModepyElementGroup
+    assert isinstance(meg, ModepyElementGroup)
 
-    meg_basis = mp.basis_for_space(meg._modepy_space, meg._modepy_shape)
+    meg_basis = mp.basis_for_space(meg.space, meg.shape)
     return mp.resampling_matrix(
             meg_basis.functions,
             grp.unit_nodes,

meshmode/discretization/visualization.py

@@ -316,12 +316,12 @@ def tensor_cell_types(self):
     def connectivity_for_element_group(self, grp):
         import modepy as mp
 
-        from meshmode.mesh import _ModepyElementGroup
+        from meshmode.mesh import ModepyElementGroup
 
-        if isinstance(grp.mesh_el_group, _ModepyElementGroup):
-            shape = grp.mesh_el_group._modepy_shape
+        if isinstance(grp.mesh_el_group, ModepyElementGroup):
+            shape = grp.mesh_el_group.shape
             space = mp.space_for_shape(shape, grp.order)
-            assert type(space) == type(grp.mesh_el_group._modepy_space)  # noqa: E721
+            assert type(space) == type(grp.mesh_el_group.space)  # noqa: E721
 
             node_tuples = mp.node_tuples_for_space(space)
 

meshmode/dof_array.py

@@ -29,7 +29,7 @@
 from contextlib import contextmanager
 from functools import partial, update_wrapper
 from numbers import Number
-from typing import Any
+from typing import Any, TypeAlias, Union
 from warnings import warn
 
 import numpy as np
@@ -69,6 +69,9 @@
 
 # {{{ DOFArray
 
+ArithType: TypeAlias = Union["DOFArray", int, float, complex, np.generic]
+
+
 @with_container_arithmetic(
         bcast_obj_array=True,
         rel_comparison=True,
@@ -373,6 +376,21 @@ def _deserialize_init_arrays_code(cls, template_instance_name, args):
         # Why tuple([...])? https://stackoverflow.com/a/48592299
         return (f"{template_instance_name}.array_context, tuple([{arg}])")
 
+    # {{{ type stubs for arithmetic, will be replaced by @with_container_arithmetic
+
+    def __neg__(self) -> DOFArray: ...
+    def __abs__(self) -> DOFArray: ...
+    def __add__(self, other: ArithType) -> DOFArray: ...
+    def __radd__(self, other: ArithType) -> DOFArray: ...
+    def __sub__(self, other: ArithType) -> DOFArray: ...
+    def __rsub__(self, other: ArithType) -> DOFArray: ...
+    def __mul__(self, other: ArithType) -> DOFArray: ...
+    def __rmul__(self, other: ArithType) -> DOFArray: ...
+    def __truediv__(self, other: ArithType) -> DOFArray: ...
+    def __rtruediv__(self, other: ArithType) -> DOFArray: ...
+
+    # }}}
+
 # }}}
 
 

meshmode/mesh/__init__.py

@@ -28,6 +28,7 @@
 from abc import ABC, abstractmethod
 from collections.abc import Callable, Collection, Hashable, Iterable, Mapping, Sequence
 from dataclasses import InitVar, dataclass, field, replace
+from functools import partial
 from typing import (
     Any,
     ClassVar,
@@ -41,13 +42,14 @@
 import numpy.linalg as la
 
 import modepy as mp
-from pytools import memoize_method
+from pytools import memoize_method, module_getattr_for_deprecations
 
 from meshmode.mesh.tools import AffineMap, optional_array_equal
 
 
 __doc__ = """
 .. autoclass:: MeshElementGroup
+.. autoclass:: ModepyElementGroup
 .. autoclass:: SimplexElementGroup
 .. autoclass:: TensorProductElementGroup
 
@@ -321,38 +323,45 @@ def make_group(cls, *args: Any, **kwargs: Any) -> MeshElementGroup:
 
 # {{{ modepy-based element group
 
+# https://stackoverflow.com/a/13624858
+class _classproperty(property):  # noqa: N801
+    def __get__(self, owner_self: Any, owner_cls: type | None = None) -> Any:
+        assert self.fget is not None
+        return self.fget(owner_cls)
+
+
 @dataclass(frozen=True, eq=False)
-class _ModepyElementGroup(MeshElementGroup):
+class ModepyElementGroup(MeshElementGroup):
     """
-    .. attribute:: _modepy_shape_cls
+    .. attribute:: modepy_shape_cls
 
         Must be set by subclasses to generate the correct shape and spaces
         attributes for the group.
 
-    .. attribute:: _modepy_shape
-    .. attribute:: _modepy_space
+    .. attribute:: shape
+    .. attribute:: space
     """
 
-    _modepy_shape_cls: ClassVar[type[mp.Shape]]
-    _modepy_shape: mp.Shape = field(repr=False)
-    _modepy_space: mp.FunctionSpace = field(repr=False)
+    shape_cls: ClassVar[type[mp.Shape]]
+    shape: mp.Shape = field(repr=False)
+    space: mp.FunctionSpace = field(repr=False)
 
     @property
     def nvertices(self) -> int:
-        return self._modepy_shape.nvertices     # pylint: disable=no-member
+        return self.shape.nvertices     # pylint: disable=no-member
 
     @property
     @memoize_method
     def _modepy_faces(self) -> Sequence[mp.Face]:
-        return mp.faces_for_shape(self._modepy_shape)
+        return mp.faces_for_shape(self.shape)
 
     @memoize_method
     def face_vertex_indices(self) -> tuple[tuple[int, ...], ...]:
         return tuple(face.volume_vertex_indices for face in self._modepy_faces)
 
     @memoize_method
     def vertex_unit_coordinates(self) -> np.ndarray:
-        return mp.unit_vertices_for_shape(self._modepy_shape).T
+        return mp.unit_vertices_for_shape(self.shape).T
 
     @classmethod
     def make_group(cls,
@@ -361,7 +370,7 @@ def make_group(cls,
                    nodes: np.ndarray,
                    *,
                    unit_nodes: np.ndarray | None = None,
-                   dim: int | None = None) -> _ModepyElementGroup:
+                   dim: int | None = None) -> ModepyElementGroup:
 
         if unit_nodes is None:
             if dim is None:
@@ -374,7 +383,7 @@ def make_group(cls,
                 raise ValueError("'dim' does not match 'unit_nodes' dimension")
 
         # pylint: disable=abstract-class-instantiated
-        shape = cls._modepy_shape_cls(dim)
+        shape = cls.shape_cls(dim)
         space = mp.space_for_shape(shape, order)
 
         if unit_nodes is None:
@@ -388,17 +397,29 @@ def make_group(cls,
                    vertex_indices=vertex_indices,
                    nodes=nodes,
                    unit_nodes=unit_nodes,
-                   _modepy_shape=shape,
-                   _modepy_space=space)
+                   shape=shape,
+                   space=space)
+
+    @_classproperty
+    def _modepy_shape_cls(cls) -> type[mp.Shape]:  # noqa: N805  # pylint: disable=no-self-argument
+        return cls.shape_cls
+
+    @property
+    def _modepy_shape(self) -> mp.Shape:
+        return self.shape
+
+    @property
+    def _modepy_space(self) -> mp.FunctionSpace:
+        return self.space
 
 # }}}
 
 
 @dataclass(frozen=True, eq=False)
-class SimplexElementGroup(_ModepyElementGroup):
+class SimplexElementGroup(ModepyElementGroup):
     r"""Inherits from :class:`MeshElementGroup`."""
 
-    _modepy_shape_cls: ClassVar[type[mp.Shape]] = mp.Simplex
+    shape_cls: ClassVar[type[mp.Shape]] = mp.Simplex
 
     @property
     @memoize_method
@@ -407,10 +428,10 @@ def is_affine(self) -> bool:
 
 
 @dataclass(frozen=True, eq=False)
-class TensorProductElementGroup(_ModepyElementGroup):
+class TensorProductElementGroup(ModepyElementGroup):
     r"""Inherits from :class:`MeshElementGroup`."""
 
-    _modepy_shape_cls: ClassVar[type[mp.Shape]] = mp.Hypercube
+    shape_cls: ClassVar[type[mp.Shape]] = mp.Hypercube
 
     @property
     def is_affine(self) -> bool:
@@ -1472,8 +1493,8 @@ def __eq__(self, other: object) -> bool:
 # {{{ node-vertex consistency test
 
 def _mesh_group_node_vertex_error(mesh: Mesh, mgrp: MeshElementGroup) -> np.ndarray:
-    if isinstance(mgrp, _ModepyElementGroup):
-        basis = mp.basis_for_space(mgrp._modepy_space, mgrp._modepy_shape).functions
+    if isinstance(mgrp, ModepyElementGroup):
+        basis = mp.basis_for_space(mgrp.space, mgrp.shape).functions
     else:
         raise TypeError(f"unsupported group type: {type(mgrp).__name__}")
 
@@ -1535,7 +1556,7 @@ def _test_node_vertex_consistency(
     :raises InconsistentVerticesError: if the vertices are not consistent.
     """
     for igrp, mgrp in enumerate(mesh.groups):
-        if isinstance(mgrp, _ModepyElementGroup):
+        if isinstance(mgrp, ModepyElementGroup):
             _test_group_node_vertex_consistency_resampling(mesh, igrp, tol=tol)
         else:
             warn("Not implemented: node-vertex consistency check for "
@@ -2182,7 +2203,7 @@ def is_affine_simplex_group(
         return True
 
     # get matrices
-    basis = mp.basis_for_space(group._modepy_space, group._modepy_shape)
+    basis = mp.basis_for_space(group.space, group.shape)
     vinv = la.inv(mp.vandermonde(basis.functions, group.unit_nodes))
     diff = mp.differentiation_matrices(
             basis.functions, basis.gradients, group.unit_nodes)
@@ -2210,4 +2231,10 @@ def is_affine_simplex_group(
 
 # }}}
 
+
+__getattr__ = partial(module_getattr_for_deprecations, __name__, {
+        "_ModepyElementGroup": ("ModepyElementGroup", ModepyElementGroup, 2026),
+        })
+
+
 # vim: foldmethod=marker

meshmode/mesh/generation.py

@@ -46,8 +46,8 @@
 
 .. autofunction:: make_curve_mesh
 
-Curve parametrizations
-^^^^^^^^^^^^^^^^^^^^^^
+Curve parameterizations
+^^^^^^^^^^^^^^^^^^^^^^^
 
 .. autofunction:: circle
 .. autofunction:: ellipse
@@ -90,7 +90,7 @@
 """
 
 
-# {{{ test curve parametrizations
+# {{{ test curve parameterizations
 
 def circle(t: np.ndarray) -> np.ndarray:
     """
@@ -1666,7 +1666,7 @@ def warp_and_refine_until_resolved(
                     n_tail_orders=1 if warped_mesh.dim > 1 else 2)
 
             basis = mp.orthonormal_basis_for_space(
-                egrp._modepy_space, egrp._modepy_shape)
+                egrp.space, egrp.shape)
             vdm_inv = la.inv(mp.vandermonde(basis.functions, egrp.unit_nodes))
 
             mapping_coeffs = np.einsum("ij,dej->dei", vdm_inv, egrp.nodes)

meshmode/mesh/processing.py

@@ -571,9 +571,9 @@ def find_volume_mesh_element_group_orientation(
         each negatively oriented element.
     """
 
-    from meshmode.mesh import _ModepyElementGroup
+    from meshmode.mesh import ModepyElementGroup
 
-    if not isinstance(grp, _ModepyElementGroup):
+    if not isinstance(grp, ModepyElementGroup):
         raise NotImplementedError(
                 "finding element orientations "
                 "only supported on "
@@ -756,7 +756,7 @@ def _get_tensor_product_element_flip_matrix_and_vertex_permutation(
 
     flipped_unit_nodes = np.einsum("ij,jn->in", unit_flip_matrix, grp.unit_nodes)
 
-    basis = mp.basis_for_space(grp._modepy_space, grp._modepy_shape)
+    basis = mp.basis_for_space(grp.space, grp.shape)
     flip_matrix = mp.resampling_matrix(
         basis.functions,
         flipped_unit_nodes,