warp_and_refine_until_resolved: support tensor product elements

meshmode/discretization/poly_element.py

@@ -564,9 +564,7 @@ def quadrature_rule(self) -> mp.Quadrature:
         quads: list[mp.Quadrature] = []
 
         if self.dim != 1:
-            nodes_tp = cast(
-                "ArrayF",
-                reshape_array_for_tensor_product_space(self.space, self._nodes))
+            nodes_tp = reshape_array_for_tensor_product_space(self.space, self._nodes)
         else:
             nodes_tp = self._nodes
 

meshmode/mesh/generation.py

@@ -33,7 +33,7 @@
 import modepy as mp
 from pytools import log_process
 
-from meshmode.mesh import Mesh, MeshElementGroup, make_mesh
+from meshmode.mesh import Mesh, MeshElementGroup, ModepyElementGroup, make_mesh
 
 
 if TYPE_CHECKING:
@@ -1856,9 +1856,7 @@ def warp_and_refine_until_resolved(
     .. versionadded:: 2018.1
     """
     import modepy as mp
-    from modepy.modal_decay import simplex_interp_error_coefficient_estimator_matrix
 
-    from meshmode.mesh import SimplexElementGroup
     from meshmode.mesh.refinement import RefinerWithoutAdjacency
 
     if isinstance(unwarped_mesh_or_refiner, RefinerWithoutAdjacency):
@@ -1893,18 +1891,18 @@ def warp_and_refine_until_resolved(
         for base_element_nr, egrp in zip(
                 warped_mesh.base_element_nrs, warped_mesh.groups,
                 strict=True):
-            if not isinstance(egrp, SimplexElementGroup):
-                raise TypeError(
-                    f"Unsupported element group type: '{type(egrp).__name__}'")
+            n_tail_orders = 1 if warped_mesh.dim > 1 else 2
 
-            interp_err_est_mat = simplex_interp_error_coefficient_estimator_matrix(
-                    egrp.unit_nodes, egrp.order,
-                    n_tail_orders=1 if warped_mesh.dim > 1 else 2)
+            if not isinstance(egrp, ModepyElementGroup):
+                raise TypeError("only modepy element groups are supported")
 
-            basis = mp.orthonormal_basis_for_space(
-                egrp.space, egrp.shape)
+            basis = mp.orthonormal_basis_for_space(egrp.space, egrp.shape)
             vdm_inv = la.inv(mp.vandermonde(basis.functions, egrp.unit_nodes))
 
+            from modepy.modal_decay import interp_error_coefficient_estimator_matrix
+            interp_err_est_mat = interp_error_coefficient_estimator_matrix(
+                                    egrp.unit_nodes, n_tail_orders, basis)
+
             mapping_coeffs = np.einsum("ij,dej->dei", vdm_inv, egrp.nodes)
             mapping_norm_2 = np.sqrt(np.sum(mapping_coeffs**2, axis=-1))
 

test/test_mesh.py

@@ -1624,6 +1624,43 @@ def test_urchin():
     mgen.generate_urchin(3, 2, 4, 1e-4)
 
 
+def test_warp_and_refine_until_resolved_tensor_product():
+    """Test that warp_and_refine_until_resolved works with TensorProductElementGroup."""
+    from dataclasses import replace
+
+    from meshmode.mesh import TensorProductElementGroup, make_mesh
+    from meshmode.mesh.generation import warp_and_refine_until_resolved
+
+    def warp_mesh(mesh):
+        def warp_coords(coords):
+            x, y = coords
+            new_x = x + 0.1 * np.sin(2 * np.pi * y)
+            new_y = y + 0.1 * np.sin(2 * np.pi * x)
+            return np.array([new_x, new_y])
+
+        groups = [
+            replace(grp, nodes=warp_coords(grp.nodes))
+            for grp in mesh.groups]
+
+        return make_mesh(
+            warp_coords(mesh.vertices),
+            groups,
+            node_vertex_consistency_tolerance=False,
+            is_conforming=mesh.is_conforming,
+        )
+
+    mesh = mgen.generate_regular_rect_mesh(
+        a=(0, 0), b=(1, 1), nelements_per_axis=(3, 3),
+        group_cls=TensorProductElementGroup, order=4)
+
+    result = warp_and_refine_until_resolved(mesh, warp_mesh, 1e-3)
+
+    # The result must still be tensor product elements
+    assert isinstance(result.groups[0], TensorProductElementGroup)
+    # The highly oscillatory warp should have triggered refinement
+    assert result.nelements > mesh.nelements
+
+
 if __name__ == "__main__":
     import sys
     if len(sys.argv) > 1: