From 7a924cd12409226aef3dcf57d94fc3da2bce9773 Mon Sep 17 00:00:00 2001 From: "Ben R. Ryan" Date: Thu, 5 Aug 2021 10:43:38 -0600 Subject: [PATCH 1/8] First pass --- src/mesh/python/mesh_types.py | 307 ++++++++++++++++++++++++++++++++++ 1 file changed, 307 insertions(+) diff --git a/src/mesh/python/mesh_types.py b/src/mesh/python/mesh_types.py index 28ce63b66..302facec5 100644 --- a/src/mesh/python/mesh_types.py +++ b/src/mesh/python/mesh_types.py @@ -287,6 +287,313 @@ 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 soft_equiv(a, b, eps = 1.e-12): + if 2.*np.fabs(a - b)/(a + b + eps) < eps: + return True + else: + return False + + def __init__(self, bounds_per_dim, num_cells, seed=0): + 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 + 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 + for region in vor.regions: + if len(region) > 0: + regions.append(region) + + 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 = deepcopy(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: + # Which cells own this ridge? + for region in regions: + if ridge[good_vertex] in region: + + 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 == 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 + if found == False: + print("ERROR No suitable boundary point found!") + sys.exit() + 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: + orig = v2 + new = v1 + else: + orig = v1 + 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: + print("ERROR Boundary edge not identified!") + sys.exit() + + # -- update remaining base values + self.num_nodes = len(vertices) + self.coordinates_per_node = vertices + self.num_faces = len(ridge_vertices) + self.num_faces_per_cell = np.zeros(self.num_cells) + for n in range(self.num_cells): + self.num_faces_per_cell = len(cells[n]) + self.num_nodes_per_face = np.zeros(self.num_faces) + 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 boudary_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 # ------------------------------------------------------------------------------------------------ # From b3cf5d52700208978bb25deb3f0510b1dd8b47c8 Mon Sep 17 00:00:00 2001 From: "Ben R. Ryan" Date: Thu, 5 Aug 2021 10:53:10 -0600 Subject: [PATCH 2/8] Bug --- src/mesh/python/mesh_types.py | 27 ++++++++++++--------------- src/mesh/python/x3d_generator.py | 10 ++++++++-- 2 files changed, 20 insertions(+), 17 deletions(-) diff --git a/src/mesh/python/mesh_types.py b/src/mesh/python/mesh_types.py index 302facec5..c6bfa08a0 100644 --- a/src/mesh/python/mesh_types.py +++ b/src/mesh/python/mesh_types.py @@ -402,10 +402,7 @@ def norm(v): bad_vertex = 1 if bad_vertex is not None: - # Which cells own this ridge? - for region in regions: - if ridge[good_vertex] in region: - + # -- 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] @@ -429,17 +426,17 @@ def norm(v): 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 + 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 diff --git a/src/mesh/python/x3d_generator.py b/src/mesh/python/x3d_generator.py index 98e82269c..4c6af2226 100755 --- a/src/mesh/python/x3d_generator.py +++ b/src/mesh/python/x3d_generator.py @@ -10,7 +10,8 @@ import argparse # -- mesh class dictionary -mesh_type_dict = {'orth_2d_mesh': mesh_types.orth_2d_mesh, 'orth_3d_mesh': mesh_types.orth_3d_mesh} +mesh_type_dict = {'orth_2d_mesh': mesh_types.orth_2d_mesh, 'orth_3d_mesh': mesh_types.orth_3d_mesh, + 'vor_2d_mesh': mesh_types.vor_2d_mesh} # ------------------------------------------------------------------------------------------------ # # -- create argument parser @@ -24,6 +25,8 @@ help='Length per dimension.') parser.add_argument('--name', type=str, default='mesh', help='Select file name (will be prefixed with x3d. and sufficed with .in).') +parser.add_argument('--num_cells', type=int, nargs=1, default=100, + help='Number of cells') # -- parse arguments from command line args = parser.parse_args() @@ -43,7 +46,10 @@ bnd_per_dim = [[args.bnd_per_dim[2 * i], args.bnd_per_dim[2 * i + 1]] for i in range(ndim)] # -- instantiate the class for the mesh type selected -mesh = mesh_type_dict[args.mesh_type](bnd_per_dim, args.num_per_dim) +if args.mesh_type in ['orth_2d_mesh', 'orth_3d_mesh']: + mesh = mesh_type_dict[args.mesh_type](bnd_per_dim, args.num_per_dim) +elif args.mesh_type in ['vor_2d_mesh']: + mesh = mesh_type_dict[args.mesh_type](bnd_per_dim, args.num_cells) # ------------------------------------------------------------------------------------------------ # # -- write out the main mesh file in x3d format From e814cb5ba01b4703b7eea532471088b9345fec19 Mon Sep 17 00:00:00 2001 From: "Ben R. Ryan" Date: Thu, 5 Aug 2021 11:18:19 -0600 Subject: [PATCH 3/8] Runs but data in bad format --- src/mesh/python/mesh_types.py | 42 +++++++++++++++++--------------- src/mesh/python/x3d_generator.py | 2 +- 2 files changed, 24 insertions(+), 20 deletions(-) diff --git a/src/mesh/python/mesh_types.py b/src/mesh/python/mesh_types.py index c6bfa08a0..978d7fdf2 100644 --- a/src/mesh/python/mesh_types.py +++ b/src/mesh/python/mesh_types.py @@ -297,13 +297,13 @@ class vor_2d_mesh(base_mesh): mesh in an unstructured format suitable for creating unstructured-mesh input files. ''' - def soft_equiv(a, b, eps = 1.e-12): - if 2.*np.fabs(a - b)/(a + b + eps) < eps: - return True - else: - return False 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) @@ -316,6 +316,7 @@ def __init__(self, bounds_per_dim, num_cells, seed=0): # -- 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 @@ -330,9 +331,12 @@ def __init__(self, bounds_per_dim, num_cells, seed=0): # -- sanitize Voronoi diagram (remove empty regions, remove vertices and # edges outside of domain points = vor.points - for region in vor.regions: - if len(region) > 0: - regions.append(region) + vertices = [] + ridge_points = [] + ridge_vertices = [] + #for region in vor.regions: + # if len(region) > 0: + # regions.append(region) new_to_old_vertex_indices = [] old_to_new_vertex_indices = [] @@ -348,7 +352,7 @@ def __init__(self, bounds_per_dim, num_cells, seed=0): old_to_new_vertex_indices.append(-1) # -- update vertex indices for ridge_vertices - ridge_vertices = deepcopy(vor.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: @@ -573,22 +577,22 @@ def norm(v): self.num_faces = len(ridge_vertices) self.num_faces_per_cell = np.zeros(self.num_cells) for n in range(self.num_cells): - self.num_faces_per_cell = len(cells[n]) + self.num_faces_per_cell = len(cells[n]) self.num_nodes_per_face = np.zeros(self.num_faces) for n in range(self.num_faces): - self.num_nodes_per_face[n] = 2 + 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 boudary_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) + 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) # ------------------------------------------------------------------------------------------------ # diff --git a/src/mesh/python/x3d_generator.py b/src/mesh/python/x3d_generator.py index 4c6af2226..c19452fe7 100755 --- a/src/mesh/python/x3d_generator.py +++ b/src/mesh/python/x3d_generator.py @@ -25,7 +25,7 @@ help='Length per dimension.') parser.add_argument('--name', type=str, default='mesh', help='Select file name (will be prefixed with x3d. and sufficed with .in).') -parser.add_argument('--num_cells', type=int, nargs=1, default=100, +parser.add_argument('--num_cells', type=int, default=100, help='Number of cells') # -- parse arguments from command line From 5209b13d62e3758ccef4deaac6ea7ae5fadf4309 Mon Sep 17 00:00:00 2001 From: "Ben R. Ryan" Date: Thu, 5 Aug 2021 14:21:40 -0600 Subject: [PATCH 4/8] Working --- src/mesh/python/mesh_types.py | 17 +++++++++++------ src/mesh/python/x3d_generator.py | 10 ++++++---- 2 files changed, 17 insertions(+), 10 deletions(-) diff --git a/src/mesh/python/mesh_types.py b/src/mesh/python/mesh_types.py index 978d7fdf2..bf34ec195 100644 --- a/src/mesh/python/mesh_types.py +++ b/src/mesh/python/mesh_types.py @@ -19,7 +19,7 @@ def __init__(self): # -- required data self.ndim = 0 # number of dimensions self.num_nodes = 0 # total number of nodes - self.coordinates_per_node = np.array([]) # coordinate array indexed by node + self.noordinates_per_node = np.array([]) # coordinate array indexed by node self.num_cells = 0 # total number of cells self.num_faces = 0 # total number of oriented faces self.num_faces_per_cell = np.array([], dtype=int) # number of faces per cell @@ -573,12 +573,17 @@ def norm(v): # -- update remaining base values self.num_nodes = len(vertices) - self.coordinates_per_node = 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] + + #vertices self.num_faces = len(ridge_vertices) - self.num_faces_per_cell = np.zeros(self.num_cells) + self.num_faces_per_cell = np.zeros(self.num_cells, dtype=int) for n in range(self.num_cells): - self.num_faces_per_cell = len(cells[n]) - self.num_nodes_per_face = np.zeros(self.num_faces) + 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 @@ -592,7 +597,7 @@ def norm(v): for node in nodes: if node not in bdy_nodes: bdy_nodes.append(node) - self.nodes_per_side.append(bdy_nodes) + self.nodes_per_side.append(bdy_nodes) # ------------------------------------------------------------------------------------------------ # diff --git a/src/mesh/python/x3d_generator.py b/src/mesh/python/x3d_generator.py index c19452fe7..4534c889f 100755 --- a/src/mesh/python/x3d_generator.py +++ b/src/mesh/python/x3d_generator.py @@ -111,6 +111,7 @@ for node in range(mesh.num_nodes): crds = [0.0, 0.0, 0.0] for idim in range(mesh.ndim): + print(node, idim, mesh.coordinates_per_node[node, idim]) crds[idim] = mesh.coordinates_per_node[node, idim] fo.write('{0:10d} {1:.15e} {2:.15e} {3:.15e}\n'.format(node + 1, crds[0], crds[1], crds[2])) fo.write('end_nodes\n') @@ -123,13 +124,14 @@ for cell in range(mesh.num_cells): for j in range(mesh.num_faces_per_cell[cell]): # -- get the full face index for the cell and local face index - face = mesh.faces_per_cell[cell, j] + #face = mesh.faces_per_cell[cell, j] + face = mesh.faces_per_cell[cell][j] # -- get the number of nodes for this face nnpf = mesh.num_nodes_per_face[face] # -- build face string face_str = '{0:10d}{1:10d}'.format(face + 1, nnpf) # -- append node indices to face string - for node in mesh.nodes_per_face[face, :]: + for node in mesh.nodes_per_face[face]:#, :]: face_str += '{0:10d}'.format(node + 1) # -- write face line fo.write(face_str + '\n') @@ -146,7 +148,7 @@ # -- build cell string cell_str = '{0:10d}{1:10d}'.format(cell + 1, nfpc) # -- append face indices to cell string - for face in mesh.faces_per_cell[cell, :]: + for face in mesh.faces_per_cell[cell]:#, :]: cell_str += '{0:10d}'.format(face + 1) # -- write cell line fo.write(cell_str + '\n') @@ -184,7 +186,7 @@ # -- assume two boundaries per dimension # -- x3d only needs a unique node list for i in range(2 * mesh.ndim): - np.savetxt(fname + '.bdy' + str(i + 1) + '.in', np.unique(mesh.nodes_per_side[i] + 1), fmt='%d') + np.savetxt(fname + '.bdy' + str(i + 1) + '.in', np.unique(np.array(mesh.nodes_per_side[i]) + 1), fmt='%d') # ------------------------------------------------------------------------------------------------ # # end of x3d_generator.py From 97086ae4d589eb86497190182261583847a021c2 Mon Sep 17 00:00:00 2001 From: "Ben R. Ryan" Date: Thu, 5 Aug 2021 14:43:01 -0600 Subject: [PATCH 5/8] Cleanup --- src/mesh/python/mesh_types.py | 5 ----- src/mesh/python/x3d_generator.py | 9 ++++----- 2 files changed, 4 insertions(+), 10 deletions(-) diff --git a/src/mesh/python/mesh_types.py b/src/mesh/python/mesh_types.py index bf34ec195..f425e4c37 100644 --- a/src/mesh/python/mesh_types.py +++ b/src/mesh/python/mesh_types.py @@ -334,9 +334,6 @@ def soft_equiv(a, b, eps = 1.e-12): vertices = [] ridge_points = [] ridge_vertices = [] - #for region in vor.regions: - # if len(region) > 0: - # regions.append(region) new_to_old_vertex_indices = [] old_to_new_vertex_indices = [] @@ -577,8 +574,6 @@ def norm(v): for n, vertex in enumerate(vertices): self.coordinates_per_node[n, 0] = vertex[0] self.coordinates_per_node[n, 1] = vertex[1] - - #vertices 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): diff --git a/src/mesh/python/x3d_generator.py b/src/mesh/python/x3d_generator.py index 4534c889f..fb8869eba 100755 --- a/src/mesh/python/x3d_generator.py +++ b/src/mesh/python/x3d_generator.py @@ -111,7 +111,6 @@ for node in range(mesh.num_nodes): crds = [0.0, 0.0, 0.0] for idim in range(mesh.ndim): - print(node, idim, mesh.coordinates_per_node[node, idim]) crds[idim] = mesh.coordinates_per_node[node, idim] fo.write('{0:10d} {1:.15e} {2:.15e} {3:.15e}\n'.format(node + 1, crds[0], crds[1], crds[2])) fo.write('end_nodes\n') @@ -124,14 +123,13 @@ for cell in range(mesh.num_cells): for j in range(mesh.num_faces_per_cell[cell]): # -- get the full face index for the cell and local face index - #face = mesh.faces_per_cell[cell, j] face = mesh.faces_per_cell[cell][j] # -- get the number of nodes for this face nnpf = mesh.num_nodes_per_face[face] # -- build face string face_str = '{0:10d}{1:10d}'.format(face + 1, nnpf) # -- append node indices to face string - for node in mesh.nodes_per_face[face]:#, :]: + for node in mesh.nodes_per_face[face]: face_str += '{0:10d}'.format(node + 1) # -- write face line fo.write(face_str + '\n') @@ -148,7 +146,7 @@ # -- build cell string cell_str = '{0:10d}{1:10d}'.format(cell + 1, nfpc) # -- append face indices to cell string - for face in mesh.faces_per_cell[cell]:#, :]: + for face in mesh.faces_per_cell[cell]: cell_str += '{0:10d}'.format(face + 1) # -- write cell line fo.write(cell_str + '\n') @@ -186,7 +184,8 @@ # -- assume two boundaries per dimension # -- x3d only needs a unique node list for i in range(2 * mesh.ndim): - np.savetxt(fname + '.bdy' + str(i + 1) + '.in', np.unique(np.array(mesh.nodes_per_side[i]) + 1), fmt='%d') + np.savetxt(fname + '.bdy' + str(i + 1) + '.in', + np.unique(np.array(mesh.nodes_per_side[i]) + 1), fmt='%d') # ------------------------------------------------------------------------------------------------ # # end of x3d_generator.py From c74a8d318b12928160165772283aabff885539bf Mon Sep 17 00:00:00 2001 From: "Ben R. Ryan" Date: Thu, 5 Aug 2021 16:04:03 -0600 Subject: [PATCH 6/8] about to autopep8 --- src/mesh/python/mesh_types.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/mesh/python/mesh_types.py b/src/mesh/python/mesh_types.py index f425e4c37..fcf5946f2 100644 --- a/src/mesh/python/mesh_types.py +++ b/src/mesh/python/mesh_types.py @@ -19,7 +19,7 @@ def __init__(self): # -- required data self.ndim = 0 # number of dimensions self.num_nodes = 0 # total number of nodes - self.noordinates_per_node = np.array([]) # coordinate array indexed by node + self.coordinates_per_node = np.array([]) # coordinate array indexed by node self.num_cells = 0 # total number of cells self.num_faces = 0 # total number of oriented faces self.num_faces_per_cell = np.array([], dtype=int) # number of faces per cell From ff1d4860e7750b5b95d61396ba30d3340844c0d2 Mon Sep 17 00:00:00 2001 From: "Ben R. Ryan" Date: Thu, 5 Aug 2021 16:09:30 -0600 Subject: [PATCH 7/8] Formatted --- src/mesh/python/mesh_types.py | 146 ++++++++++++++++++---------------- 1 file changed, 77 insertions(+), 69 deletions(-) diff --git a/src/mesh/python/mesh_types.py b/src/mesh/python/mesh_types.py index fcf5946f2..1cef6626b 100644 --- a/src/mesh/python/mesh_types.py +++ b/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 @@ -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 @@ -299,8 +302,8 @@ class vor_2d_mesh(base_mesh): ''' 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: + def soft_equiv(a, b, eps=1.e-12): + if 2. * np.fabs(a - b) / (a + b + eps) < eps: return True else: return False @@ -322,8 +325,8 @@ def soft_equiv(a, b, eps = 1.e-12): # -- 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() + 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) @@ -371,9 +374,10 @@ def soft_equiv(a, b, eps = 1.e-12): # 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))) + return np.sqrt(np.dot(np.array(v), np.array(v))) + def norm(v): - return np.array(v)/mag(v) + return np.array(v) / mag(v) origin = np.array(origin) direction = np.array(norm(direction)) pt1 = np.array(pt1) @@ -382,10 +386,10 @@ def norm(v): 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) + 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 + return origin + t1 * direction else: return None @@ -406,10 +410,10 @@ def norm(v): # -- 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] + 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]] + 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], @@ -428,49 +432,55 @@ def norm(v): # -- 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: + 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]]) + [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 + ridge_vertices[n][bad_vertex] = len(vertices) - 1 found = True - if found == False: + 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): + 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: + 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]]) + [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 + ridge_vertices[n][bad_vertex] = len(vertices) - 1 found = True - if found == False: - print("ERROR No suitable boundary point found!") - sys.exit() + assert (found is not False), 'No suitable boundary point found!' continue # -- add boundary edges, ignoring corner vertices for now @@ -480,48 +490,46 @@ def norm(v): v1 = ridge[0] v2 = ridge[1] if v1 >= new_vertices_start_index or v2 >= new_vertices_start_index: - if v1 >= new_vertices_start_index: - orig = v2 - new = v1 - else: - orig = v1 - 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 + 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]) + 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]) + 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]) + 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]) + 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) @@ -538,11 +546,12 @@ def norm(v): 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] + 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) + 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]): @@ -565,15 +574,14 @@ def norm(v): soft_equiv(vertices[v_indices[1]][1], ymax)): boundary_edges['yr'].append(ridge_idx) else: - print("ERROR Boundary edge not identified!") - sys.exit() + 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.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): From 4bc562215804030d51c506a0c813855677ecbed4 Mon Sep 17 00:00:00 2001 From: "Ben R. Ryan" Date: Thu, 5 Aug 2021 16:13:19 -0600 Subject: [PATCH 8/8] Modify base array types --- src/mesh/python/mesh_types.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/mesh/python/mesh_types.py b/src/mesh/python/mesh_types.py index 1cef6626b..bc7f7cdde 100644 --- a/src/mesh/python/mesh_types.py +++ b/src/mesh/python/mesh_types.py @@ -25,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 # ------------------------------------------------------------------------------------------------ #