Ex2Ab_modules.py

import sys
import re
import os
import pathlib
import math
import numpy as np
import time
import tkinter 
from tkinter import filedialog
from netCDF4 import Dataset

# Creates GUI for selecting the mesh file.
def get_mesh(filename):

    if filename is None:
        # Function for file selection.
        def chooseFile():
            main_win.mesh_file = filedialog.askopenfilename(
                parent = main_win, initialdir = '\\', 
                title = 'Please select a directory', 
                filetypes = (
                ('Exodus II', '.exo'),
                ('Exodus II', '.e'),
                ('Exodus II', '.g'),
                ('All files', '*.*')))
        
        # Function for closing selection window.
        def end():
            main_win.destroy()

        main_win = tkinter.Tk()
        main_win.geometry('250x200')
        main_win.mesh_file = ''
        b_choose_file = tkinter.Button(
            main_win, text = 'Choose Mesh File', width = 20, height = 3, 
            command = chooseFile)
        b_choose_file.place(x = 50, y = 30)
        b_choose_file.width = 100
        b_done = tkinter.Button(
            main_win, text = 'Done', width = 20, height = 3, 
            command = end)
        b_done.place(x = 50, y = 110)
        b_done.width = 100
        main_win.mainloop()

        valid_mesh = True
        if (not main_win.mesh_file):
            print('Ex2Ab has been terminated by user')
            sys.exit()
        
        try:
            rootgrp = Dataset(main_win.mesh_file, 'r', 
                            format = 'NETCDF3_64BIT_OFFSET')
        except:
            valid_mesh = False
        mesh = [main_win.mesh_file, valid_mesh]
    else:
        valid_mesh = True
        try:
            rootgrp = Dataset(filename, 'r', 
                            format = 'NETCDF3_64BIT_OFFSET')
        except:
            valid_mesh = False
        mesh = [filename, valid_mesh]
    return(mesh)

# Function for the runtime in seconds.
def timing(start_time):
    runtime = (time.time() - start_time)
    rounded_runtime = round(runtime, 2)
    return rounded_runtime


# This is to check that all blocks have compatible elements. Cubit doesn't seem
# to have pentahedral elements like Abaqus. If other codes writing out to 
# Exodus II tag the elements differently, these labels need to be updated.
def check_elements(blocks, rootgrp, el_types):
    abort = 0
    for i in range(blocks):
        a = str(rootgrp.get_variables_by_attributes(
                name = 'connect' + str(i + 1)))
        el_check = re.search('elem_type: (.+)\\n', a)
        el_type = el_check.group(1)
        el_types.append(el_type)
        if ((el_type != 'HEX8') and (el_type != 'TETRA4') 
            and (el_type != 'TETRA10') and (el_type != 'HEX20') 
            and (el_type != 'HEX') and (el_type != 'TETRA')):
            print('**ERROR** Block ' + str(i + 1) 
                  + ' has incompatible elements! Type: ' + str(el_type))
            abort = 1
    if (abort == 1):
        print('You should rethink your meshing choices...' 
              + 'and perhaps your life choices\n')
        sys.exit()
        
# Name files based on source mesh.
def name_files(file):
    ext = pathlib.Path(file).suffix
    name = file[:-len(ext)]
    # Check if named file already exists. Name will be appended with "_#" to 
    # avoid overwriting.
    exist = os.path.isfile(name + '.inp')
    h = 0
    while exist is True:
        h = h + 1
        job = name + '_' + str(h)
        exist = os.path.isfile(job + '.inp')
    return job

# Get attributes and material definitions.
def get_materials(blocks, rootgrp, mats):
    # Check if attributes were defined for blocks. Attributes are used to 
    # specify density.
    att = True
    try: ATT = rootgrp.variables['attrib1']
    except: att = False
    if (att == True):
        print('Density attribute was defined for blocks')
    if (att == False):
        print('Default densities are being used')
    
    for i in range(blocks):
        if (att == True):
            ATT = rootgrp.variables['attrib' + str(i + 1)]
            mats[i] = str(ATT[0]).strip('[]')
        if (att == False):
            # Default densities are arbitrarily set to 1g/cc.
            mats[i] = -1.0 
        if (blocks <= 100):
            print('Cell ' + str(i + 1) + ' density: ' 
                   + format((mats[i]), '1.10E'))
    return(att)

# Get the node coordinates for the mesh.
def get_nodes(rootgrp):
    exodus = True
    try: xyz = np.array(rootgrp.variables['coord'])
    except: exodus = False
    # This can probably be combined with previous try statement.
    # Thought coord was [n, 3], but this ensures it is treated that
    # way even if it's [3, n].
    if (exodus == True and len(xyz[:, 0]) == 3):
        xyz = xyz.T
    if (exodus == False):
        len_node = len(rootgrp.variables['coordx'])
        xyz = np.zeros([len_node, 3])
        xyz[:, 0] = np.array(rootgrp.variables['coordx'])
        xyz[:, 1] = np.array(rootgrp.variables['coordy'])
        xyz[:, 2] = np.array(rootgrp.variables['coordz'])
    return(xyz)

# Set up and write the MCNP skeletal input file.
def mcnp(rootgrp, blocks, att, xyz, job, mats):
    # Check if node sets exist.
    reflect = True
    try: ns = rootgrp.variables['ns_prop1']
    except: reflect = False
    # Get number of node sets, if they exist.
    if (reflect == True):
        sets = len(ns)
    if (reflect == True):
        print(str(sets) 
              + ' reflecting boundaries will be defined from node sets')
    
    # Find max and min xyz positions to find model centroid.
    dx = np.max(xyz[:, 0]) - np.min(xyz[:, 0])
    dy = np.max(xyz[:, 1]) - np.min(xyz[:, 1])
    dz = np.max(xyz[:, 2]) - np.min(xyz[:, 2])
    x0 = np.min(xyz[:, 0]) + (0.5 * dx)
    y0 = np.min(xyz[:, 1]) + (0.5 * dy)
    z0 = np.min(xyz[:, 2]) + (0.5 * dz)
    d = np.sqrt((dx**2) + (dy**2) + (dz**2))
    # Radius is slightly larger than needed to contain model.
    rad = 1.01 * (d / 2)
    
    # This section defines surfaces to create reflecting boundaries. This is 
    # only done if the user specified node sets. The node set operations will 
    # only work properly for planar surfaces defined on the exterior of the 
    # model. Considering the plausible uses of reflecting boundaries, this 
    # should not be a problem.
    if (reflect == True):
        # Stores xyz vector for each node set.
        plane = np.zeros([sets, 3]) 
        # Stores distance for each vector.
        dist = np.zeros([sets]) 
        p_max = np.zeros([sets])
        p_min = np.zeros([sets])
        # Offset distance for each reflecting boundary.
        offset = np.zeros([sets]) 
        for i in range(sets):
            p_mult = np.zeros([3, 3])
            nodes = np.array(rootgrp.variables['node_ns' + str(i + 1)])
            nnodes = len(nodes)
            mid = math.floor(nnodes / 2.0)
            points = np.zeros([3])
            # Store the first, middle, and last node from set. These are the 3 
            # points used to define a plane.
            points[0] = nodes[0]
            points[1] = nodes[mid]
            points[2] = nodes[nnodes - 1]
            coord = np.zeros([3, 3])
            # Retrieve the coordinates of the 3 points.
            for p in range(3):
                point = int(points[p])
                coord[p, 0] = np.float(xyz[ (point - 1), 0])
                coord[p, 1] = np.float(xyz[ (point - 1), 1])
                coord[p, 2] = np.float(xyz[ (point - 1), 2])
            # Xyz vectors are found for points 1 to 2 and points 1 to 3.
            v1 = np.zeros([3])
            v2 = np.zeros([3])
            for v in range(3):
                v1[v] = coord[1, v] - coord[0, v]
                v2[v] = coord[2, v] - coord[0, v]
            # Cross product gives normal vector of the plane.
            nv = np.cross(v1, v2) 
            D = 0
            # Calculate the plane distance. Distance is D in Ax+By+Cz=D.
            for n in range(3):
                plane[i][n] = nv[n]
                D = D + (coord[0, n] * plane[i, n])
            p_max = np.amax(plane, axis = 1)
            p_min = np.amin(plane, axis = 1)
            # Normalize the plane equation coefficients.
            for n in range(3):
                if (abs(p_max[i]) > abs(p_min[i])):
                    plane[i, n] = plane[i, n] / p_max[i]
                    dist[i] = D / p_max[i]
                if (abs(p_max[i]) < abs(p_min[i])):
                    plane[i, n] = plane[i, n] / p_min[i]
                    dist[i] = D / p_min[i]
            # Fill array for testing plane position after applying offset. Want 
            # to find the plane that is beyond the model boundaries.
            for a in range(3):
                p_mult[0, a] = plane[i, a] * dist[i]
                p_mult[1, a] = plane[i, a] * (dist[i] + 1e-5)
                p_mult[2, a] = plane[i, a] * (dist[i] - 1e-5)
            # Find distances from model center to planes 
            dpc = np.zeros([3])
            for a in range(3):
                dpc[a] = np.sqrt(((x0 - p_mult[a, 0])**2) 
                                 + ((y0 - p_mult[a, 1])**2) 
                                 + ((z0 - p_mult[a, 2])**2))
            # Compare offset plane distance to initial plane distance. Choose 
            # offset that moves plane further from the model center.
            if (dpc[1] > dpc[0]):
                offset[i] = 1e-5
            if (dpc[2] > dpc[0]):
                offset[i] = -1e-5
    
    # Begin writing the MCNP skeletal input to file.
    mcnp = open(job + '.mcnpinp', 'w')
    mcnp.write('c These are the pseudo-cells for running MCNP with a UM\n' 
               + 'c Created from file: ' + str(job) + '.inp\nc \nc \n' 
               + 'c PSEUDO-CELLS\n' 
               + 'c Each of these cells is one block (element set) '
               + 'from your UM\n')
    if (att == False):
        mcnp.write('c Please add accurate material densities for each ' 
                   + 'pseudo-cell\n')
        
    # Pseudo-cells are created from each block.
    for i in range(blocks):
        mcnp.write(format((i + 1),'0>2d') + '        ' + str(i + 1)
                   + '     ' + format(mats[i], '1.10E') + '  0  u=1\n')
    
    # Background cell.
    mcnp.write('c This is the mesh universe background cell,' 
               + ' it is void by default\n' 
               + 'c It fills in any space between your UM geometry' 
               + ' and the universe boundary surface\n' 
               + format((blocks + 1), '0>2d') 
               + '        0                        0'
               + '  u=1\n')
    
    # Legacy cells.
    # Cell containing universe.
    mcnp.write('c \nc LEGACY CELLS\n' 
               + 'c This cell contains the mesh universe with all of your UM' 
               + ' geometry\n' 
               + format((blocks + 2), '0>2d') 
               + '        0    -99')
    
    # Write additional cell definitions if node sets were used.
    if (reflect == True): 
        for h in range(sets):
            if (offset[h] > 0):
                mcnp.write(' -' + str(98 - h))
            if (offset[h] < 0):
                mcnp.write(' ' + str(98 - h))
    
    mcnp.write('  fill=1\n')
    
    # Creates cell between UM geometry and reflecting boundaries.
    if (reflect == True):
        mcnp.write('c This is the space between your UM universe boundary' 
                   + ' and the reflecting boundaries\n' 
                   + format((blocks + 3), '0>2d') + '        0    -99')
        for h in range(sets):
            if (offset[h] > 0):
                mcnp.write(' -' + str(98 - sets - h))
            if (offset[h] < 0):
                mcnp.write(' ' + str(98 - sets - h))
        mcnp.write(' #' + format((blocks + 2), '0>2d') + '\n')
    
    # Termination region outside universe.
    mcnp.write('c This cell is the space outside the UM universe\n' 
               + format((blocks + 4), '0>2d') 
               + '        0     #' 
               + format((blocks + 2), '0>2d') 
               + ' #' + format((blocks + 3), '0>2d')
               + '\n\nc \n' 
               + 'c SURFACES\n')
    
    if (reflect == True):
        mcnp.write('c Additional surfaces for reflecting boundaries were' 
                   + ' defined from side sets\n' 
                   + 'c All reflecting surfaces are intentionally offset from' 
                   + ' the UM universe boundary\nc \n')
    # Create surfaces based on node sets.
        for k in range(sets):
            p1 = [ format(plane[k, 0], '1.10E'), 
                   format(plane[k, 1], '1.10E'),
                   format(plane[k, 2], '1.10E'),
                   format((dist[k] + offset[k]), '1.10E') ]
            p2 = [ format(plane[k, 0], '1.10E'), 
                   format(plane[k, 1], '1.10E'),
                   format(plane[k, 2], '1.10E'),
                   format((dist[k] + (2 * offset[k])), '1.10E') ]
            mcnp.write(str(98 - k) + '   P    '  
                        + '{: >17} {: >17} {: >17} {: >17}'.format(*p1) + '\n'
                        + '*' + str(98 - sets - k) + '  P    ')
            mcnp.write('{: >17} {: >17} {: >17} {: >17}'.format(*p2) + '\n')
    
    # Surface 99 is always output as a sphere that encompasses all UM geometry.
    s1 = [ format(x0, '1.10E'), 
           format(y0, '1.10E'),
           format(z0, '1.10E'),
           format(rad, '1.10E') ]
    mcnp.write('99   SPH  ' 
               + '{: >17} {: >17} {: >17} {: >17}'.format(*s1)
               + '\n\nc \n'
               + 'c DATA CARDS\n' 
               + 'c These cards link MCNP cell definitions to the' 
               + ' appropriate mesh file\n' 
               + 'c The number on embed must match the universe number' 
               + ' in the cell cards\n' 
               + 'embed1 meshgeo=abaqus\n' 
               + '       mgeoin=' + str(job) + '.inp\n' 
               + '       meeout=' + str(job) + '.eeout\n' 
               + '       length=1.0\n' 
               + '       background=' + str(blocks + 1) + '\n' 
               + 'c Set filetype=ASCII if human readable output is desired\n' 
               + '       filetype=binary\n' 
               + '       matcell=\n')
    
    # Setting material and pseudo-cell numbers.
    for i in range(blocks):
        mcnp.write('       ' + str(i + 1) + ' ' + str(i + 1) + '\n')
    
    # Setting up skeletal importance card.
    mcnp.write('c\nc \nc \n' 
               + 'c Default particle type is N\n' 
               + 'IMP:N 1 ')
    # Importances assume all cells are 1 except the outermost cell where 
    # particles terminate.
    if (reflect == True):
        mcnp.write(str(blocks + 2) + 'R 0\n')
    if (reflect == False):
        mcnp.write(str(blocks + 1) + 'R 0\n')
    
    mcnp.close()
    
# Set up and write the Abaqus inp file
def abaqus(rootgrp, blocks, len_node, el_types, mats, job, start_time, xyz):
    q = open(job + '.inp', 'w')
    q.write('*Heading\n' 
            + ' each element will become a different pseudo-cell\n' 
            + '** Job name: ' + str(job) + ' Model Name: ' + str(job) 
            + '\n' 
            + '** Generated by: Ex2Ab because Cubit does it WRONG\n'
            + '*Preprint, echo=NO, model=NO, history=NO, contact=NO\n' 
            + '**\n' 
            + '** PARTS\n' 
            + '**\n' 
            + '*Part, name=PART-DEFAULT' 
            + '\n' 
            + '*Node\n')
    
    print(str(len_node) + ' nodes detected, Runtime: ' 
          + str(timing(start_time)) + ' sec')
    
    # Write xyz node coordinate list to file.
    outList = []
    for n in range(len_node):
        line = (str(n + 1) + ' ' + '{0:.16e}'.format(xyz[n, 0]) + ',' 
                                    + '{0:.16e}'.format(xyz[n, 1])  + ',' 
                                    + '{0:.16e}'.format(xyz[n, 2])  + '\n')
        outList.append(line)
    q.writelines(outList)

    # List all elements and their nodes by block.
    L = 0
    k = 0
    elType = 'start'
    for i in range(blocks):
        outList = []
        t = np.array(rootgrp.variables['connect' + str(i + 1)])
        col = len(t[0])
        # Checks if element type has been written.
        # If not, the type code is written.
        if(elType != el_types[i]):
            q.write('*Element, type=C3D' + str(col) + '\n')
            # Only 1st and 2nd order tets, pents, and hexs are supported by 
            # MCNP. Cubit does not have pents like Abaqus does, so really only 
            # hexs and tets are supported if Cubit is used.
            elType = el_types[i]
        L = L + k
        k = len(t)
        if (blocks <= 100):
            print('Block ' + format((i + 1), '0>2d') 
                + ' has ' + format((k), '0>6d') + ' ' 
                + str(el_types[i]) + ' elements, Runtime: ' 
                + str(timing(start_time)) + ' sec')
        for j in range(k):
            #q.write(str(j + L + 1) + ',')
            # 2nd order hexs will be output on 2 lines of 15 and 5 nodes.
            if (col == 20):
            # Abaqus uses a different node order than exodus II does for HEX20 
            # elements. The order goes: vertices (the HEX8 nodes), then looping 
            # around faces listing the mid-edge nodes.
                line = ( str(j + L + 1) + ',' + str(t[j, 0]) + ','
                          + str(t[j, 1]) + ',' + str(t[j, 2]) + ','
                          + str(t[j, 3]) + ',' + str(t[j, 4]) + ','
                          + str(t[j, 5]) + ',' + str(t[j, 6]) + ','
                          + str(t[j, 7]) + ',' + str(t[j, 8]) + ','
                          + str(t[j, 9]) + ',' + str(t[j, 10]) + ','
                          + str(t[j, 11]) + ',' + str(t[j, 16]) + ','
                          + str(t[j, 17]) + ',' + str(t[j, 18]) + ',\n     '
                          + str(t[j, 19]) + ',' + str(t[j, 12]) + ','
                          + str(t[j, 13]) + ',' + str(t[j, 14]) + ','
                          + str(t[j, 15]) + '\n')
            # All other elements list all nodes on 1 line
            if (col != 20):
                #t[j, :].tofile(q, sep = ',', format = '%i')
                #q.write('\n')
                line = str(j + L + 1) + ','
                for c in range(col):
                    if (c != col - 1):
                        line = line + str(t[j, c]) + ','
                    else:
                        line = line + str(t[j, c]) + '\n'
            outList.append(line)
        q.writelines(outList)
                
    print(str(L + k) + ' elements detected, Runtime: ' 
          + str(timing(start_time)) + ' sec')
    
    # List node sets. Hard-coded as only 1 set. Not sure for when this would 
    # change
    # since models should be a single part. This doesn't count node sets used 
    # for reflecting boundaries.
    q.write('*Nset, nset=ALLNODES, generate\n' + '1,' + str(len_node) 
            + ',1\n')
    L = 0
    k = 0
    
    # Some basic meshes created with libMesh did not have the 'elem_map' field.
    # This checks if it's there. If not, sequential numbers are assumed.
    el_map = True
    try: el_num = np.array(rootgrp.variables['elem_map'])
    except: el_map = False
    
    # List element sets.
    outList = []
    for i in range(blocks):
        #outList = []
        t = np.array(rootgrp.variables['connect' + str(i + 1)])
        L = L + k
        k = len(t)
        #q.write('*Elset, elset=SET_MATERIAL_TALLY_' + str(i + 1) + '\n')
        outList.append('*Elset, elset=SET_MATERIAL_TALLY_' + str(i + 1) + '\n')
        for j in range(k):
            if (el_map == True):
                #q.write(str(el_num[j + L]) + ',')
                #line = line + str(el_num[j + L]) + ','
                outList.append(str(el_num[j + L]) + ',')
            if (el_map == False):
                #q.write(str(j + L + 1) + ',')
                #line = line + str(j + L + 1) + ','
                outList.append(str(j + L + 1) + ',')
            if (((j + 1) % 16 == 0) and ((j + 1) != k)):
                #q.write('\n')
                #line = line + '\n'
                outList.append('\n')
        outList.append('\n')
    q.writelines(outList)
    
    # List sections of the model.
    for i in range(blocks):
        t = rootgrp.variables['connect' + str(i + 1)]
        q.write('** Section: Section-' + str(i + 1) + '-SET_MATERIAL_TALLY_'
                + str(i + 1) + '\n'
                + '*Solid Section, elset=SET_MATERIAL_TALLY_' + str(i + 1) 
                + ', material=Material_block_' + str(i + 1) + '\n,\n')
    
    # Write assembly information.
    q.write('*End Part\n**\n**\n'
            + '** ASSEMBLY\n**\n'
            + '*Assembly, name=Assembly\n**\n' 
            + '*Instance, name=PART-DEFAULT_1, part=PART-DEFAULT\n' 
            + '*End Instance\n**\n*End Assembly\n**\n')
    
    # List materials and densities.
    q.write('** MATERIALS\n**\n')
    for i in range(blocks):
        t = rootgrp.variables['connect' + str(i + 1)]
        q.write('*Material, name=Material_block_' + str(i + 1) 
                + '\n*Density\n' 
                + str(mats[i]) + ',\n')
    
    # List other information.
    # There are some hardcoded numbers here. Does not affect MCNP. May only 
    # matter for Abaqus.
    q.write('** ------------------------------'
            + '----------------------------------\n**\n'
            + '** STEP: DefaultSet\n**\n'
            + '*Step, name=DefaultSet, nlgeom=NO\n' 
            + '*Static\n1., 1., 1e-05, 1.\n**\n' 
            + '** OUTPUT REQUESTS\n**\n'
            + '*Restart, write, frequency=0\n**\n' 
            + '** FIELD OUTPUT: F-Output-1\n**\n' 
            + '*Output, field, variable=PRESELECT\n**\n'
            + '** HISTORY OUTPUT: H-Output-1\n**\n'
            + '*Output, history, variable=PRESELECT\n*End Step')
    
    q.close()