Merge branch 'main' into pytato-array-context-transforms

Author: thomasgibson

.ci/install-downstream.sh

@@ -14,7 +14,7 @@ cd "$DOWNSTREAM_PROJECT"
 echo "*** $DOWNSTREAM_PROJECT version: $(git rev-parse --short HEAD)"
 
 transfer_requirements_git_urls ../requirements.txt ./requirements.txt
-sed -i "/egg=meshmode/ c git+file://$(readlink -f ..)#egg=meshmode" requirements.txt
+edit_requirements_txt_for_downstream_in_subdir
 
 export CONDA_ENVIRONMENT=.test-conda-env-py3.yml
 
@@ -23,10 +23,9 @@ export PYTEST_ADDOPTS="-k 'not (slowtest or octave or mpi)'"
 
 if [[ "$DOWNSTREAM_PROJECT" = "mirgecom" ]]; then
     # can't turn off MPI in mirgecom
-    sudo apt-get update
-    sudo apt-get install openmpi-bin libopenmpi-dev
     export CONDA_ENVIRONMENT=conda-env.yml
     export CISUPPORT_PARALLEL_PYTEST=no
+    echo "- mpi4py" >> "$CONDA_ENVIRONMENT"
 else
     sed -i "/mpi4py/ d" requirements.txt
 fi

README.rst

@@ -11,7 +11,29 @@ meshmode: High-Order Meshes and Discontinuous Function Spaces
     :alt: Python Package Index Release Page
     :target: https://pypi.org/project/meshmode/
 
+Meshmode provides the "boring bits" of high-order unstructured discretization,
+for simplices (triangles, tetrahedra) and tensor products (quads, hexahedra).
+Features:
+
+- 1/2/3D, line/surface/volume discretizations in each, curvilinear supported.
+- "Everything is a (separate) discretization." (mesh boundaries are, element surfaces are,
+  refined versions of the same mesh are) "Connections" transfer information
+  between discretizations.
+- Periodic connectivity.
+- Mesh partitioning (not just) for distributed execution (e.g. via MPI).
+- Interpolatory, quadrature (overintegration), and modal element-local discretizations.
+- Independent of execution environment (GPU/CPU, numpy, ...)
+  via `array contexts <https://github.com/inducer/arraycontext/>`__.
+- Simple mesh refinement (via bisection). Adjacency currently only
+  maintained if uniform.
+- Input from Gmsh, Visualization to Vtk (both high-order curvilinear).
+- Easy data exchange with `Firedrake <https://www.firedrakeproject.org/>`__.
+
+Meshmode emerged as the shared discretization layer for `pytential
+<https://github.com/inducer/pytential/>`__ (layer potentials) and `grudge
+<https://github.com/inducer/grudge>`__ (discontinous Galerkin).
+
+Places on the web related to meshmode:
+
 * `Source code on Github <https://github.com/inducer/meshmode>`_
 * `Documentation <https://documen.tician.de/meshmode>`_
-
-.. TODO

meshmode/discretization/visualization.py

@@ -1016,6 +1016,10 @@ def show_scalar_in_matplotlib_3d(self, field, **kwargs):
 
         nodes = self._vis_nodes_numpy()
         field = _resample_to_numpy(self.connection, self.vis_discr, field)
+        if not np.can_cast(field.dtype, float, "same_kind"):
+            raise ValueError(
+                    f"fields of dtype {field.dtype} are not supported: "
+                    "cannot be converted to float")
 
         assert nodes.shape[0] == self.vis_discr.ambient_dim
 

meshmode/mesh/generation.py

@@ -646,28 +646,61 @@ def generate_torus_and_cycle_vertices(
         r_major: float, r_minor: float,
         n_major: int = 20, n_minor: int = 10, order: int = 1,
         node_vertex_consistency_tolerance: Optional[Union[float, bool]] = None,
-        unit_nodes: Optional[np.ndarray] = None):
+        unit_nodes: Optional[np.ndarray] = None,
+        group_cls: Optional[type] = None,
+        ):
+    from meshmode.mesh import SimplexElementGroup, TensorProductElementGroup
+    if group_cls is None:
+        group_cls = SimplexElementGroup
+
     a = r_major
     b = r_minor
     u, v = np.mgrid[0:2*np.pi:2*np.pi/n_major, 0:2*np.pi:2*np.pi/n_minor]
 
     # https://web.archive.org/web/20160410151837/https://www.math.hmc.edu/~gu/curves_and_surfaces/surfaces/torus.html  # noqa
-    x = np.cos(u)*(a+b*np.cos(v))
-    y = np.sin(u)*(a+b*np.cos(v))
-    z = b*np.sin(v)
-    vertices = (np.vstack((x[np.newaxis], y[np.newaxis], z[np.newaxis]))
+    x = np.cos(u) * (a + b*np.cos(v))
+    y = np.sin(u) * (a + b*np.cos(v))
+    z = b * np.sin(v)
+
+    vertices = (
+            np.vstack((x[np.newaxis], y[np.newaxis], z[np.newaxis]))
             .transpose(0, 2, 1).copy().reshape(3, -1))
 
     def idx(i, j):
         return (i % n_major) + (j % n_minor) * n_major
-    vertex_indices = ([(idx(i, j), idx(i+1, j), idx(i, j+1))
-            for i in range(n_major) for j in range(n_minor)]
-            + [(idx(i+1, j), idx(i+1, j+1), idx(i, j+1))
-            for i in range(n_major) for j in range(n_minor)])
+
+    if issubclass(group_cls, SimplexElementGroup):
+        # NOTE: this makes two triangles from a the square like
+        #   (i, j+1)    (i+1, j+1)
+        #       o---------o
+        #       | \       |
+        #       |   \     |
+        #       |     \   |
+        #       |       \ |
+        #       o---------o
+        #   (i, j)      (i+1, j)
+
+        vertex_indices = ([
+            (idx(i, j), idx(i+1, j), idx(i, j+1))
+            for i in range(n_major) for j in range(n_minor)
+            ] + [
+            (idx(i+1, j), idx(i+1, j+1), idx(i, j+1))
+            for i in range(n_major) for j in range(n_minor)
+            ])
+    elif issubclass(group_cls, TensorProductElementGroup):
+        # NOTE: this should match the order of the points in modepy
+        vertex_indices = [
+            (idx(i, j), idx(i+1, j), idx(i, j+1), idx(i+1, j+1))
+            for i in range(n_major) for j in range(n_minor)
+            ]
+    else:
+        raise TypeError(f"unsupported 'group_cls': {group_cls}")
 
     vertex_indices = np.array(vertex_indices, dtype=np.int32)
-    grp = make_group_from_vertices(vertices, vertex_indices, order,
-            unit_nodes=unit_nodes)
+    grp = make_group_from_vertices(
+            vertices, vertex_indices, order,
+            unit_nodes=unit_nodes,
+            group_cls=group_cls)
 
     # ambient_dim, nelements, nunit_nodes
     nodes = grp.nodes.copy()
@@ -693,12 +726,9 @@ def idx(i, j):
     x = np.sum(nodes*rvec, axis=0) - a
 
     minor_theta = np.arctan2(nodes[2], x)
-
-    nodes[0] = np.cos(major_theta)*(a+b*np.cos(
-        minor_theta))
-    nodes[1] = np.sin(major_theta)*(a+b*np.cos(
-        minor_theta))
-    nodes[2] = b*np.sin(minor_theta)
+    nodes[0] = np.cos(major_theta) * (a + b*np.cos(minor_theta))
+    nodes[1] = np.sin(major_theta) * (a + b*np.cos(minor_theta))
+    nodes[2] = b * np.sin(minor_theta)
 
     from meshmode.mesh import Mesh
     return (
@@ -716,7 +746,8 @@ def generate_torus(
         r_major: float, r_minor: float,
         n_major: int = 20, n_minor: int = 10, order: int = 1,
         node_vertex_consistency_tolerance: Optional[Union[float, bool]] = None,
-        unit_nodes: Optional[np.ndarray] = None):
+        unit_nodes: Optional[np.ndarray] = None,
+        group_cls: Optional[type] = None):
     r"""Generate a torus.
 
     .. figure:: images/torus.png
@@ -755,10 +786,12 @@ def generate_torus(
     :returns: a :class:`~meshmode.mesh.Mesh` of a torus.
 
     """
+
     mesh, _, _ = generate_torus_and_cycle_vertices(
             r_major, r_minor, n_major, n_minor, order,
             node_vertex_consistency_tolerance=node_vertex_consistency_tolerance,
-            unit_nodes=unit_nodes)
+            unit_nodes=unit_nodes,
+            group_cls=group_cls)
 
     return mesh
 

test/test_mesh.py

@@ -879,7 +879,7 @@ def test_quad_mesh_3d(mesh_name, order=3, visualize=False):
 # {{{ test_cube_icosahedron
 
 @pytest.mark.parametrize("order", [2, 3])
-def test_cube_icosphere(actx_factory, order, visualize=True):
+def test_cube_icosphere(actx_factory, order, visualize=False):
     mesh = mgen.generate_sphere(
             r=1.0, order=order,
             group_cls=TensorProductElementGroup,
@@ -898,6 +898,29 @@ def test_cube_icosphere(actx_factory, order, visualize=True):
 # }}}
 
 
+# {{{ test_tensor_torus
+
+@pytest.mark.parametrize("order", [3, 4])
+def test_tensor_torus(actx_factory, order, visualize=False):
+    mesh = mgen.generate_torus(
+            r_major=10.0, r_minor=5,
+            n_major=24, n_minor=12,
+            order=order,
+            group_cls=TensorProductElementGroup,
+            )
+
+    if not visualize:
+        return
+
+    from meshmode.mesh.visualization import vtk_visualize_mesh
+    actx = actx_factory()
+    vtk_visualize_mesh(actx, mesh,
+            f"quad_torus_order_{order:03d}.vtu",
+            vtk_high_order=False, overwrite=True)
+
+# }}}
+
+
 # {{{ mesh boundary gluing
 
 @pytest.mark.parametrize("use_tree", [False, True])

test/test_visualization.py

@@ -184,7 +184,8 @@ def test_visualizers(actx_factory, dim, group_cls):
 
     if mesh.dim == 2 and is_simplex:
         try:
-            vis.show_scalar_in_matplotlib_3d(f, do_show=False)
+            # NOTE: matplotlib only supports real fields
+            vis.show_scalar_in_matplotlib_3d(actx.np.real(f), do_show=False)
         except ImportError:
             logger.info("matplotlib not available")