I am going from a triangular mesh coarsened to a polygonal mesh by taking shared edges between two triangles and connecting them. I am finding an error with the order of connection. I don’t know if this is a problem with my generation or my coarsening algorithm. Here is both pieces of code:
Coarsening Algorithm:
from collections import defaultdict
import numpy as np
def find_shared_edges(cell_to_node):
"""Find shared edges between triangles and store triangle indices that share the edges."""
edge_to_triangles = defaultdict(list)
shared_edges = []
# Loop over each triangle in the mesh
for i, triangle in enumerate(cell_to_node):
edges = get_edges(triangle)
for edge in edges:
edge_key = tuple(sorted(edge)) # Sorting ensures consistency when comparing edges
edge_to_triangles[edge_key].append(i)
# Identify shared edges between two triangles
for edge, triangles in edge_to_triangles.items():
if len(triangles) == 2: # Shared edge between exactly two triangles
shared_edges.append((triangles[0], triangles[1], edge))
return shared_edges
def get_edges(triangle):
"""Return the edges of a triangle, represented as a list of tuples of vertex coordinates."""
return [(triangle[i], triangle[(i + 1) % 3]) for i in range(3)]
def merge_triangles(tri1, tri2, shared_edge):
"""Merge two triangles into one by combining their shared edge and the remaining vertices,
following the specified order: shared edge, non-shared vertex, shared edge, remaining vertex."""
# Get the vertices that are not part of the shared edge
tri1_remaining = [v for v in tri1 if v not in shared_edge]
tri2_remaining = [v for v in tri2 if v not in shared_edge]
# Ensure the order of the new triangle: shared edge, non-shared vertex, shared edge, remaining vertex
# shared_edge consists of two vertices, tri1_remaining and tri2_remaining each contain one vertex
new_triangle = [shared_edge[0], tri1_remaining[0], shared_edge[1], tri2_remaining[0]]
return new_triangle
def coarsen_mesh(cell_to_node):
"""Coarsen the mesh by merging triangles with shared edges."""
used_cells = set() # Track merged triangles
shared_edges = find_shared_edges(cell_to_node) # Find shared edges between triangles
coarsened_mesh = [] # Store the new mesh of merged triangles
for i, j, shared_edge in shared_edges:
if i not in used_cells and j not in used_cells:
# Merge the triangles into one new triangle
merged_triangle = merge_triangles(cell_to_node[i], cell_to_node[j], shared_edge)
# Add the merged triangle to the coarsened mesh
coarsened_mesh.append(merged_triangle)
# Mark the original triangles as merged
used_cells.add(i)
used_cells.add(j)
# Add remaining triangles that were not merged
for i, triangle in enumerate(cell_to_node):
if i not in used_cells:
coarsened_mesh.append(triangle)
return coarsened_mesh
# Example input: A list of triangles with coordinates (example data)
# Coarsen the mesh
coarsened_cell_to_node = coarsen_mesh(cell_to_node)
# Convert coarsened mesh to a NumPy array and print its shape
coarsened_cell_to_node = np.array(coarsened_cell_to_node, dtype=object) # Use dtype=object for non-homogeneous arrays
print(coarsened_cell_to_node)
# Output the coarsened mesh
#print("Coarsened Mesh:")
#for triangle in coarsened_mesh:
# print(triangle)
Mesh Generation:
# Define geometric and topological dimensions
gdim = 2
shape = "quadrilateral"
degree = 1
# Define the cell and domain
cell = ufl.Cell(shape, geometric_dimension=gdim)
domain = ufl.Mesh(ufl.VectorElement("Lagrange", cell, degree))
# Create the mesh
mesh = dolfinx.mesh.create_mesh(MPI.COMM_WORLD, coarsed_test, new_coords, domain)
tdim = mesh.topology.dim
fdim = tdim - 1
mesh.topology.create_connectivity(fdim, tdim)
# Check if PyVista backend is available