8.4  Parsing a mesh

import numpy as np
import mechanicskit as mk

# Parse nodes and elements from file
node_coords = []
connectivity = []

with open('../assets/100x100mm_8_linear_trias.pc', 'r') as file:
    for line in file:
        if line.startswith('NODE'): # Collect node coordinates
            parts = line.split()
            NID = int(parts[1])
            NX = float(parts[2])
            NY = float(parts[3])
            NZ = float(parts[4])
            node_coords.append([NX, NY, NZ])
        elif line.startswith('SHELL'): # Collect element connectivity (1-based)
            parts = line.split()
            EID = int(parts[1])
            N1 = int(parts[3])
            N2 = int(parts[4])
            N3 = int(parts[5])
            connectivity.append([N1, N2, N3])

# Create Mesh object with 1-based indexing
mesh = mk.Mesh(node_coords, connectivity)

print(f"Mesh info: {mesh.n_nodes} nodes, {mesh.n_elements} elements")
print(f"Element type: {mesh.element_type}")

# Compute total area
area = 0
for iel, nodes, coords in mesh.iter_elements():
    # iter_elements() yields: element_number, node_IDs, node_coordinates
    # Extract the three node coordinates for this triangular element
    n1 = coords[0]  # first node coordinates
    n2 = coords[1]  # second node coordinates
    n3 = coords[2]  # third node coordinates

    # Compute edge vectors and area
    u = n2 - n1  # triangle edge vector 1
    v = n3 - n1  # triangle edge vector 2
    w = np.cross(u, v)  # perpendicular to triangle
    area_e = 0.5 * np.linalg.norm(w)  # area of current element
    area += area_e  # add to total area

print(f"Total computed area: {area:.2f} mm²")
print(f"Expected area: {100*100:.2f} mm² (100mm × 100mm square)")
print(f"Verification: Error = {abs(area - 10000):.6f} mm²")

# Parse circle mesh from .k file
node_dict = {}  # Map node ID to coordinates
connectivity_circle = []

# Track which section we're in
in_node_section = False
in_element_section = False

with open('../assets/unit_circle_cirk_pattern.k', 'r') as file:
    for line in file:
        # Check for section headers
        if line.startswith('*NODE'):
            in_node_section = True
            in_element_section = False
            continue
        elif line.startswith('*ELEMENT_SHELL'):
            in_node_section = False
            in_element_section = True
            continue
        elif line.startswith('*'):
            # New section started, stop parsing
            in_node_section = False
            in_element_section = False
            continue

        # Parse node data (fixed-width format)
        if in_node_section and line.strip():
            # Fixed-width columns: 8 chars for ID, 16 chars each for X, Y, Z
            NID = int(line[0:8])
            NX = float(line[8:24])
            NY = float(line[24:40])
            NZ = float(line[40:56])
            node_dict[NID] = [NX, NY, NZ]

        # Parse element data (space-separated)
        elif in_element_section and line.strip():
            parts = line.split()
            EID = int(parts[0])
            # Column 1 is material/property ID, columns 2-4 are the three unique node IDs
            # (Column 5 repeats the last node for triangular elements)
            N1 = int(parts[2])
            N2 = int(parts[3])
            N3 = int(parts[4])
            connectivity_circle.append([N1, N2, N3])

# Build sequential node coordinate list and update connectivity to match
# Map original node IDs to sequential indices (1-based for Mesh class)
node_id_to_index = {nid: idx+1 for idx, nid in enumerate(sorted(node_dict.keys()))}
node_coords_circle = [node_dict[nid] for nid in sorted(node_dict.keys())]

# Remap connectivity from original node IDs to sequential indices
connectivity_remapped = [[node_id_to_index[n1], node_id_to_index[n2], node_id_to_index[n3]]
                         for n1, n2, n3 in connectivity_circle]

# Create Mesh object with 1-based indexing
mesh_circle = mk.Mesh(node_coords_circle, connectivity_remapped)

print(f"Circle mesh: {mesh_circle.n_nodes} nodes, {mesh_circle.n_elements} elements")
print(f"Node ID range: {min(node_dict.keys())} to {max(node_dict.keys())}")

# Prepare data for visualization using Mesh
P_circle = mesh_circle.nodes[:, :2]  # Extract x, y coordinates from Mesh

# Visualize
fig, ax = plt.subplots(figsize=(7, 7))

# Plot mesh using mesh.elements (1-based connectivity) directly
mk.patch('Faces', mesh_circle.elements, 'Vertices', P_circle,
         'FaceColor', 'white',
         'EdgeColor', 'black',
         'EdgeAlpha', 0.3,
         'LineWidth', 1.0,
         ax=ax)

# Plot nodes
ax.plot(P_circle[:, 0], P_circle[:, 1], 'o', color='blue', markersize=4)

# Add grid for verification
ax.grid(True, alpha=0.5)
ax.set_aspect('equal')
ax.set_xlim(-1.2, 1.2)
ax.set_ylim(-1.2, 1.2)

# Draw unit circle outline for reference
theta = np.linspace(0, 2*np.pi, 100)
ax.plot(np.cos(theta), np.sin(theta), 'r--', linewidth=2, alpha=0.5, label='Unit circle')

ax.set_xlabel('x', fontsize=12)
ax.set_ylabel('y', fontsize=12)
ax.set_title(f'Unit Circle Mesh ({mesh_circle.n_elements} elements)\n' +
             f'{mesh_circle.n_nodes} nodes',
             fontsize=14)
ax.legend()
plt.show()