Add random Voronoi mesh to mesh generators

src/mesh/python/mesh_types.py

@@ -15,6 +15,7 @@ class base_mesh:
     Base class with data required by mesh file generators.
     This contains no methods for creating the data members.
     '''
+
     def __init__(self):
         # -- required data
         self.ndim = 0  # number of dimensions
@@ -24,9 +25,9 @@ def __init__(self):
         self.num_faces = 0  # total number of oriented faces
         self.num_faces_per_cell = np.array([], dtype=int)  # number of faces per cell
         self.num_nodes_per_face = np.array([], dtype=int)  # number of nodes per face
-        self.faces_per_cell = np.array([[]], dtype=int)  # face indexes per cell
-        self.nodes_per_face = np.array([[]], dtype=int)  # node indexes per face
-        self.nodes_per_side = [np.array([[]], dtype=int)]  # list of arrays of node per bdy face
+        self.faces_per_cell = [np.array([], dtype=int)]  # face indexes per cell
+        self.nodes_per_face = [np.array([], dtype=int)]  # node indexes per face
+        self.nodes_per_side = [[np.array([], dtype=int)]]  # list of arrays of node per bdy face
 
 
 # ------------------------------------------------------------------------------------------------ #
@@ -37,6 +38,7 @@ class orth_2d_mesh(base_mesh):
     This class generates an orthogonally structured mesh in an unstructured
     format suitable for creating unstructured-mesh input files.
     '''
+
     def __init__(self, bounds_per_dim, num_cells_per_dim):
 
         # -- short-cuts
@@ -125,6 +127,7 @@ class orth_3d_mesh(base_mesh):
     This class generates an orthogonally structured mesh in an unstructured
     format suitable for creating unstructured-mesh input files.
     '''
+
     def __init__(self, bounds_per_dim, num_cells_per_dim):
 
         # -- short-cuts
@@ -287,6 +290,319 @@ def __init__(self, bounds_per_dim, num_cells_per_dim):
                                nodes_per_side_zlow, nodes_per_side_zhig]
 
 
+# ------------------------------------------------------------------------------------------------ #
+# Voronoi 2D mesh type
+class vor_2d_mesh(base_mesh):
+    '''
+    Class for randomly generated 2D Voronoi mesh data.
+    This class generates a set of random nodes inside a bounding box, and then
+    constructs a mesh from the Voronoi diagram of these points. It casts this
+    mesh in an unstructured format suitable for creating unstructured-mesh input
+    files.
+    '''
+
+    def __init__(self, bounds_per_dim, num_cells, seed=0):
+        def soft_equiv(a, b, eps=1.e-12):
+            if 2. * np.fabs(a - b) / (a + b + eps) < eps:
+                return True
+            else:
+                return False
+        from scipy.spatial import Voronoi
+        np.random.seed(seed)
+
+        # -- convenience variables
+        xmin = bounds_per_dim[0][0]
+        xmax = bounds_per_dim[0][1]
+        ymin = bounds_per_dim[1][0]
+        ymax = bounds_per_dim[1][1]
+
+        # -- number of dimensions
+        ndim = len(bounds_per_dim)
+        self.ndim = ndim
+        self.num_cells = num_cells
+        assert (ndim == 2), 'ndim != 2, exiting...'
+
+        # -- randomly distribute grid-generating points inside bounding box
+        points = np.zeros([num_cells, 2])
+        for n in range(num_cells):
+            points[n, 0] = xmin + (xmax - xmin) * np.random.uniform()
+            points[n, 1] = ymin + (ymax - ymin) * np.random.uniform()
+
+        # -- create base Voronoi diagram
+        vor = Voronoi(points)
+
+        # -- sanitize Voronoi diagram (remove empty regions, remove vertices and
+        #    edges outside of domain
+        points = vor.points
+        vertices = []
+        ridge_points = []
+        ridge_vertices = []
+
+        new_to_old_vertex_indices = []
+        old_to_new_vertex_indices = []
+        for n, vertex in enumerate(vor.vertices):
+            x = vertex[0]
+            y = vertex[1]
+            if not (x < xmin or x > xmax or y < ymin or y > ymax):
+                vertices.append(vertex.tolist())
+                new_to_old_vertex_indices.append(n)
+            if n in new_to_old_vertex_indices:
+                old_to_new_vertex_indices.append(new_to_old_vertex_indices.index(n))
+            else:
+                old_to_new_vertex_indices.append(-1)
+
+        # -- update vertex indices for ridge_vertices
+        ridge_vertices = vor.ridge_vertices
+        for n in range(len(ridge_vertices)):
+            for m in range(2):
+                if ridge_vertices[n][m] != -1:
+                    ridge_vertices[n][m] = old_to_new_vertex_indices[ridge_vertices[n][m]]
+
+        # -- remove ridges from ridge_vertices and ridge_points if [-1,-1] i.e.
+        #    entirely outside of bounding box
+        ridge_points = vor.ridge_points.tolist()
+        tmp_ridge_vertices = []
+        tmp_ridge_points = []
+        for n in range(len(ridge_vertices)):
+            if not (ridge_vertices[n][0] == -1 and ridge_vertices[n][1] == -1):
+                tmp_ridge_vertices.append(ridge_vertices[n])
+                tmp_ridge_points.append(ridge_points[n])
+        ridge_vertices = tmp_ridge_vertices
+        ridge_points = tmp_ridge_points
+
+        # -- helper function for logic identifying which direction to find
+        #    boundary intersections in
+        def ray_segment_intersection(origin, direction, pt1, pt2):
+            def mag(v):
+                return np.sqrt(np.dot(np.array(v), np.array(v)))
+
+            def norm(v):
+                return np.array(v) / mag(v)
+            origin = np.array(origin)
+            direction = np.array(norm(direction))
+            pt1 = np.array(pt1)
+            pt2 = np.array(pt2)
+
+            v1 = origin - pt1
+            v2 = pt2 - pt1
+            v3 = np.array([-direction[1], direction[0]])
+            t1 = np.cross(v2, v1) / np.dot(v2, v3)
+            t2 = np.dot(v1, v3) / np.dot(v2, v3)
+            if t1 >= 0. and t2 >= 0. and t2 <= 1.:
+                return origin + t1 * direction
+            else:
+                return None
+
+        # -- add boundary vertices for ridges with a -1 (i.e. missing) vertex,
+        #    with some rather odd logic for ensuring that we cast rays from
+        #    existing vertices towards the correct boundary
+        new_vertices_start_index = len(vertices)
+        for n, ridge in enumerate(ridge_vertices):
+            bad_vertex = None
+            if ridge[0] == -1:
+                good_vertex = 1
+                bad_vertex = 0
+            if ridge[1] == -1:
+                good_vertex = 0
+                bad_vertex = 1
+
+            if bad_vertex is not None:
+                # -- determine which cells own this ridge
+                point1 = points[ridge_points[n][0]]
+                point2 = points[ridge_points[n][1]]
+                midpoint = [(point1[0] + point2[0]) / 2, (point1[1] + point2[1]) / 2]
+                distances = np.zeros(4)
+                bpt1s = [[xmin, xmin], [xmax, xmax], [xmin, xmax], [xmin, xmax]]
+                bpt2s = [[ymin, ymax], [ymin, ymax], [ymin, ymin], [ymax, ymax]]
+
+                origin = [vertices[ridge[good_vertex]][0], vertices[ridge[good_vertex]][1]]
+                direction = [midpoint[0] - vertices[ridge[good_vertex]][0],
+                             midpoint[1] - vertices[ridge[good_vertex]][1]]
+                found = False
+                for bc in range(4):
+                    pt1 = [bpt1s[bc][0], bpt2s[bc][0]]
+                    pt2 = [bpt1s[bc][1], bpt2s[bc][1]]
+                    # -- check reverse directions only if forward directions don't work
+                    for d in range(1):
+                        if d == 1:
+                            direction[0] = -direction[0]
+                            direction[1] = -direction[1]
+                        newpt = ray_segment_intersection(origin, direction, pt1, pt2)
+                        if newpt is not None:
+                            # -- check for intersection with other faces
+                            bad = False
+                            for tmpridge in ridge_vertices:
+                                if (tmpridge[0] == len(vertices) - 1 or
+                                    tmpridge[1] == len(vertices) - 1 or
+                                    tmpridge[0] == ridge[good_vertex] or
+                                    tmpridge[1] == ridge[good_vertex] or
+                                    tmpridge[0] == -1 or tmpridge[1] == -1):
+                                    continue
+                                tstpt = ray_segment_intersection(origin, direction,
+                                                                 [vertices[tmpridge[0]][0],
+                                                                  vertices[tmpridge[0]][1]],
+                                                                 [vertices[tmpridge[1]][0],
+                                                                     vertices[tmpridge[1]][1]])
+                                if tstpt is not None:
+                                    bad = True
+                                    break
+                            if not bad:
+                                vertices.append([newpt[0], newpt[1]])
+                                ridge_vertices[n][bad_vertex] = len(vertices) - 1
+                                found = True
+                if found is False:
+                    direction[0] = -direction[0]
+                    direction[1] = -direction[1]
+                    for bc in range(4):
+                        pt1 = [bpt1s[bc][0], bpt2s[bc][0]]
+                        pt2 = [bpt1s[bc][1], bpt2s[bc][1]]
+                        for d in range(1, 2):
+                            newpt = ray_segment_intersection(origin, direction, pt1, pt2)
+                            if newpt is not None:
+                                # -- check for intersection with other faces
+                                bad = False
+                                for tmpridge in ridge_vertices:
+                                    if (tmpridge[0] == len(vertices) - 1 or
+                                        tmpridge[1] == len(vertices) - 1 or
+                                        tmpridge[0] == ridge[good_vertex] or
+                                        tmpridge[1] == ridge[good_vertex] or
+                                        tmpridge[0] == -1 or tmpridge[1] == -1):
+                                        continue
+                                    tstpt = ray_segment_intersection(origin, direction,
+                                                                     [vertices[tmpridge[0]][0],
+                                                                      vertices[tmpridge[0]][1]],
+                                                                     [vertices[tmpridge[1]][0],
+                                                                         vertices[tmpridge[1]][1]])
+                                    if tstpt is not None:
+                                        bad = True
+                                        break
+                                if not bad:
+                                    vertices.append([newpt[0], newpt[1]])
+                                    ridge_vertices[n][bad_vertex] = len(vertices) - 1
+                                    found = True
+                    assert (found is not False), 'No suitable boundary point found!'
+                continue
+
+        # -- add boundary edges, ignoring corner vertices for now
+        points_to_work_on = []
+        points_to_work_on_vertices = []
+        for n, ridge in enumerate(ridge_vertices):
+            v1 = ridge[0]
+            v2 = ridge[1]
+            if v1 >= new_vertices_start_index or v2 >= new_vertices_start_index:
+                if v1 >= new_vertices_start_index:
+                    new = v1
+                else:
+                    new = v2
+                for point in ridge_points[n]:
+                    if point not in points_to_work_on:
+                        points_to_work_on.append(point)
+                        points_to_work_on_vertices.append([new, -1])
+                    else:
+                        idx = points_to_work_on.index(point)
+                        points_to_work_on_vertices[idx][1] = new
+
+        # -- add corners
+        for verts in points_to_work_on_vertices:
+            v1 = vertices[verts[0]]
+            v2 = vertices[verts[1]]
+            # -- bottom left
+            if (soft_equiv(v1[0], xmin) or soft_equiv(v2[0], xmin)) and \
+               (soft_equiv(v1[1], ymin) or soft_equiv(v2[1], ymin)):
+                vertices.append([xmin, ymin])
+                ridge_vertices.append([verts[0], len(vertices) - 1])
+                ridge_vertices.append([verts[1], len(vertices) - 1])
+            # -- top left
+            elif (soft_equiv(v1[0], xmin) or soft_equiv(v2[0], xmin)) and \
+                 (soft_equiv(v1[1], ymax) or soft_equiv(v2[1], ymax)):
+                vertices.append([xmin, ymax])
+                ridge_vertices.append([verts[0], len(vertices) - 1])
+                ridge_vertices.append([verts[1], len(vertices) - 1])
+            # -- bottom right
+            elif (soft_equiv(v1[0], xmax) or soft_equiv(v2[0], xmax)) and \
+                 (soft_equiv(v1[1], ymin) or soft_equiv(v2[1], ymin)):
+                vertices.append([xmax, ymin])
+                ridge_vertices.append([verts[0], len(vertices) - 1])
+                ridge_vertices.append([verts[1], len(vertices) - 1])
+            # -- top right
+            elif (soft_equiv(v1[0], xmax) or soft_equiv(v2[0], xmax)) and \
+                 (soft_equiv(v1[1], ymax) or soft_equiv(v2[1], ymax)):
+                vertices.append([xmax, ymax])
+                ridge_vertices.append([verts[0], len(vertices) - 1])
+                ridge_vertices.append([verts[1], len(vertices) - 1])
+            # -- not a corner
+            else:
+                ridge_vertices.append(verts)
+
+        # -- assign ridge vertices to regions and boundaries
+        cells = []
+        cell_nodes = []
+        boundary_edges = {}
+        boundary_edges['xl'] = []
+        boundary_edges['xr'] = []
+        boundary_edges['yl'] = []
+        boundary_edges['yr'] = []
+        for n, point in enumerate(points):
+            cell_nodes.append(n)
+            cells.append([])
+        for ridge_idx, v_indices in enumerate(ridge_vertices):
+            midpoint = [(vertices[v_indices[0]][0] + vertices[v_indices[1]][0]) / 2,
+                        (vertices[v_indices[0]][1] + vertices[v_indices[1]][1]) / 2]
+            distances = np.zeros(len(cell_nodes))
+            for n, node in enumerate(cell_nodes):
+                distances[n] = np.sqrt((points[n][0] - midpoint[0])**2 +
+                                       (points[n][1] - midpoint[1])**2)
+            indices = np.argsort(distances)
+            distances = np.sort(distances)
+            if soft_equiv(distances[0], distances[1]):
+                # -- not a boundary
+                cells[indices[0]].append(ridge_idx)
+                cells[indices[1]].append(ridge_idx)
+            else:
+                # -- a boundary
+                cells[indices[0]].append(ridge_idx)
+                if (soft_equiv(vertices[v_indices[0]][0], xmin) and
+                    soft_equiv(vertices[v_indices[1]][0], xmin)):
+                    boundary_edges['xl'].append(ridge_idx)
+                elif (soft_equiv(vertices[v_indices[0]][0], xmax) and
+                      soft_equiv(vertices[v_indices[1]][0], xmax)):
+                    boundary_edges['xr'].append(ridge_idx)
+                elif (soft_equiv(vertices[v_indices[0]][1], ymin) and
+                      soft_equiv(vertices[v_indices[1]][1], ymin)):
+                    boundary_edges['yl'].append(ridge_idx)
+                elif (soft_equiv(vertices[v_indices[0]][1], ymax) and
+                      soft_equiv(vertices[v_indices[1]][1], ymax)):
+                    boundary_edges['yr'].append(ridge_idx)
+                else:
+                    assert (False), 'Boundary edge not identified'
+
+        # -- update remaining base values
+        self.num_nodes = len(vertices)
+        self.coordinates_per_node = np.zeros([self.num_nodes, 2])
+        for n, vertex in enumerate(vertices):
+            self.coordinates_per_node[n, 0] = vertex[0]
+            self.coordinates_per_node[n, 1] = vertex[1]
+        self.num_faces = len(ridge_vertices)
+        self.num_faces_per_cell = np.zeros(self.num_cells, dtype=int)
+        for n in range(self.num_cells):
+            self.num_faces_per_cell[n] = len(cells[n])
+        self.num_nodes_per_face = np.zeros(self.num_faces, dtype=int)
+        for n in range(self.num_faces):
+            self.num_nodes_per_face[n] = 2
+        self.faces_per_cell = cells
+        self.nodes_per_face = ridge_vertices
+        self.nodes_per_side = []
+        for n in range(4):
+            bdy_key = list(boundary_edges.keys())[n]
+            bdy_nodes = []
+            for bdy in boundary_edges[bdy_key]:
+                nodes = ridge_vertices[bdy]
+                for node in nodes:
+                    if node not in bdy_nodes:
+                        bdy_nodes.append(node)
+            self.nodes_per_side.append(bdy_nodes)
+
+
 # ------------------------------------------------------------------------------------------------ #
 # end of mesh_types.py
 # ------------------------------------------------------------------------------------------------ #