Conversion of ABAQUS .inp mesh into xdmf with surface data

[JohannesMenning]

Hello Everyone,

I want to convert a 3D ABAQUS Mesh with definition of element and node sets into the xdmf format.

It works perfectly, using e.g. the following Code, when I just want to convert the mesh and preserve the information of the subdomains:

    import numpy as np
    import meshio
    filename='MME-Cube'
    inp_file_path = filename+".inp"
    
    counter = 1
    mesh = meshio.read(inp_file_path)
    cells = mesh.get_cells_type("hexahedron")
    cell_sets = mesh.cell_sets_dict
    cell_data = np.zeros(len(cells))
    for marker, set in cell_sets.items():
        for type, entities in set.items():
            if marker == 'Set-domain':
                    cell_data[entities] = int(counter)
                    counter += 1
    
    out_mesh = meshio.Mesh(points=mesh.points, cells={
                           "hexahedron": cells}, cell_data={"cell_data": [cell_data]})
    meshio.write(filename+".xdmf", out_mesh)

However, when I try to add the information about the surfaces, I get stuck.

I have read several posts, e.g. this post or this or this.

I guess, the problem is, that usually, a 2D surface mesh is required, which is not provided by ABAQUS.

This is the Code which I tried so far, however when loading the xdmf file in dolfinx I can not access "cell_tags" or "facet_tags"

Here is my 2nd Code:

import meshio
import numpy as np

def read_inp(file_path):
    nodes = []
    elements = {}
    node_sets = {}
    element_sets = {}
    
    with open(file_path, 'r') as file:
        lines = file.readlines()
    
    for i, line in enumerate(lines):
        if line.startswith('*Node'):
            j = i + 1
            while j < len(lines) and not lines[j].startswith('*'):
                parts = lines[j].strip().split(',')
                nodes.append([float(p) for p in parts[1:]])
                j += 1
        elif line.startswith('*Element'):
            element_type = line.strip().split('=')[1]
            elements[element_type] = []
            j = i + 1
            while j < len(lines) and not lines[j].startswith('*'):
                parts = lines[j].strip().split(',')
                elements[element_type].append([int(p) for p in parts[1:]])
                j += 1
        elif line.startswith('*Nset'):
            set_name = line.strip().split('=')[1]
            node_sets[set_name] = []
            j = i + 1
            while j < len(lines) and not lines[j].startswith('*'):
                node_sets[set_name].extend([int(p) for p in lines[j].strip().split(',')])
                j += 1
        elif line.startswith('*Elset'):
            set_name = line.strip().split('=')[1]
            element_sets[set_name] = []
            j = i + 1
            while j < len(lines) and not lines[j].startswith('*'):
                element_sets[set_name].extend([int(p) for p in lines[j].strip().split(',')])
                j += 1
    
    return nodes, elements, node_sets, element_sets

def convert_to_meshio(nodes, elements, node_sets, element_sets, filename):
    points = np.array(nodes)
    cells = []
    cell_data = {"cell_tags": []}
    facet_tags = np.zeros(len(points), dtype=int)

    element_type_mapping = {
        "C3D4": "tetra",
        "C3D8R": "hexahedron",
        "C3D8": "hexahedron",
        "C3D10": "tetra10",
        "C3D15": "tetra15",
        # Add more element type conversions as needed
    }

    for element_type, element_list in elements.items():
        cell_type = element_type_mapping.get(element_type)
        if cell_type:
            cells.append((cell_type, np.array(element_list) - 1))  # Adjust for zero-based indexing

    # Initialize tags for each cell type
    for cell_type, cell_array in cells:
        cell_tags = np.zeros(len(cell_array), dtype=int)

    subdomain_marker_dict = {}
    count_subdomains = 1


    ######################################################################################
    mesh_abaqus = meshio.read(filename+'.inp') 
    # Populate the marker dictionaries for cell_data
    for key, value in mesh_abaqus.cell_sets_dict.items():
        if key.startswith('Set-domain'):
            subdomain_marker_dict[key] = count_subdomains
            count_subdomains += 1

    cell_data_subdomains = np.zeros(len(cell_array), dtype=int)


    # Assign marker data
    for marker, abq_set in mesh_abaqus.cell_sets_dict.items():
        for abq_celltype, entities in abq_set.items():
            if abq_celltype == 'hexahedron':
                if marker in subdomain_marker_dict:
                    cell_data_subdomains[entities] = int(subdomain_marker_dict[marker])
    cell_data = {"cell_tags": [cell_data_subdomains]}
    
    ####
    # Creating facet tags from node sets
    lines_surfaces = []
    for set_name, set_nodes in node_sets.items():
        if not set_name.startswith('Set-domain'):
            tag = list(node_sets.keys()).index(set_name) + 100 # to distinguish surface from cell markers
            lines_surfaces.append(f"{set_name} = {tag}")
            for node in set_nodes:
                facet_tags[node - 1] = tag  # Adjust for zero-based indexing

    with open(filename + "_facets.txt", 'w') as file:
        for line in lines_surfaces:
            file.write(line+'\n')

    facet_data = {"facet_tags": facet_tags}

    mesh = meshio.Mesh(points, cells, cell_data=cell_data, point_data=facet_data)
    meshio.write(filename+".xdmf", mesh, file_format="xdmf", data_format="HDF")



if __name__ == "__main__":

    
    filename='MME-Cube'
    inp_file_path = filename+".inp"

    nodes, elements, node_sets, element_sets = read_inp(inp_file_path)
    mesh = convert_to_meshio(nodes, elements, node_sets, element_sets, filename)

And Here is a minimal input file:

*Heading
** Job name: MME-Cube Model name: Model-1
** Generated by: Abaqus/CAE Unofficial Packaging Version
*Preprint, echo=NO, model=NO, history=NO, contact=NO
**
** PARTS
**
*Part, name=MME-Cube
*Node
      1,           0.,          15.,          20.
      2,           0.,          7.5,          20.
      3,           0.,           0.,          20.
      4,           0.,          15.,    13.333333
      5,           0.,          7.5,    13.333333
      6,           0.,           0.,    13.333333
      7,           0.,          15.,   6.66666651
      8,           0.,          7.5,   6.66666651
      9,           0.,           0.,   6.66666651
     10,           0.,          15.,           0.
     11,           0.,          7.5,           0.
     12,           0.,           0.,           0.
     13,         -7.5,          15.,          20.
     14,         -7.5,          7.5,          20.
     15,         -7.5,           0.,          20.
     16,         -7.5,          15.,    13.333333
     17,         -7.5,          7.5,    13.333333
     18,         -7.5,           0.,    13.333333
     19,         -7.5,          15.,   6.66666651
     20,         -7.5,          7.5,   6.66666651
     21,         -7.5,           0.,   6.66666651
     22,         -7.5,          15.,           0.
     23,         -7.5,          7.5,           0.
     24,         -7.5,           0.,           0.
     25,         -15.,          15.,          20.
     26,         -15.,          7.5,          20.
     27,         -15.,           0.,          20.
     28,         -15.,          15.,    13.333333
     29,         -15.,          7.5,    13.333333
     30,         -15.,           0.,    13.333333
     31,         -15.,          15.,   6.66666651
     32,         -15.,          7.5,   6.66666651
     33,         -15.,           0.,   6.66666651
     34,         -15.,          15.,           0.
     35,         -15.,          7.5,           0.
     36,         -15.,           0.,           0.
*Element, type=C3D8R
 1, 13, 14, 17, 16,  1,  2,  5,  4
 2, 14, 15, 18, 17,  2,  3,  6,  5
 3, 16, 17, 20, 19,  4,  5,  8,  7
 4, 17, 18, 21, 20,  5,  6,  9,  8
 5, 19, 20, 23, 22,  7,  8, 11, 10
 6, 20, 21, 24, 23,  8,  9, 12, 11
 7, 25, 26, 29, 28, 13, 14, 17, 16
 8, 26, 27, 30, 29, 14, 15, 18, 17
 9, 28, 29, 32, 31, 16, 17, 20, 19
10, 29, 30, 33, 32, 17, 18, 21, 20
11, 31, 32, 35, 34, 19, 20, 23, 22
12, 32, 33, 36, 35, 20, 21, 24, 23
*Nset, nset=Set-top
  1,  2,  3, 13, 14, 15, 25, 26, 27
*Elset, elset=Set-top
 1, 2, 7, 8
*Nset, nset=Set-bottom
 10, 11, 12, 22, 23, 24, 34, 35, 36
*Elset, elset=Set-bottom
  5,  6, 11, 12
*Nset, nset=Set-domain, generate
  1,  36,   1
*Elset, elset=Set-domain, generate
  1,  12,   1
*End Part
**  
**
** ASSEMBLY
**
*Assembly, name=Assembly
**  
*Instance, name=MME-Cube-1, part=MME-Cube
*End Instance
**  
*End Assembly

I appreciate any tips how I could solve this problem.

Best regards
Martin