convert_nastran_to_gmsh.py

import numpy as np

def read_nastran_mesh(file_path):
    nodes = {}
    elements = {}
    element_tags = {}
    physical_names = {}
    physical_names_dim = {}

    with open(file_path, "r") as file:
        lines = file.readlines()

        # Flag variables to determine whether to read nodes or elements
        reading_nodes = False
        reading_elements = False
        reading_physical_names = False

        num_nodes = 0
        num_elements = 0
        num_physical_names = 0

        for line in lines:
            if line.startswith("$ Node cards"):
                reading_nodes = True
                reading_elements = False
                reading_physical_names = False
                continue
            elif line.startswith("$ Element cards"):
                reading_elements = True
                reading_nodes = False
                reading_physical_names = False
                continue
            elif line.startswith("$ Property cards"):
                reading_elements = False
                reading_nodes = False
                reading_physical_names = True
                continue
            elif line.startswith("$ Material cards"):
                reading_elements = False
                reading_nodes = False
                reading_physical_names = False
                continue

            if reading_nodes:
                if line.startswith("GRID"):
                    parts = line.split()
                    node_id = int(parts[1])
                    x = float(parts[3])
                    y = float(parts[4])
                    z = float(parts[5]) if len(parts) > 3 else 0.0
                    nodes[node_id] = (x, y, z)
            elif reading_elements:
                if line.startswith("CQUAD") or\
                   line.startswith("CROD"):
                    parts = line.split()
                    element_id = int(parts[1])
                    element_nodes = [int(parts[i]) for i in range(3, len(parts))]
                    elements[element_id] = element_nodes
                    if element_id in element_tags:
                        element_tags[element_id].append(int(parts[2]))
                    else:
                        element_tags[element_id] = [int(parts[2])]
            elif reading_physical_names:
                if line.startswith("$ Name:"):
                    num_physical_names += 1
                    parts = line.split()
                    physical_name = parts[2]
                    physical_names[num_physical_names] = physical_name

                if line.startswith("PROD"):
                    pdim = 1
                    parts = line.split()
                    physical_id = int(parts[1])
                    physical_names_dim[physical_id] = pdim
                elif line.startswith("PSHELL"):
                    pdim = 2
                    parts = line.split()
                    physical_id = int(parts[1])
                    physical_names_dim[physical_id] = pdim

    # sort the nodes 
    sorted_nodes = dict(sorted(nodes.items()))


    return sorted_nodes, elements, element_tags, physical_names, physical_names_dim


def write_gmsh_mesh(file_path, nodes, elements, mesh_format,
                    physical_names, element_tags, physical_dim):
    with open(file_path, 'w') as file:
        # Write mesh format
        file.write("$MeshFormat\n")
        file.write(f"{mesh_format}\n")
        file.write("$EndMeshFormat\n")

        # Write physical names
        file.write("$PhysicalNames\n")
        file.write(f"{len(physical_names)}\n")
        for phys_id, phys_name in physical_names.items():
        #for i, (type_id, number, name) in enumerate(physical_names):
            #file.write(f"{type_id} {number} {name}\n")
            phys_dim = physical_dim[phys_id]
            file.write(f"{phys_dim} {phys_id} {phys_name}\n")
        file.write("$EndPhysicalNames\n")

        # Write nodes
        file.write("$Nodes\n")
        file.write(f"{len(nodes.items())}\n")
        for node_id, coordinates in nodes.items():
            file.write(f"{node_id} {coordinates[0]} {coordinates[1]} {coordinates[2]}\n")
        file.write("$EndNodes\n")

        # Write elements
        file.write("$Elements\n")
        file.write(f"{len(elements)}\n")
        for elem_id, element_nodes in elements.items():
        #for i, (element_nodes, element_tag) in enumerate(zip(elements, element_tags)):
            if len(element_nodes) == 2:
                   element_type = 1
            elif len(element_nodes) == 4:
                   element_type = 3
            elif len(element_nodes) == 8:
                   element_type = 5
            element_type = 3 if len(element_nodes) == 4 else 5  # Determine element type based on number of nodes
            tags = element_tags[elem_id]
            num_tags = len(tags)
            file.write(f"{elem_id} {element_type} {num_tags} {' '.join(map(str, tags))} {' '.join(map(str, element_nodes))}\n")
        file.write("$EndElements\n")

def sort_lexicographic(element_nodes, nodes):
    # Sort hexahedral element nodes based on their coordinates
    # This maintains the original gmsh order, but puts the first node
    # at the minimum x, y, z coordinate
    sorted_nodes = sorted(element_nodes,
                          key=lambda node_id: (nodes[node_id][2]))

    # sort the min and max z seperately for y
    sorted_nodes[:4] = sorted(sorted_nodes[:4],
                          key=lambda node_id: float(nodes[node_id][1]))
    sorted_nodes[4:] = sorted(sorted_nodes[4:],
                          key=lambda node_id: float(nodes[node_id][1]))

    # sort by x depending on which y extrema the nodes live on
    sorted_nodes[:2] = sorted(sorted_nodes[:2],
                          key=lambda node_id: float(nodes[node_id][0]))
    sorted_nodes[2:4] = sorted(sorted_nodes[2:4],
                          key=lambda node_id: float(nodes[node_id][0]))
    sorted_nodes[4:6] = sorted(sorted_nodes[4:6],
                          key=lambda node_id: float(nodes[node_id][0]))
    sorted_nodes[6:] = sorted(sorted_nodes[6:],
                          key=lambda node_id: float(nodes[node_id][0]))

    return sorted_nodes

def sort_gmsh(element_nodes, nodes):
    # Sort hexahedral element nodes based on their coordinates
    # This maintains the original gmsh order, but puts the first node
    # at the minimum x, y, z coordinate
    sorted_nodes = sorted(element_nodes,
                          key=lambda node_id: (nodes[node_id][2]))

    # sort the min and max z seperately for y
    sorted_nodes[:4] = sorted(sorted_nodes[:4],
                          key=lambda node_id: float(nodes[node_id][1]))
    sorted_nodes[4:] = sorted(sorted_nodes[4:],
                          key=lambda node_id: float(nodes[node_id][1]))

    # sort by x depending on which y extrema the nodes live on
    sorted_nodes[:2] = sorted(sorted_nodes[:2],
                          key=lambda node_id: float(nodes[node_id][0]))
    sorted_nodes[2:4] = sorted(sorted_nodes[2:4],
                          key=lambda node_id: float(nodes[node_id][0]), reverse=True)
    sorted_nodes[4:6] = sorted(sorted_nodes[4:6],
                          key=lambda node_id: float(nodes[node_id][0]))
    sorted_nodes[6:] = sorted(sorted_nodes[6:],
                          key=lambda node_id: float(nodes[node_id][0]), reverse=True)

    return sorted_nodes

def sort_gmsh2(element_nodes, nodes):
    # Sort hexahedral element nodes based on their coordinates
    # This maintains the original gmsh order, but puts the first node
    # at the minimum x, y, z coordinate
    # first sort based on x
    sorted_nodes = sorted(element_nodes,
                          key=lambda node_id: (nodes[node_id][0]))

    # sort the min and max x nodes seperately for y
    sorted_nodes[:4] = sorted(sorted_nodes[:4],
                          key=lambda node_id: float(nodes[node_id][1]))
    sorted_nodes[4:] = sorted(sorted_nodes[4:],
                          key=lambda node_id: float(nodes[node_id][1]))

    # sort by z depending on which (x,y) extrema the nodes live on
    sorted_nodes[:2] = sorted(sorted_nodes[:2],
                          key=lambda node_id: float(nodes[node_id][2]))
    sorted_nodes[2:4] = sorted(sorted_nodes[2:4],
                          key=lambda node_id: float(nodes[node_id][2]), reverse=True)
    sorted_nodes[4:6] = sorted(sorted_nodes[4:6],
                          key=lambda node_id: float(nodes[node_id][2]))
    sorted_nodes[6:] = sorted(sorted_nodes[6:],
                          key=lambda node_id: float(nodes[node_id][2]), reverse=True)

    return sorted_nodes

def sort_gmsh3(element_nodes, nodes):
    # Sort hexahedral element nodes based on their coordinates
    # This maintains the original gmsh order, but puts the first node
    # at the minimum x, y, z coordinate
    # first sort based on x
    #print(f"{element_nodes=}")
    sorted_nodes = element_nodes.copy()
    sorted_nodes[1] = element_nodes[4]
    sorted_nodes[4] = element_nodes[1]
    sorted_nodes[7] = element_nodes[2]
    sorted_nodes[2] = element_nodes[7]
    #print(f"{element_nodes=}")
    #print(f"{sorted_nodes=}")

    return sorted_nodes

# Example usage:
file_path = "coarse.msh"  # Change this to your desired output file path
original_file_path = "coarse.bdf"  # Path to the original Gmsh mesh file
nodes, elements, element_tags, physical_names, physical_dim = read_nastran_mesh(original_file_path)

"""
# detect any backwards elements
for element in elements:
    if len(element) == 8:  # Hex elements only
        r = np.asarray(nodes[element[1]]) - np.asarray(nodes[element[0]])
        r_hat = r/np.linalg.norm(r)
        s = np.asarray(nodes[element[2]]) - np.asarray(nodes[element[1]])
        s_hat = s/np.linalg.norm(s)
        t = np.asarray(nodes[element[4]]) - np.asarray(nodes[element[0]])
        t_hat = t/np.linalg.norm(t)
        t_candidate = np.cross(r_hat, s_hat)
        t_candidate_hat = t_candidate/np.linalg.norm(t_candidate)
        #if (t_candidate_hat != t_hat).all:
        #if not np.allclose(t_candidate_hat, t_hat):
        if np.dot(t_hat, t_candidate_hat) < 0.:
            print(f"{s_hat=}, {r_hat=}, {t_hat=}, {t_candidate_hat=}")

# sort the element verticies for processing by mirgecom/meshmode
sorted_elements = []
for element in elements:
    if len(element) == 8:  # Hex elements only
        sorted_element = sort_gmsh3(element, nodes)
        sorted_element = sort_gmsh3(element, nodes)
    else:
        sorted_element = element
    sorted_elements.append(sorted_element)

elements = sorted_elements
"""

# This is the mesh format meshmode expects
mesh_format = "2.2 0 8"

print("Nodes:")
for node_id, coordinates in nodes.items():
    print(f"Node {node_id}: {coordinates}")

print("\nElements:")
for element_id, element_nodes in elements.items():
    print(f"Element Nodes for Element {element_id}: {element_nodes}")

print("\nElement Tags:")
for element_id, tags in element_tags.items():
    print(f"Element Tags for Element {element_id}: {tags}")

print("\nPhysical Names:")
for name_id, name in physical_names.items():
    print(f"Physical name {name_id}: {name}")

write_gmsh_mesh(file_path, nodes, elements, mesh_format, physical_names, element_tags, physical_dim)