From 934ec897f621134645cf8a774b65a46066f5eabe Mon Sep 17 00:00:00 2001
From: Ellis <eshoag2@illinois.edu>
Date: Fri, 10 Feb 2017 16:59:22 -0600
Subject: [PATCH 01/12] Initial partition_mesh code

---
 meshmode/mesh/processing.py | 51 +++++++++++++++++++++++++++++++++++++
 1 file changed, 51 insertions(+)

diff --git a/meshmode/mesh/processing.py b/meshmode/mesh/processing.py
index ff74ac81..27d3ef6e 100644
--- a/meshmode/mesh/processing.py
+++ b/meshmode/mesh/processing.py
@@ -39,6 +39,57 @@ __doc__ = """
 .. autofunction:: affine_map
 """
 
+def partition_mesh(mesh, part_per_element, part_nr):
+    """
+    :arg mesh: A :class:`meshmode.mesh.Mesh` to be partitioned.
+    :arg part_per_element: A :class:`numpy.ndarray` containing one
+        integer per element of *mesh* indicating which part of the
+        partitioned mesh the element is to become a part of.
+    :arg part_nr: The part number of the mesh to return.
+
+    :returns: A tuple ``(part_mesh, part_to_global)``, where *part_mesh*
+        is a :class:`meshmode.mesh.Mesh` that is a partition of mesh, and
+        *part_to_global* is a :class:`numpy.ndarray` mapping element
+        numbers on *part_mesh* to ones in *mesh*.
+    """
+
+    queried_elems = np.where(part_per_element == part_nr)[0]
+
+    '''
+    Here I assume that all elements are taken from the 0th group.
+    Will there be more groups? How do I handle this?
+    My guess is that len(part_per_element) will be equal to the sum of the
+    number of elements in each group. And that the first element of groups[1]
+    comes just after the last element in groups[0].
+    '''
+    mesh_group = mesh.groups[0]
+
+    new_vertex_indices = mesh_group.vertex_indices[queried_elems]
+
+    # A sorted np.array of vertex indices we need in our new mesh (without duplicates).
+    required_vertex_indices = np.array(list(set(new_vertex_indices.ravel())))
+
+    new_vertices = np.zeros((mesh.ambient_dim, len(required_vertex_indices)))
+    for dim in range(mesh.ambient_dim):
+        new_vertices[dim] = mesh.vertices[dim][required_vertex_indices]
+
+    # We need to update our indices to be in range [0, len(required_vertex_indices)].
+    for i in range(len(new_vertex_indices)):
+        for j in range(len(new_vertex_indices[0])):
+            new_vertex_indices[i, j] = np.where(required_vertex_indices == new_vertex_indices[i, j])[0]
+
+    new_nodes = np.zeros((mesh.ambient_dim, len(queried_elems), mesh_group.nunit_nodes))
+    for i in range(mesh.ambient_dim):
+        for j in range(len(queried_elems)):
+            new_nodes[i, j, :] = mesh_group.nodes[i, queried_elems[j], :]
+
+    from meshmode.mesh import MeshElementGroup, Mesh
+
+    mesh_element_group = MeshElementGroup(mesh_group.order, new_vertex_indices, new_nodes, element_nr_base=mesh_group.element_nr_base, node_nr_base=mesh_group.node_nr_base, unit_nodes=mesh_group.unit_nodes, dim=mesh_group.dim)
+
+    part_mesh = Mesh(new_vertices, [mesh_element_group])
+
+    return (part_mesh, queried_elems)
 
 # {{{ orientations
 
-- 
GitLab


From 0c1dc43bce315ae8c0cfaeb70ef8e5f01e9d55a5 Mon Sep 17 00:00:00 2001
From: Ellis <eshoag2@illinois.edu>
Date: Fri, 10 Feb 2017 18:51:55 -0600
Subject: [PATCH 02/12] added partition_mesh test

---
 test/test_meshmode.py | 23 +++++++++++++++++++++++
 1 file changed, 23 insertions(+)

diff --git a/test/test_meshmode.py b/test/test_meshmode.py
index 3dabe0cf..eeda06b5 100644
--- a/test/test_meshmode.py
+++ b/test/test_meshmode.py
@@ -47,6 +47,29 @@ import pytest
 import logging
 logger = logging.getLogger(__name__)
 
+# {{{ partition_mesh
+
+@pytest.mark.parametrize("mesh_type", ["cloverleaf", "starfish"])
+@pytest.mark.parametrize("npoints", [10, 1000])
+def test_partition_mesh(mesh_type, npoints):
+    from meshmode.mesh.generation import make_curve_mesh, cloverleaf, starfish
+
+    if mesh_type == "cloverleaf":
+        mesh = make_curve_mesh(cloverleaf, np.linspace(0, 1, npoints), order=3)
+    elif mesh_type == "starfish":
+        mesh = make_curve_mesh(starfish, np.linspace(0, 1, npoints), order=3)
+
+    #TODO: Create an actuall adjacency list from the mesh.
+    adjacency_list = mesh.nodal_adjacency
+
+    from pymetis import part_graph
+    part_per_element = np.array(part_graph(3, adjacency=adjacency_list))
+
+    from meshmode.mesh.processing import partition_mesh
+    (part_mesh, part_to_global) = partition_mesh(mesh, part_per_element, 0)
+
+# }}}
+
 
 # {{{ circle mesh
 
-- 
GitLab


From 90346bedeec0398a1952a511a6958b7e4f26b1cf Mon Sep 17 00:00:00 2001
From: Ellis <eshoag2@illinois.edu>
Date: Fri, 10 Feb 2017 20:01:56 -0600
Subject: [PATCH 03/12] partition mesh almost supports multiple element groups

---
 meshmode/mesh/processing.py | 64 +++++++++++++++++++++++--------------
 test/test_meshmode.py       | 17 ++++++++--
 2 files changed, 54 insertions(+), 27 deletions(-)

diff --git a/meshmode/mesh/processing.py b/meshmode/mesh/processing.py
index 27d3ef6e..7f520f44 100644
--- a/meshmode/mesh/processing.py
+++ b/meshmode/mesh/processing.py
@@ -31,6 +31,7 @@ import modepy as mp
 
 
 __doc__ = """
+.. autofunction:: partition_mesh
 .. autofunction:: find_volume_mesh_element_orientations
 .. autofunction:: perform_flips
 .. autofunction:: find_bounding_box
@@ -39,6 +40,7 @@ __doc__ = """
 .. autofunction:: affine_map
 """
 
+
 def partition_mesh(mesh, part_per_element, part_nr):
     """
     :arg mesh: A :class:`meshmode.mesh.Mesh` to be partitioned.
@@ -52,45 +54,59 @@ def partition_mesh(mesh, part_per_element, part_nr):
         *part_to_global* is a :class:`numpy.ndarray` mapping element
         numbers on *part_mesh* to ones in *mesh*.
     """
+    assert len(part_per_element) == mesh.nelements
+
+    queried_elems = np.where(np.array(part_per_element) == part_nr)[0]
+
+    num_groups = len(mesh.groups)
+    new_indices = []
+    new_nodes = []
+
+    # The set of vertex indices we need.
+    index_set = set()
 
-    queried_elems = np.where(part_per_element == part_nr)[0]
+    start_idx = 0
+    for group_nr in range(num_groups):
+        mesh_group = mesh.groups[group_nr]
+        #end_idx = start_idx + np.where(queried_elems >= mesh_group.nelements - start_idx)[0] - 1
+        end_idx = len(queried_elems)
 
-    '''
-    Here I assume that all elements are taken from the 0th group.
-    Will there be more groups? How do I handle this?
-    My guess is that len(part_per_element) will be equal to the sum of the
-    number of elements in each group. And that the first element of groups[1]
-    comes just after the last element in groups[0].
-    '''
-    mesh_group = mesh.groups[0]
+        new_indices.append(mesh_group.vertex_indices[queried_elems[start_idx:end_idx]])
 
-    new_vertex_indices = mesh_group.vertex_indices[queried_elems]
+        new_nodes.append(np.zeros((mesh.ambient_dim, len(queried_elems), mesh_group.nunit_nodes)))
+        for i in range(mesh.ambient_dim):
+            for j in range(start_idx, end_idx):
+                new_nodes[group_nr][i, j, :] = mesh_group.nodes[i, queried_elems[j], :]
+
+        index_set |= set(new_indices[group_nr].ravel())
+
+        start_idx = end_idx
 
     # A sorted np.array of vertex indices we need in our new mesh (without duplicates).
-    required_vertex_indices = np.array(list(set(new_vertex_indices.ravel())))
+    required_indices = np.array(list(index_set))
 
-    new_vertices = np.zeros((mesh.ambient_dim, len(required_vertex_indices)))
+    new_vertices = np.zeros((mesh.ambient_dim, len(required_indices)))
     for dim in range(mesh.ambient_dim):
-        new_vertices[dim] = mesh.vertices[dim][required_vertex_indices]
+        new_vertices[dim] = mesh.vertices[dim][required_indices]
 
-    # We need to update our indices to be in range [0, len(required_vertex_indices)].
-    for i in range(len(new_vertex_indices)):
-        for j in range(len(new_vertex_indices[0])):
-            new_vertex_indices[i, j] = np.where(required_vertex_indices == new_vertex_indices[i, j])[0]
-
-    new_nodes = np.zeros((mesh.ambient_dim, len(queried_elems), mesh_group.nunit_nodes))
-    for i in range(mesh.ambient_dim):
-        for j in range(len(queried_elems)):
-            new_nodes[i, j, :] = mesh_group.nodes[i, queried_elems[j], :]
+    # We need to update our indices to be in range [0, len(mesh_group.nelements)].
+    for group_nr in range(num_groups):
+        for i in range(len(new_indices[group_nr])):
+            for j in range(len(new_indices[group_nr][0])):
+                new_indices[group_nr][i, j] = np.where(required_indices == new_indices[group_nr][i, j])[0]
 
     from meshmode.mesh import MeshElementGroup, Mesh
 
-    mesh_element_group = MeshElementGroup(mesh_group.order, new_vertex_indices, new_nodes, element_nr_base=mesh_group.element_nr_base, node_nr_base=mesh_group.node_nr_base, unit_nodes=mesh_group.unit_nodes, dim=mesh_group.dim)
+    new_mesh_groups = []
+    for group_nr in range(num_groups):
+        mesh_group = mesh.groups[group_nr]
+        new_mesh_groups.append(MeshElementGroup(mesh_group.order, new_indices[group_nr], new_nodes[group_nr], unit_nodes=mesh_group.unit_nodes, dim=mesh_group.dim))
 
-    part_mesh = Mesh(new_vertices, [mesh_element_group])
+    part_mesh = Mesh(new_vertices, new_mesh_groups)
 
     return (part_mesh, queried_elems)
 
+
 # {{{ orientations
 
 def find_volume_mesh_element_group_orientation(vertices, grp):
diff --git a/test/test_meshmode.py b/test/test_meshmode.py
index eeda06b5..b5b2cde3 100644
--- a/test/test_meshmode.py
+++ b/test/test_meshmode.py
@@ -51,8 +51,19 @@ logger = logging.getLogger(__name__)
 
 @pytest.mark.parametrize("mesh_type", ["cloverleaf", "starfish"])
 @pytest.mark.parametrize("npoints", [10, 1000])
-def test_partition_mesh(mesh_type, npoints):
-    from meshmode.mesh.generation import make_curve_mesh, cloverleaf, starfish
+def test_partition_mesh(mesh_type=None, npoints=None):
+    from meshmode.mesh.generation import generate_torus
+
+    my_mesh = generate_torus(2, 1, n_outer=2, n_inner=2)
+
+    part_per_element = np.array([0, 1, 2, 1, 1, 2, 1, 0])
+
+    print(my_mesh.nelements, my_mesh.groups[0].nelements)
+
+    #(part_mesh, part_to_global) = partition_mesh(my_mesh, part_per_element, 1)
+    from meshmode.mesh.processing import partition_mesh
+    partition_mesh(my_mesh, part_per_element, 1)
+    """from meshmode.mesh.generation import make_curve_mesh, cloverleaf, starfish
 
     if mesh_type == "cloverleaf":
         mesh = make_curve_mesh(cloverleaf, np.linspace(0, 1, npoints), order=3)
@@ -67,7 +78,7 @@ def test_partition_mesh(mesh_type, npoints):
 
     from meshmode.mesh.processing import partition_mesh
     (part_mesh, part_to_global) = partition_mesh(mesh, part_per_element, 0)
-
+"""
 # }}}
 
 
-- 
GitLab


From 28691b122cf04ec22250131b29a278899d8198f0 Mon Sep 17 00:00:00 2001
From: Ellis <eshoag2@illinois.edu>
Date: Sun, 12 Feb 2017 12:15:39 -0600
Subject: [PATCH 04/12] partition_mesh handles multiple MeshElementGroups

---
 meshmode/mesh/processing.py | 39 +++++++++++++++++++++++++++++--------
 test/test_meshmode.py       | 28 +++++++++++++++-----------
 2 files changed, 48 insertions(+), 19 deletions(-)

diff --git a/meshmode/mesh/processing.py b/meshmode/mesh/processing.py
index 7f520f44..05ce910e 100644
--- a/meshmode/mesh/processing.py
+++ b/meshmode/mesh/processing.py
@@ -56,6 +56,7 @@ def partition_mesh(mesh, part_per_element, part_nr):
     """
     assert len(part_per_element) == mesh.nelements
 
+    # Contains the indices of the elements requested.
     queried_elems = np.where(np.array(part_per_element) == part_nr)[0]
 
     num_groups = len(mesh.groups)
@@ -68,39 +69,61 @@ def partition_mesh(mesh, part_per_element, part_nr):
     start_idx = 0
     for group_nr in range(num_groups):
         mesh_group = mesh.groups[group_nr]
-        #end_idx = start_idx + np.where(queried_elems >= mesh_group.nelements - start_idx)[0] - 1
+
+        # Find the index of first element in the next group
         end_idx = len(queried_elems)
+        for idx in range(len(queried_elems)):
+            if queried_elems[idx] - start_idx > mesh_group.nelements:
+                end_idx = idx
+                break
 
-        new_indices.append(mesh_group.vertex_indices[queried_elems[start_idx:end_idx]])
+        elems = queried_elems[start_idx:end_idx]
+        new_indices.append(mesh_group.vertex_indices[elems])
 
-        new_nodes.append(np.zeros((mesh.ambient_dim, len(queried_elems), mesh_group.nunit_nodes)))
+        new_nodes.append(
+            np.zeros((mesh.ambient_dim, len(queried_elems), mesh_group.nunit_nodes)))
         for i in range(mesh.ambient_dim):
             for j in range(start_idx, end_idx):
-                new_nodes[group_nr][i, j, :] = mesh_group.nodes[i, queried_elems[j], :]
+                elems = queried_elems[j]
+                new_nodes[group_nr][i, j, :] = mesh_group.nodes[i, elems, :]
 
         index_set |= set(new_indices[group_nr].ravel())
 
         start_idx = end_idx
 
-    # A sorted np.array of vertex indices we need in our new mesh (without duplicates).
+    # A sorted np.array of vertex indices we need (without duplicates).
     required_indices = np.array(list(index_set))
 
     new_vertices = np.zeros((mesh.ambient_dim, len(required_indices)))
     for dim in range(mesh.ambient_dim):
         new_vertices[dim] = mesh.vertices[dim][required_indices]
 
-    # We need to update our indices to be in range [0, len(mesh_group.nelements)].
+    # Our indices need to be in range [0, len(mesh_group.nelements)].
     for group_nr in range(num_groups):
         for i in range(len(new_indices[group_nr])):
             for j in range(len(new_indices[group_nr][0])):
-                new_indices[group_nr][i, j] = np.where(required_indices == new_indices[group_nr][i, j])[0]
+                original_index = new_indices[group_nr][i, j]
+                new_indices[group_nr][i, j] = np.where(
+                    required_indices == original_index)[0]
+
+    """
+    print("mesh vertices: ", mesh.vertices)
+    print("mesh indices: ", mesh.groups[0].vertex_indices)
+    print("mesh nodes: ", mesh.groups[0].nodes)
+    print("queried_elems: ", queried_elems)
+    print("indices: ", new_indices[0])
+    print("nodes: ", new_nodes[0])
+    print("vertices: ", new_vertices)
+    """
 
     from meshmode.mesh import MeshElementGroup, Mesh
 
     new_mesh_groups = []
     for group_nr in range(num_groups):
         mesh_group = mesh.groups[group_nr]
-        new_mesh_groups.append(MeshElementGroup(mesh_group.order, new_indices[group_nr], new_nodes[group_nr], unit_nodes=mesh_group.unit_nodes, dim=mesh_group.dim))
+        new_mesh_groups.append(
+            MeshElementGroup(mesh_group.order, new_indices[group_nr],
+                new_nodes[group_nr], unit_nodes=mesh_group.unit_nodes))
 
     part_mesh = Mesh(new_vertices, new_mesh_groups)
 
diff --git a/test/test_meshmode.py b/test/test_meshmode.py
index b5b2cde3..99a9284e 100644
--- a/test/test_meshmode.py
+++ b/test/test_meshmode.py
@@ -47,38 +47,44 @@ import pytest
 import logging
 logger = logging.getLogger(__name__)
 
+
 # {{{ partition_mesh
 
-@pytest.mark.parametrize("mesh_type", ["cloverleaf", "starfish"])
-@pytest.mark.parametrize("npoints", [10, 1000])
-def test_partition_mesh(mesh_type=None, npoints=None):
+#@pytest.mark.parametrize("mesh_type", ["cloverleaf", "starfish"])
+#@pytest.mark.parametrize("npoints", [10, 1000])
+def test_partition_mesh(mesh_type, npoints):
     from meshmode.mesh.generation import generate_torus
-
     my_mesh = generate_torus(2, 1, n_outer=2, n_inner=2)
 
     part_per_element = np.array([0, 1, 2, 1, 1, 2, 1, 0])
 
-    print(my_mesh.nelements, my_mesh.groups[0].nelements)
-
-    #(part_mesh, part_to_global) = partition_mesh(my_mesh, part_per_element, 1)
     from meshmode.mesh.processing import partition_mesh
-    partition_mesh(my_mesh, part_per_element, 1)
-    """from meshmode.mesh.generation import make_curve_mesh, cloverleaf, starfish
+    (part_mesh, part_to_global) = partition_mesh(my_mesh, part_per_element, 1)
+
+
+    """
+    from meshmode.mesh.generation import make_curve_mesh, cloverleaf, starfish
 
     if mesh_type == "cloverleaf":
         mesh = make_curve_mesh(cloverleaf, np.linspace(0, 1, npoints), order=3)
     elif mesh_type == "starfish":
         mesh = make_curve_mesh(starfish, np.linspace(0, 1, npoints), order=3)
 
+    #from meshmode.mesh.visualization import draw_2d_mesh
+    #draw_2d_mesh(mesh)
+    #import matplotlib.pyplot as pt
+    #pt.show()
+
     #TODO: Create an actuall adjacency list from the mesh.
-    adjacency_list = mesh.nodal_adjacency
+    adjacency_list = mesh.facial_adjacency_groups
+    print(adjacency_list)
 
     from pymetis import part_graph
     part_per_element = np.array(part_graph(3, adjacency=adjacency_list))
 
     from meshmode.mesh.processing import partition_mesh
     (part_mesh, part_to_global) = partition_mesh(mesh, part_per_element, 0)
-"""
+    """
 # }}}
 
 
-- 
GitLab


From 9d674052d3ea669f0509995938c49f64604771f8 Mon Sep 17 00:00:00 2001
From: ellis <eshoag2@illinois.edu>
Date: Sun, 12 Feb 2017 18:02:18 -0600
Subject: [PATCH 05/12] partition_mesh fixes and tests

---
 meshmode/mesh/processing.py | 16 ++------
 test/test_meshmode.py       | 73 ++++++++++++++++++++-----------------
 2 files changed, 42 insertions(+), 47 deletions(-)

diff --git a/meshmode/mesh/processing.py b/meshmode/mesh/processing.py
index 05ce910e..02a7e3bc 100644
--- a/meshmode/mesh/processing.py
+++ b/meshmode/mesh/processing.py
@@ -54,7 +54,7 @@ def partition_mesh(mesh, part_per_element, part_nr):
         *part_to_global* is a :class:`numpy.ndarray` mapping element
         numbers on *part_mesh* to ones in *mesh*.
     """
-    assert len(part_per_element) == mesh.nelements
+    assert len(part_per_element) == mesh.nelements, "part_per_element must have shape (mesh.nelements,)"
 
     # Contains the indices of the elements requested.
     queried_elems = np.where(np.array(part_per_element) == part_nr)[0]
@@ -106,23 +106,13 @@ def partition_mesh(mesh, part_per_element, part_nr):
                 new_indices[group_nr][i, j] = np.where(
                     required_indices == original_index)[0]
 
-    """
-    print("mesh vertices: ", mesh.vertices)
-    print("mesh indices: ", mesh.groups[0].vertex_indices)
-    print("mesh nodes: ", mesh.groups[0].nodes)
-    print("queried_elems: ", queried_elems)
-    print("indices: ", new_indices[0])
-    print("nodes: ", new_nodes[0])
-    print("vertices: ", new_vertices)
-    """
-
-    from meshmode.mesh import MeshElementGroup, Mesh
+    from meshmode.mesh import SimplexElementGroup, Mesh
 
     new_mesh_groups = []
     for group_nr in range(num_groups):
         mesh_group = mesh.groups[group_nr]
         new_mesh_groups.append(
-            MeshElementGroup(mesh_group.order, new_indices[group_nr],
+            SimplexElementGroup(mesh_group.order, new_indices[group_nr],
                 new_nodes[group_nr], unit_nodes=mesh_group.unit_nodes))
 
     part_mesh = Mesh(new_vertices, new_mesh_groups)
diff --git a/test/test_meshmode.py b/test/test_meshmode.py
index 99a9284e..740b917b 100644
--- a/test/test_meshmode.py
+++ b/test/test_meshmode.py
@@ -50,41 +50,46 @@ logger = logging.getLogger(__name__)
 
 # {{{ partition_mesh
 
-#@pytest.mark.parametrize("mesh_type", ["cloverleaf", "starfish"])
-#@pytest.mark.parametrize("npoints", [10, 1000])
-def test_partition_mesh(mesh_type, npoints):
-    from meshmode.mesh.generation import generate_torus
-    my_mesh = generate_torus(2, 1, n_outer=2, n_inner=2)
+@pytest.mark.parameterize("mesh_type", ["torus", "box"])
+def test_partition_mesh(mesh_type):
+    if mesh_type == "torus":
+        from meshmode.mesh.generation import generate_torus
+        my_mesh = generate_torus(2, 1, n_outer=2, n_inner=2)
+
+        part_per_element = np.array([0, 1, 2, 1, 1, 2, 1, 0])
+
+        from meshmode.mesh.processing import partition_mesh
+        (part_mesh, part_to_global) = partition_mesh(my_mesh, part_per_element, 0)
+        assert part_mesh.nelements == 2
+
+        (part_mesh, part_to_global) = partition_mesh(my_mesh, part_per_element, 1)
+        assert part_mesh.nelements == 4
+
+        (part_mesh, part_to_global) = partition_mesh(my_mesh, part_per_element, 2)
+        assert part_mesh.nelements == 2
+    elif mesh_type == "box":
+        from meshmode.mesh.generation import generate_box_mesh
+        seg = np.linspace(0, 1, 10)
+        axis_coords = (seg, seg, seg)
+        mesh = generate_box_mesh(axis_coords)
+
+        adjacency_list = np.zeros((mesh.nelements,), dtype=set)
+        for elem in range(mesh.nelements):
+            adjacency_list[elem] = set()
+            for n in range(mesh.nodal_adjacency.neighbors_starts[elem], mesh.nodal_adjacency.neighbors_starts[elem + 1]):
+                adjacency_list[elem].add(mesh.nodal_adjacency.neighbors[n])
+
+        from pymetis import part_graph
+        (_, part_per_element) = part_graph(3, adjacency=adjacency_list)
+
+        from meshmode.mesh.processing import partition_mesh
+        (part_mesh0, _) = partition_mesh(mesh, np.array(part_per_element), 0)
+        (part_mesh1, _) = partition_mesh(mesh, np.array(part_per_element), 1)
+        (part_mesh2, _) = partition_mesh(mesh, np.array(part_per_element), 2)
+
+        assert mesh.nelements == (part_mesh0.nelements
+            + part_mesh1.nelements + part_mesh2.nelements)
 
-    part_per_element = np.array([0, 1, 2, 1, 1, 2, 1, 0])
-
-    from meshmode.mesh.processing import partition_mesh
-    (part_mesh, part_to_global) = partition_mesh(my_mesh, part_per_element, 1)
-
-
-    """
-    from meshmode.mesh.generation import make_curve_mesh, cloverleaf, starfish
-
-    if mesh_type == "cloverleaf":
-        mesh = make_curve_mesh(cloverleaf, np.linspace(0, 1, npoints), order=3)
-    elif mesh_type == "starfish":
-        mesh = make_curve_mesh(starfish, np.linspace(0, 1, npoints), order=3)
-
-    #from meshmode.mesh.visualization import draw_2d_mesh
-    #draw_2d_mesh(mesh)
-    #import matplotlib.pyplot as pt
-    #pt.show()
-
-    #TODO: Create an actuall adjacency list from the mesh.
-    adjacency_list = mesh.facial_adjacency_groups
-    print(adjacency_list)
-
-    from pymetis import part_graph
-    part_per_element = np.array(part_graph(3, adjacency=adjacency_list))
-
-    from meshmode.mesh.processing import partition_mesh
-    (part_mesh, part_to_global) = partition_mesh(mesh, part_per_element, 0)
-    """
 # }}}
 
 
-- 
GitLab


From 7f07c08cdc45f0e17b6b72eae6317b10c458836a Mon Sep 17 00:00:00 2001
From: ellis <eshoag2@illinois.edu>
Date: Sun, 12 Feb 2017 23:51:34 -0600
Subject: [PATCH 06/12] spelling

---
 meshmode/mesh/processing.py |  3 ++-
 test/test_meshmode.py       | 17 +++++++++--------
 2 files changed, 11 insertions(+), 9 deletions(-)

diff --git a/meshmode/mesh/processing.py b/meshmode/mesh/processing.py
index 02a7e3bc..00b7a8bf 100644
--- a/meshmode/mesh/processing.py
+++ b/meshmode/mesh/processing.py
@@ -54,7 +54,8 @@ def partition_mesh(mesh, part_per_element, part_nr):
         *part_to_global* is a :class:`numpy.ndarray` mapping element
         numbers on *part_mesh* to ones in *mesh*.
     """
-    assert len(part_per_element) == mesh.nelements, "part_per_element must have shape (mesh.nelements,)"
+    assert len(part_per_element) == mesh.nelements, (
+        "part_per_element must have shape (mesh.nelements,)")
 
     # Contains the indices of the elements requested.
     queried_elems = np.where(np.array(part_per_element) == part_nr)[0]
diff --git a/test/test_meshmode.py b/test/test_meshmode.py
index 740b917b..4cdffd85 100644
--- a/test/test_meshmode.py
+++ b/test/test_meshmode.py
@@ -50,7 +50,7 @@ logger = logging.getLogger(__name__)
 
 # {{{ partition_mesh
 
-@pytest.mark.parameterize("mesh_type", ["torus", "box"])
+@pytest.mark.parametrize("mesh_type", ["torus", "box"])
 def test_partition_mesh(mesh_type):
     if mesh_type == "torus":
         from meshmode.mesh.generation import generate_torus
@@ -59,14 +59,14 @@ def test_partition_mesh(mesh_type):
         part_per_element = np.array([0, 1, 2, 1, 1, 2, 1, 0])
 
         from meshmode.mesh.processing import partition_mesh
-        (part_mesh, part_to_global) = partition_mesh(my_mesh, part_per_element, 0)
-        assert part_mesh.nelements == 2
+        (part_mesh0, _) = partition_mesh(my_mesh, part_per_element, 0)
+        (part_mesh1, _) = partition_mesh(my_mesh, part_per_element, 1)
+        (part_mesh2, _) = partition_mesh(my_mesh, part_per_element, 2)
 
-        (part_mesh, part_to_global) = partition_mesh(my_mesh, part_per_element, 1)
-        assert part_mesh.nelements == 4
+        assert part_mesh0.nelements == 2
+        assert part_mesh1.nelements == 4
+        assert part_mesh2.nelements == 2
 
-        (part_mesh, part_to_global) = partition_mesh(my_mesh, part_per_element, 2)
-        assert part_mesh.nelements == 2
     elif mesh_type == "box":
         from meshmode.mesh.generation import generate_box_mesh
         seg = np.linspace(0, 1, 10)
@@ -76,7 +76,8 @@ def test_partition_mesh(mesh_type):
         adjacency_list = np.zeros((mesh.nelements,), dtype=set)
         for elem in range(mesh.nelements):
             adjacency_list[elem] = set()
-            for n in range(mesh.nodal_adjacency.neighbors_starts[elem], mesh.nodal_adjacency.neighbors_starts[elem + 1]):
+            starts = mesh.nodal_adjacency.neighbors_starts
+            for n in range(starts[elem], starts[elem + 1]):
                 adjacency_list[elem].add(mesh.nodal_adjacency.neighbors[n])
 
         from pymetis import part_graph
-- 
GitLab


From 68c870eec9e622a41e4579c92bc287acfac6eaa5 Mon Sep 17 00:00:00 2001
From: ellis <eshoag2@illinois.edu>
Date: Tue, 14 Feb 2017 12:20:53 -0600
Subject: [PATCH 07/12] partition_mesh works with multiple MeshElementGroups

---
 meshmode/mesh/processing.py | 45 +++++++++++++++++++++++--------------
 test/test_meshmode.py       | 23 ++++++++++---------
 2 files changed, 41 insertions(+), 27 deletions(-)

diff --git a/meshmode/mesh/processing.py b/meshmode/mesh/processing.py
index 00b7a8bf..28d146b7 100644
--- a/meshmode/mesh/processing.py
+++ b/meshmode/mesh/processing.py
@@ -66,30 +66,39 @@ def partition_mesh(mesh, part_per_element, part_nr):
 
     # The set of vertex indices we need.
     index_set = set()
-
+    skip_groups = []
+    num_prev_elems = 0
     start_idx = 0
     for group_nr in range(num_groups):
         mesh_group = mesh.groups[group_nr]
 
         # Find the index of first element in the next group
         end_idx = len(queried_elems)
-        for idx in range(len(queried_elems)):
-            if queried_elems[idx] - start_idx > mesh_group.nelements:
+        for idx in range(start_idx, len(queried_elems)):
+            if queried_elems[idx] - num_prev_elems >= mesh_group.nelements:
                 end_idx = idx
                 break
 
-        elems = queried_elems[start_idx:end_idx]
+        if start_idx == end_idx:
+            skip_groups.append(group_nr)
+            new_indices.append(np.array([]))
+            new_nodes.append(np.array([]))
+            num_prev_elems += mesh_group.nelements
+            continue
+
+        elems = queried_elems[start_idx:end_idx] - num_prev_elems
         new_indices.append(mesh_group.vertex_indices[elems])
 
         new_nodes.append(
-            np.zeros((mesh.ambient_dim, len(queried_elems), mesh_group.nunit_nodes)))
+            np.zeros((mesh.ambient_dim, end_idx - start_idx, mesh_group.nunit_nodes)))
         for i in range(mesh.ambient_dim):
             for j in range(start_idx, end_idx):
-                elems = queried_elems[j]
-                new_nodes[group_nr][i, j, :] = mesh_group.nodes[i, elems, :]
+                elems = queried_elems[j] - num_prev_elems
+                new_nodes[group_nr][i, j - start_idx, :] = mesh_group.nodes[i, elems, :]
 
         index_set |= set(new_indices[group_nr].ravel())
 
+        num_prev_elems += mesh_group.nelements
         start_idx = end_idx
 
     # A sorted np.array of vertex indices we need (without duplicates).
@@ -99,22 +108,24 @@ def partition_mesh(mesh, part_per_element, part_nr):
     for dim in range(mesh.ambient_dim):
         new_vertices[dim] = mesh.vertices[dim][required_indices]
 
-    # Our indices need to be in range [0, len(mesh_group.nelements)].
+    # Our indices need to be in range [0, len(mesh.nelements)].
     for group_nr in range(num_groups):
-        for i in range(len(new_indices[group_nr])):
-            for j in range(len(new_indices[group_nr][0])):
-                original_index = new_indices[group_nr][i, j]
-                new_indices[group_nr][i, j] = np.where(
-                    required_indices == original_index)[0]
+        if not group_nr in skip_groups:
+            for i in range(len(new_indices[group_nr])):
+                for j in range(len(new_indices[group_nr][0])):
+                    original_index = new_indices[group_nr][i, j]
+                    new_indices[group_nr][i, j] = np.where(
+                        required_indices == original_index)[0]
 
     from meshmode.mesh import SimplexElementGroup, Mesh
 
     new_mesh_groups = []
     for group_nr in range(num_groups):
-        mesh_group = mesh.groups[group_nr]
-        new_mesh_groups.append(
-            SimplexElementGroup(mesh_group.order, new_indices[group_nr],
-                new_nodes[group_nr], unit_nodes=mesh_group.unit_nodes))
+        if not group_nr in skip_groups:
+            mesh_group = mesh.groups[group_nr]
+            new_mesh_groups.append(
+                SimplexElementGroup(mesh_group.order, new_indices[group_nr],
+                    new_nodes[group_nr], unit_nodes=mesh_group.unit_nodes))
 
     part_mesh = Mesh(new_vertices, new_mesh_groups)
 
diff --git a/test/test_meshmode.py b/test/test_meshmode.py
index 4cdffd85..30b5e001 100644
--- a/test/test_meshmode.py
+++ b/test/test_meshmode.py
@@ -50,7 +50,7 @@ logger = logging.getLogger(__name__)
 
 # {{{ partition_mesh
 
-@pytest.mark.parametrize("mesh_type", ["torus", "box"])
+@pytest.mark.parametrize("mesh_type", ["torus", "boxes"])
 def test_partition_mesh(mesh_type):
     if mesh_type == "torus":
         from meshmode.mesh.generation import generate_torus
@@ -67,11 +67,13 @@ def test_partition_mesh(mesh_type):
         assert part_mesh1.nelements == 4
         assert part_mesh2.nelements == 2
 
-    elif mesh_type == "box":
-        from meshmode.mesh.generation import generate_box_mesh
-        seg = np.linspace(0, 1, 10)
-        axis_coords = (seg, seg, seg)
-        mesh = generate_box_mesh(axis_coords)
+    elif mesh_type == "boxes":
+        from meshmode.mesh.generation import generate_regular_rect_mesh
+        mesh1 = generate_regular_rect_mesh(a=(0, 0, 0), b=(1, 1, 1), n=(5, 5, 5))
+        mesh2 = generate_regular_rect_mesh(a=(2, 2, 2), b=(3, 3, 3), n=(5, 5, 5))
+
+        from meshmode.mesh.processing import merge_disjoint_meshes
+        mesh = merge_disjoint_meshes([mesh1, mesh2])
 
         adjacency_list = np.zeros((mesh.nelements,), dtype=set)
         for elem in range(mesh.nelements):
@@ -81,12 +83,13 @@ def test_partition_mesh(mesh_type):
                 adjacency_list[elem].add(mesh.nodal_adjacency.neighbors[n])
 
         from pymetis import part_graph
-        (_, part_per_element) = part_graph(3, adjacency=adjacency_list)
+        (_, p) = part_graph(3, adjacency=adjacency_list)
+        part_per_element = np.array(p)
 
         from meshmode.mesh.processing import partition_mesh
-        (part_mesh0, _) = partition_mesh(mesh, np.array(part_per_element), 0)
-        (part_mesh1, _) = partition_mesh(mesh, np.array(part_per_element), 1)
-        (part_mesh2, _) = partition_mesh(mesh, np.array(part_per_element), 2)
+        (part_mesh0, _) = partition_mesh(mesh, part_per_element, 0)
+        (part_mesh1, _) = partition_mesh(mesh, part_per_element, 1)
+        (part_mesh2, _) = partition_mesh(mesh, part_per_element, 2)
 
         assert mesh.nelements == (part_mesh0.nelements
             + part_mesh1.nelements + part_mesh2.nelements)
-- 
GitLab


From ac1f0817c61011db434d9c188517b4cb766a6a0f Mon Sep 17 00:00:00 2001
From: ellis <eshoag2@illinois.edu>
Date: Wed, 15 Feb 2017 14:29:24 -0600
Subject: [PATCH 08/12] pymetis dependency

---
 meshmode/mesh/processing.py | 10 ++++++----
 requirements.txt            |  3 +++
 2 files changed, 9 insertions(+), 4 deletions(-)

diff --git a/meshmode/mesh/processing.py b/meshmode/mesh/processing.py
index 28d146b7..8a6f0464 100644
--- a/meshmode/mesh/processing.py
+++ b/meshmode/mesh/processing.py
@@ -90,11 +90,13 @@ def partition_mesh(mesh, part_per_element, part_nr):
         new_indices.append(mesh_group.vertex_indices[elems])
 
         new_nodes.append(
-            np.zeros((mesh.ambient_dim, end_idx - start_idx, mesh_group.nunit_nodes)))
+            np.zeros(
+                (mesh.ambient_dim, end_idx - start_idx, mesh_group.nunit_nodes)))
         for i in range(mesh.ambient_dim):
             for j in range(start_idx, end_idx):
                 elems = queried_elems[j] - num_prev_elems
-                new_nodes[group_nr][i, j - start_idx, :] = mesh_group.nodes[i, elems, :]
+                new_idx = j - start_idx
+                new_nodes[group_nr][i, new_idx, :] = mesh_group.nodes[i, elems, :]
 
         index_set |= set(new_indices[group_nr].ravel())
 
@@ -110,7 +112,7 @@ def partition_mesh(mesh, part_per_element, part_nr):
 
     # Our indices need to be in range [0, len(mesh.nelements)].
     for group_nr in range(num_groups):
-        if not group_nr in skip_groups:
+        if group_nr not in skip_groups:
             for i in range(len(new_indices[group_nr])):
                 for j in range(len(new_indices[group_nr][0])):
                     original_index = new_indices[group_nr][i, j]
@@ -121,7 +123,7 @@ def partition_mesh(mesh, part_per_element, part_nr):
 
     new_mesh_groups = []
     for group_nr in range(num_groups):
-        if not group_nr in skip_groups:
+        if group_nr not in skip_groups:
             mesh_group = mesh.groups[group_nr]
             new_mesh_groups.append(
                 SimplexElementGroup(mesh_group.order, new_indices[group_nr],
diff --git a/requirements.txt b/requirements.txt
index 6ed102e6..f5bb3464 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -9,3 +9,6 @@ git+https://github.com/inducer/loopy.git
 git+https://github.com/inducer/boxtree.git
 git+https://github.com/inducer/sumpy.git
 git+https://github.com/inducer/pytential.git
+
+# requires pymetis for tests for partition_mesh
+git+https://github.com/inducer/pymetis.git
-- 
GitLab


From a8a93eca3fde16cb43cae5d7f2fb12683c1cd2ec Mon Sep 17 00:00:00 2001
From: ellis <eshoag2@illinois.edu>
Date: Wed, 15 Feb 2017 22:18:15 -0600
Subject: [PATCH 09/12] separated test_partition_mesh

---
 test/test_meshmode.py | 85 +++++++++++++++++++++----------------------
 1 file changed, 42 insertions(+), 43 deletions(-)

diff --git a/test/test_meshmode.py b/test/test_meshmode.py
index 30b5e001..aca8cd4c 100644
--- a/test/test_meshmode.py
+++ b/test/test_meshmode.py
@@ -50,49 +50,48 @@ logger = logging.getLogger(__name__)
 
 # {{{ partition_mesh
 
-@pytest.mark.parametrize("mesh_type", ["torus", "boxes"])
-def test_partition_mesh(mesh_type):
-    if mesh_type == "torus":
-        from meshmode.mesh.generation import generate_torus
-        my_mesh = generate_torus(2, 1, n_outer=2, n_inner=2)
-
-        part_per_element = np.array([0, 1, 2, 1, 1, 2, 1, 0])
-
-        from meshmode.mesh.processing import partition_mesh
-        (part_mesh0, _) = partition_mesh(my_mesh, part_per_element, 0)
-        (part_mesh1, _) = partition_mesh(my_mesh, part_per_element, 1)
-        (part_mesh2, _) = partition_mesh(my_mesh, part_per_element, 2)
-
-        assert part_mesh0.nelements == 2
-        assert part_mesh1.nelements == 4
-        assert part_mesh2.nelements == 2
-
-    elif mesh_type == "boxes":
-        from meshmode.mesh.generation import generate_regular_rect_mesh
-        mesh1 = generate_regular_rect_mesh(a=(0, 0, 0), b=(1, 1, 1), n=(5, 5, 5))
-        mesh2 = generate_regular_rect_mesh(a=(2, 2, 2), b=(3, 3, 3), n=(5, 5, 5))
-
-        from meshmode.mesh.processing import merge_disjoint_meshes
-        mesh = merge_disjoint_meshes([mesh1, mesh2])
-
-        adjacency_list = np.zeros((mesh.nelements,), dtype=set)
-        for elem in range(mesh.nelements):
-            adjacency_list[elem] = set()
-            starts = mesh.nodal_adjacency.neighbors_starts
-            for n in range(starts[elem], starts[elem + 1]):
-                adjacency_list[elem].add(mesh.nodal_adjacency.neighbors[n])
-
-        from pymetis import part_graph
-        (_, p) = part_graph(3, adjacency=adjacency_list)
-        part_per_element = np.array(p)
-
-        from meshmode.mesh.processing import partition_mesh
-        (part_mesh0, _) = partition_mesh(mesh, part_per_element, 0)
-        (part_mesh1, _) = partition_mesh(mesh, part_per_element, 1)
-        (part_mesh2, _) = partition_mesh(mesh, part_per_element, 2)
-
-        assert mesh.nelements == (part_mesh0.nelements
-            + part_mesh1.nelements + part_mesh2.nelements)
+def test_partition_torus_mesh():
+    from meshmode.mesh.generation import generate_torus
+    my_mesh = generate_torus(2, 1, n_outer=2, n_inner=2)
+
+    part_per_element = np.array([0, 1, 2, 1, 1, 2, 1, 0])
+
+    from meshmode.mesh.processing import partition_mesh
+    (part_mesh0, _) = partition_mesh(my_mesh, part_per_element, 0)
+    (part_mesh1, _) = partition_mesh(my_mesh, part_per_element, 1)
+    (part_mesh2, _) = partition_mesh(my_mesh, part_per_element, 2)
+
+    assert part_mesh0.nelements == 2
+    assert part_mesh1.nelements == 4
+    assert part_mesh2.nelements == 2
+
+
+def test_partition_boxes_mesh():
+    from meshmode.mesh.generation import generate_regular_rect_mesh
+    mesh1 = generate_regular_rect_mesh(a=(0, 0, 0), b=(1, 1, 1), n=(5, 5, 5))
+    mesh2 = generate_regular_rect_mesh(a=(2, 2, 2), b=(3, 3, 3), n=(5, 5, 5))
+
+    from meshmode.mesh.processing import merge_disjoint_meshes
+    mesh = merge_disjoint_meshes([mesh1, mesh2])
+
+    adjacency_list = np.zeros((mesh.nelements,), dtype=set)
+    for elem in range(mesh.nelements):
+        adjacency_list[elem] = set()
+        starts = mesh.nodal_adjacency.neighbors_starts
+        for n in range(starts[elem], starts[elem + 1]):
+            adjacency_list[elem].add(mesh.nodal_adjacency.neighbors[n])
+
+    from pymetis import part_graph
+    (_, p) = part_graph(3, adjacency=adjacency_list)
+    part_per_element = np.array(p)
+
+    from meshmode.mesh.processing import partition_mesh
+    (part_mesh0, _) = partition_mesh(mesh, part_per_element, 0)
+    (part_mesh1, _) = partition_mesh(mesh, part_per_element, 1)
+    (part_mesh2, _) = partition_mesh(mesh, part_per_element, 2)
+
+    assert mesh.nelements == (part_mesh0.nelements
+        + part_mesh1.nelements + part_mesh2.nelements)
 
 # }}}
 
-- 
GitLab


From 42e569a244f0e1d028e0082dc80c5710fb537b94 Mon Sep 17 00:00:00 2001
From: ellis <eshoag2@illinois.edu>
Date: Wed, 15 Feb 2017 22:27:59 -0600
Subject: [PATCH 10/12] more general group types

---
 meshmode/mesh/processing.py | 7 ++++++-
 1 file changed, 6 insertions(+), 1 deletion(-)

diff --git a/meshmode/mesh/processing.py b/meshmode/mesh/processing.py
index 8a6f0464..27cb9a03 100644
--- a/meshmode/mesh/processing.py
+++ b/meshmode/mesh/processing.py
@@ -53,6 +53,11 @@ def partition_mesh(mesh, part_per_element, part_nr):
         is a :class:`meshmode.mesh.Mesh` that is a partition of mesh, and
         *part_to_global* is a :class:`numpy.ndarray` mapping element
         numbers on *part_mesh* to ones in *mesh*.
+
+    .. versionadded:: 2017.1
+
+    .. warning:: Interface is not final. Connectivity between elements
+        across groups needs to be added.
     """
     assert len(part_per_element) == mesh.nelements, (
         "part_per_element must have shape (mesh.nelements,)")
@@ -126,7 +131,7 @@ def partition_mesh(mesh, part_per_element, part_nr):
         if group_nr not in skip_groups:
             mesh_group = mesh.groups[group_nr]
             new_mesh_groups.append(
-                SimplexElementGroup(mesh_group.order, new_indices[group_nr],
+                type(mesh_group)(mesh_group.order, new_indices[group_nr],
                     new_nodes[group_nr], unit_nodes=mesh_group.unit_nodes))
 
     part_mesh = Mesh(new_vertices, new_mesh_groups)
-- 
GitLab


From 39643d25d1c5bead0d11a8fc4eac9ba834ae6f9f Mon Sep 17 00:00:00 2001
From: ellis <eshoag2@illinois.edu>
Date: Wed, 15 Feb 2017 23:14:58 -0600
Subject: [PATCH 11/12] variable sizes in test_partition_boxes_mesh

---
 test/test_meshmode.py | 17 +++++++++--------
 1 file changed, 9 insertions(+), 8 deletions(-)

diff --git a/test/test_meshmode.py b/test/test_meshmode.py
index aca8cd4c..a793743c 100644
--- a/test/test_meshmode.py
+++ b/test/test_meshmode.py
@@ -67,9 +67,11 @@ def test_partition_torus_mesh():
 
 
 def test_partition_boxes_mesh():
+    n = 5
+    num_parts = 7
     from meshmode.mesh.generation import generate_regular_rect_mesh
-    mesh1 = generate_regular_rect_mesh(a=(0, 0, 0), b=(1, 1, 1), n=(5, 5, 5))
-    mesh2 = generate_regular_rect_mesh(a=(2, 2, 2), b=(3, 3, 3), n=(5, 5, 5))
+    mesh1 = generate_regular_rect_mesh(a=(0, 0, 0), b=(1, 1, 1), n=(n, n, n))
+    mesh2 = generate_regular_rect_mesh(a=(2, 2, 2), b=(3, 3, 3), n=(n, n, n))
 
     from meshmode.mesh.processing import merge_disjoint_meshes
     mesh = merge_disjoint_meshes([mesh1, mesh2])
@@ -82,16 +84,15 @@ def test_partition_boxes_mesh():
             adjacency_list[elem].add(mesh.nodal_adjacency.neighbors[n])
 
     from pymetis import part_graph
-    (_, p) = part_graph(3, adjacency=adjacency_list)
+    (_, p) = part_graph(num_parts, adjacency=adjacency_list)
     part_per_element = np.array(p)
 
     from meshmode.mesh.processing import partition_mesh
-    (part_mesh0, _) = partition_mesh(mesh, part_per_element, 0)
-    (part_mesh1, _) = partition_mesh(mesh, part_per_element, 1)
-    (part_mesh2, _) = partition_mesh(mesh, part_per_element, 2)
+    new_meshes = [
+        partition_mesh(mesh, part_per_element, i)[0] for i in range(num_parts)]
 
-    assert mesh.nelements == (part_mesh0.nelements
-        + part_mesh1.nelements + part_mesh2.nelements)
+    assert mesh.nelements == np.sum(
+        [new_meshes[i].nelements for i in range(num_parts)])
 
 # }}}
 
-- 
GitLab


From 058fd557d23a01e6ac2ec827518663e6efb7f3ce Mon Sep 17 00:00:00 2001
From: ellis <eshoag2@illinois.edu>
Date: Thu, 16 Feb 2017 10:21:48 -0600
Subject: [PATCH 12/12] possible optimizations for required_indices

---
 meshmode/mesh/processing.py | 15 ++++++++++-----
 1 file changed, 10 insertions(+), 5 deletions(-)

diff --git a/meshmode/mesh/processing.py b/meshmode/mesh/processing.py
index 27cb9a03..37e4ac26 100644
--- a/meshmode/mesh/processing.py
+++ b/meshmode/mesh/processing.py
@@ -70,7 +70,11 @@ def partition_mesh(mesh, part_per_element, part_nr):
     new_nodes = []
 
     # The set of vertex indices we need.
-    index_set = set()
+    # NOTE: There are two methods for producing required_indices.
+    #   Optimizations may come from further exploring these options.
+    #index_set = np.array([], dtype=int)
+    index_sets = np.array([], dtype=set)
+
     skip_groups = []
     num_prev_elems = 0
     start_idx = 0
@@ -103,13 +107,15 @@ def partition_mesh(mesh, part_per_element, part_nr):
                 new_idx = j - start_idx
                 new_nodes[group_nr][i, new_idx, :] = mesh_group.nodes[i, elems, :]
 
-        index_set |= set(new_indices[group_nr].ravel())
+        #index_set = np.append(index_set, new_indices[group_nr].ravel())
+        index_sets = np.append(index_sets, set(new_indices[group_nr].ravel()))
 
         num_prev_elems += mesh_group.nelements
         start_idx = end_idx
 
     # A sorted np.array of vertex indices we need (without duplicates).
-    required_indices = np.array(list(index_set))
+    #required_indices = np.unique(np.sort(index_set))
+    required_indices = np.array(list(set.union(*index_sets)))
 
     new_vertices = np.zeros((mesh.ambient_dim, len(required_indices)))
     for dim in range(mesh.ambient_dim):
@@ -124,8 +130,6 @@ def partition_mesh(mesh, part_per_element, part_nr):
                     new_indices[group_nr][i, j] = np.where(
                         required_indices == original_index)[0]
 
-    from meshmode.mesh import SimplexElementGroup, Mesh
-
     new_mesh_groups = []
     for group_nr in range(num_groups):
         if group_nr not in skip_groups:
@@ -134,6 +138,7 @@ def partition_mesh(mesh, part_per_element, part_nr):
                 type(mesh_group)(mesh_group.order, new_indices[group_nr],
                     new_nodes[group_nr], unit_nodes=mesh_group.unit_nodes))
 
+    from meshmode.mesh import Mesh
     part_mesh = Mesh(new_vertices, new_mesh_groups)
 
     return (part_mesh, queried_elems)
-- 
GitLab