tetgentofenics.py

#!/usr/bin/env python3
"""

Usage: ./tetgentofenics.py filename

    Note: This cript will look for the TetGen files filename.node and 
    filename.ele. If not found, an error messgae will be displayed.
    The Tetegn face file is not necessary.
    A new file called 'filename_from_tetgen.xml' will be created if all 
    goes well. Make sure you don't have a file with the same name you want to 
    keep before running this. 
 
    READING IN FENICS AND VISUALIZE
    ------------------------------------
    This script is particularly useful for multi-domains (i.e., multi-parts)
    geometries. To import the generated mesh in fenics and visualize the 
    geometry with the various domains do something like (assuming 
    output_mesh.xml was generated by this script)
    
    from fenics import *
    from dolfin import *
    
    mesh = Mesh("filename_from_tetgen.xml")
    mf = MeshFunction("size_t", mesh, 3, mesh.domains())
    domains_file =  File('output_mesh.pvd') 
    domains_file << mf

    Then, open output_mesh.pvd with Paraview and see the geometry and the 
    various domains.
    
    This script was tested with
    - Fenics/Dolfin 2019.1.0
    - Solidworks 2019 Education Edition
    - Python 3.6.10 on Ubuntu
    
@author: Alberto Corrias (Github: albertocorrias)
"""
import sys

if (len(sys.argv) != 2):
    raise Exception('Only filename expected. Usage: ./tetgentofenics.py filename')

tetgen_root_filename = sys.argv[1]

input_node_filename = tetgen_root_filename + '.node';
input_elem_filename = tetgen_root_filename + '.ele';
output_xml_filename = tetgen_root_filename + '_from_tetgen.xml';

try:
    input_node_file = open(input_node_filename, "r")
except:
    print('ERROR: Unable to open input file ' + input_node_filename)
    exit();

try:    
    input_ele_file = open(input_elem_filename, "r")
except:
    print('ERROR: Unable to open input file ' + input_elem_filename)
    exit();

try:    
    output_file = open(output_xml_filename, "w")
except:
    print('ERROR: Unable to open  output file')
    exit();    
    
#Headers and footers
header = '<?xml version="1.0"?> \n \
<dolfin xmlns:dolfin="http://fenicsproject.org"> \n \
<mesh celltype="tetrahedron" dim="3"> \n';

footer = '<data /> </mesh> \n \
</dolfin>\n'
output_file.write(header);
#########################

nodes=0;
elems = 0;
domain_tags = []
counter = 0
for line in input_node_file:
    split_line = line.split();
    if (split_line[0] != '#'):
        if counter==0:
            num_nodes = str(split_line[0])
            node_string_header = '<vertices size="' +num_nodes + '">\n';
            output_file.write(node_string_header)
        else:
            node_number = str(split_line[0])
            x = str(split_line[1])
            y = str(split_line[2])
            z = str(split_line[3])
            xml_node_line = '<vertex index="' + node_number+ \
            '" x="' + x +\
            '" y="' + y +\
            '" z="' + z +'" />\n'
            output_file.write(xml_node_line)
        counter = counter +1
    
node_string_footer = '</vertices> \n'
output_file.write(node_string_footer)
input_node_file.close()

###############Elements part
counter = 0;
num_elems = 0;
attributes = [];
for line in input_ele_file:
    split_line = line.split();
    if (split_line[0] != '#'):
        if counter==0:
            num_elems = split_line[0]
            num_elem_str = str(num_elems)
            elem_string_header = '<cells size="' + num_elem_str + '">\n';
            output_file.write(elem_string_header)
        else:
            ele_number = str(split_line[0])
            node_1 = str(split_line[1])
            node_2 = str(split_line[2])
            node_3 = str(split_line[3]) 
            node_4 = str(split_line[4]) 
            attributes.append(split_line[5])
            xml_elem_line = '<tetrahedron index="'+ele_number +\
            '" v0="' + node_1 +\
            '" v1="' + node_2 +\
            '" v2="' + node_3 +\
            '" v3="' + node_4 + '" />\n'
            output_file.write(xml_elem_line)
        counter = counter + 1
    
elem_string_footer = '</cells> \n'
output_file.write(elem_string_footer)
input_ele_file.close();

domain_string_header = '  <domains> \n \
<mesh_value_collection name="f" type="uint" dim="3" size="' + str(num_elems) + '">\n'
output_file.write(domain_string_header)

for el_idx in range(0,int(num_elems)):
    domain_element_string = '<value cell_index="' + str(el_idx) +\
    '" local_entity="0" value="'+str(attributes[el_idx])+'" />\n'
    
    output_file.write(domain_element_string)
    
    
domain_string_footer = '</mesh_value_collection> \n </domains>';
output_file.write(domain_string_footer)
#Footers and close
output_file.write(footer)
output_file.close();