plot pder for the elements too

Author: alexfikl

test/test_symbolic.py

@@ -334,19 +334,32 @@ def test_node_reduction(actx_factory):
 
 # {{{ test_stretch_factor
 
+# {{{ twisted mesh
+
 def make_twisted_mesh(order, cls):
+    #  2       3             5
+    #   o------o-----------o
+    #   |      |           |
+    #   |      |           |
+    #   |      |           |
+    #   o------o-----------o
+    #  0       1            4
+    #
+    #
     vertices = np.array([
         [-1, -1, 0], [1, -1, 0], [-1, 1, 0], [1, 1, 0],
+        [4, -1, 0], [4, 1, 0],
         ], dtype=np.float64).T
 
     import meshmode.mesh as mm
     if issubclass(cls, mm.SimplexElementGroup):
         vertex_indices = np.array([
             (0, 1, 2), (1, 3, 2),
+            (1, 4, 3), (4, 5, 3),
             ], dtype=np.int32)
     elif issubclass(cls, mm.TensorProductElementGroup):
         vertex_indices = np.array([
-            (0, 1, 2, 3),
+            (0, 1, 2, 3), (1, 4, 3, 5),
             ], dtype=np.int32)
     else:
         raise ValueError
@@ -358,13 +371,10 @@ def make_twisted_mesh(order, cls):
             group_cls=cls)
 
     def wobble(x):
-        rx, ry = 2, 0.5
-        theta = np.pi / 4
-
         result = np.empty_like(x)
-        result[0] = rx * (np.cos(theta) * x[0] - np.sin(theta) * x[1])
-        result[1] = np.sin(ry * (np.sin(theta) * x[0] + np.cos(theta) * x[1]))
-        result[2] = x[2]
+        result[0] = x[0] + 0.5 * np.sin(x[1])
+        result[1] = x[1] + 0.5 * np.cos(x[0])
+        result[2] = x[2] + 0.5 * np.cos(x[1]) + 0.5 * np.sin(x[0])
         # result[2] = np.sin(x[1]) * np.sin(x[0])
 
         return result
@@ -373,6 +383,10 @@ def wobble(x):
     mesh = mm.Mesh(vertices, [grp], is_conforming=True)
     return map_mesh(mesh, wobble)
 
+# }}}
+
+
+# {{{ torus elements
 
 def make_torus_mesh(order, cls, a=2.0, b=1.0, n_major=12, n_minor=6):
     u, v = np.mgrid[0:2*np.pi:2*np.pi/n_major, 0:2*np.pi:2*np.pi/n_minor]
@@ -421,6 +435,10 @@ def idx(i, j):
 
     return mm.Mesh(vertices, [grp.copy(nodes=nodes)], is_conforming=True)
 
+# }}}
+
+
+# {{{ gmsh sphere
 
 def make_gmsh_sphere(order: int, cls: type):
     from meshmode.mesh.io import ScriptSource
@@ -460,6 +478,10 @@ def make_gmsh_sphere(order: int, cls: type):
             target_unit="MM",
             )
 
+# }}}
+
+
+# {{{ gmsh torus
 
 def make_gmsh_torus(order: int, cls: type):
     from meshmode.mesh.io import ScriptSource
@@ -499,6 +521,10 @@ def make_gmsh_torus(order: int, cls: type):
             target_unit="MM",
             )
 
+# }}}
+
+
+# {{{ symbolic
 
 def metric_from_form1(form1, metric_type: str):
     from pytential.symbolic.primitives import _small_mat_eigenvalues
@@ -540,10 +566,12 @@ def make_quad_stretch_factors(ambient_dim: int, metric_type: str):
 
     return metric_from_form1(form1, metric_type)
 
+# }}}
+
 
 @pytest.mark.parametrize("order", [4, 8])
 def test_stretch_factor(actx_factory, order,
-        mesh_name="torus", metric_type="eigs",
+        mesh_name="torus", metric_type="singvals",
         visualize=False):
     logging.basicConfig(level=logging.INFO)
     actx = actx_factory()
@@ -598,7 +626,8 @@ def print_bounds(x, name):
                 name, i,
                 actx.to_numpy(actx.np.min(si))[()],
                 actx.to_numpy(actx.np.min(si))[()]
-                ))
+                ), end=" ")
+        print()
 
     print_bounds(s, "s")
     print_bounds(q, "q")
@@ -610,15 +639,20 @@ def print_bounds(x, name):
 
     # {{{ plot vtk
 
+    s_pder = bind(simplex_discr, sym.parametrization_derivative_matrix(3, 2))(actx)
+    q_pder = bind(quad_discr, sym.parametrization_derivative_matrix(3, 2))(actx)
+
     from meshmode.discretization.visualization import make_visualizer
     vis = make_visualizer(actx, simplex_discr, order, force_equidistant=True)
     vis.write_vtk_file(f"simplex_{suffix}.vtu",
-            [(f"s{i}", si) for i, si in enumerate(s)],
+            [(f"s{i}", si) for i, si in enumerate(s)]
+            + [(f"J_{i}_{j}", pder) for (i, j), pder in np.ndenumerate(s_pder)],
             overwrite=True, use_high_order=True)
 
     vis = make_visualizer(actx, quad_discr, order, force_equidistant=True)
     vis.write_vtk_file(f"quad_{suffix}.vtu",
-            [(f"q{i}", qi) for i, qi in enumerate(q)],
+            [(f"q{i}", qi) for i, qi in enumerate(q)]
+            + [(f"J_{i}_{j}", pder) for (i, j), pder in np.ndenumerate(q_pder)],
             overwrite=True, use_high_order=True)
 
     # }}}