Refiner: Make 1d element refinement work

meshmode/mesh/refinement/__init__.py

@@ -73,17 +73,19 @@ class Refiner(object):
     # {{{ constructor
 
     def __init__(self, mesh):
-        from meshmode.mesh.tesselate import tesselatetet, tesselatetri
+        from meshmode.mesh.tesselate import tesselateseg, tesselatetet, tesselatetri
         self.lazy = False
         self.seen_tuple = {}
         self.group_refinement_records = []
+        seg_node_tuples, seg_result = tesselateseg()
         tri_node_tuples, tri_result = tesselatetri()
         tet_node_tuples, tet_result = tesselatetet()
         #quadrilateral_node_tuples = [
         #print tri_result, tet_result
-        self.simplex_node_tuples = [None, None, tri_node_tuples, tet_node_tuples]
+        self.simplex_node_tuples = [
+                None, seg_node_tuples, tri_node_tuples, tet_node_tuples]
         # Dimension-parameterized tesselations for refinement
-        self.simplex_result = [None, None, tri_result, tet_result]
+        self.simplex_result = [None, seg_result, tri_result, tet_result]
         #print tri_node_tuples, tri_result
         #self.simplex_node_tuples, self.simplex_result = tesselatetet()
         self.last_mesh = mesh
@@ -594,8 +596,7 @@ class Refiner(object):
                     midpoint_vertices = []
                     vertex_indices = grp.vertex_indices[iel_grp]
                     #if simplex
-                    if len(grp.vertex_indices[iel_grp]) == len(self.last_mesh.vertices)+1:
-
+                    if len(vertex_indices) == grp.dim + 1:
                         # {{{ Get midpoints for all pairs of vertices
 
                         for i in range(len(vertex_indices)):

meshmode/mesh/refinement/utils.py

@@ -80,7 +80,7 @@ def check_nodal_adj_against_geometry(mesh, tol=1e-12):
                     nearby_origin_vertex = mesh.vertices[
                             :, nearby_grp.vertex_indices[nearby_iel][0]]  # noqa
                     transformation = np.empty(
-                            (len(mesh.vertices), nvertices_per_element-1))
+                            (mesh.ambient_dim, mesh.ambient_dim))
                     vertex_transformed = vertex - nearby_origin_vertex
 
                     for inearby_vertex_index, nearby_vertex_index in enumerate(

meshmode/mesh/tesselate.py

@@ -6,6 +6,12 @@ def add_tuples(a, b):
     return tuple(ac+bc for ac, bc in zip(a, b))
 
 
+def tesselateseg():
+    node_tuples = [(0,), (1,), (2,)]
+    result = [(0, 1), (1, 2)]
+    return [node_tuples, result]
+
+
 def tesselatetri():
     result = []
 

test/test_refinement.py

@@ -31,7 +31,7 @@ from pyopencl.tools import (  # noqa
         pytest_generate_tests_for_pyopencl
         as pytest_generate_tests)
 from meshmode.mesh.generation import (  # noqa
-        generate_icosahedron, generate_box_mesh)
+        generate_icosahedron, generate_box_mesh, make_curve_mesh, ellipse)
 from meshmode.mesh.refinement.utils import check_nodal_adj_against_geometry
 from meshmode.mesh.refinement import Refiner
 
@@ -80,6 +80,15 @@ def uniform_refine_flags(mesh):
     #     partial(random_refine_flags, 0.4),
     #     3),
 
+    ("3_to_1_ellipse_unif",
+        partial(
+            make_curve_mesh,
+            partial(ellipse, 3),
+            np.linspace(0, 1, 20),
+            order=2),
+        uniform_refine_flags,
+        4),
+
     ("rect2d_rand",
         partial(generate_box_mesh, (
             np.linspace(0, 1, 3),
@@ -138,6 +147,7 @@ def test_refinement(case_name, mesh_gen, flag_gen, num_generations):
     PolynomialEquidistantGroupFactory
     ])
 @pytest.mark.parametrize(("mesh_name", "dim", "mesh_pars"), [
+    ("circle", 1, [20, 30, 40]),
     ("blob", 2, [1e-1, 8e-2, 5e-2]),
     ("warp", 2, [4, 5, 6]),
     ("warp", 3, [4, 5, 6]),
@@ -171,7 +181,12 @@ def test_refinement_connection(
     for mesh_par in mesh_pars:
         # {{{ get mesh
 
-        if mesh_name == "blob":
+        if mesh_name == "circle":
+            assert dim == 1
+            h = 1 / mesh_par
+            mesh = make_curve_mesh(
+                partial(ellipse, 1), np.linspace(0, 1, mesh_par + 1), order=1)
+        elif mesh_name == "blob":
             assert dim == 2
             h = mesh_par
             mesh = gen_blob_mesh(h)

[Matt Wala]

@inducer, I'm not 100% what this line is doing, but after changing it the tests pass. Can you take a look to see if my change makes sense?

[Andreas Klöckner]

That builds a transformation from (roughly) barycentric coordinates (minus the first node) to global space. What you've done here isn't quite right fix. Try changing the solve down below to a lstsq(...)[0].

[Matt Wala]

Done in cb97c26a. When using lstsq we also need to check if the vertex lies in the space spanned by the element's vertices.