abaqusToOpenFOAM.py

'''
Run like:

abaqus cae noGUI=abaqusToOpenFOAM.py -- path/to/cae/model model-name
# So the polymesh will be upper diagonal
renumberMesh -overwrite
# Make sure all is OK
checkMesh
# If model built in mm
transformPoints -scale '(1e-3 1e-3 1e-3)'
'''
from collections import OrderedDict
import os
import errno
import sys

from abaqus import *
from abaqusConstants import *
from caeModules import *

'''
The owner and neighbour files each contain a list of cell numbers. The position in that list is the face number. For example:

(
0 // this is face 0
0 // this is face 1
0 // this is face 2
1 // this is face 3
1 // this is face 4
... etc.

So, face 2's "owner" is cell 0, face 4's "owner" is cell 1, and so on.
Note that the owner file is longer than the neighbour file... that's because the owner file also lists boundary faces, whereas the neighbour file does not.
'''

nodes = OrderedDict()
elements = {}
faces = OrderedDict()

def mkdir_p(path):
    try:
        os.makedirs(path)
    except OSError as exc:  # Python >2.5
        if exc.errno == errno.EEXIST and os.path.isdir(path):
            pass
        else:
            raise

PATH = 'constant/polyMesh'
mkdir_p(PATH)

face_connectivity = {FC3D4: {FACE1:[0, 2, 1],
                             FACE2:[0, 1, 3],
                             FACE3:[1, 2, 3],
                             FACE4:[0, 3, 2]},
                     FC3D6: {FACE1:[0, 2, 1],
                             FACE2:[3, 4, 5],
                             FACE3:[0, 1, 4, 3],
                             FACE4:[1, 2, 5, 4],
                             FACE5:[0, 3, 5, 2]}}

FOOTER = '// ************************************************************************* //'

class Node(object):
    def __init__(self, index, node):
        self.index = index
        self.coordinates = node.coordinates
        self.label = node.label
        self.instanceName = node.instanceName

    def __repr__(self):
        return '(Node label {}, index {})'.format(self.label, self.index)

    def __hash__(self):
        return hash(repr(self.index))


class Element(object):
    def __init__(self, element, index):
        self.cell_index = index
        self.label = element.label
        self.type = element.type
        self.instanceName = element.instanceName
        self.connectivity = element.connectivity
        self.getElemFaces = element.getElemFaces

    def __repr__(self):
        return '(Element label {}, cell index {})'.format(self.label, self.instanceName)

    def __hash__(self):
        return hash((self.label, self.instanceName))


class Boundary(object):
    def __init__(self, name, wall):
        self.name = name
        self.wall = wall


class Face(object):
    def __init__(self, face):
        self.index = None
        self.face_id = face.face
        self.nodes = face.getNodes()
        self.neighbour = None
        for element in face.getElements():
            if element.label != face.label:
                self.neighbour = element
            else:
                self.owner = element
        self.boundary = Boundary('default', False)


def write_points(nodes):
    header = r'''/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  4.0                                   |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       vectorField;
    location    "constant/polyMesh";
    object      points;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
    '''
    with open(os.path.join(PATH,'points'), 'w') as f:
        f.write(header)
        f.write('{}\n'.format(len(nodes)))
        f.write('(\n')
        for node in nodes.itervalues():
            f.write('// Node label {}\n'.format(node.label))
            f.write('({} {} {})\n'.format(*node.coordinates))
        f.write(')\n')
        f.write(FOOTER + '\n')


def write_faces(faces):
    header = r'''/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  4.0                                   |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       faceList;
    location    "constant/polyMesh";
    object      faces;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
'''
    with open(os.path.join(PATH,'faces'), 'w') as f:
        f.write(header)
        f.write('{}\n'.format(len(faces)))
        f.write('(\n')
        for face in faces.itervalues():
            f.write('{}({}) // Element label {}\n'.format(len(face.nodes), 
                                                            " ".join([str(nodes[(n.label, n.instanceName)].index) for n in face.nodes]), 
                                                            face.owner.label))
        f.write(')\n')
        f.write(FOOTER + '\n')


def write_owner(faces):
    header = r'''/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  4.0                                   |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       labelList;
    location    "constant/polyMesh";
    object      owner;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
'''
    with open(os.path.join(PATH,'owner'), 'w') as f:
        num_neighbour = 0
        f.write(header)
        f.write('{}\n'.format(len(faces)))
        f.write('(\n')
        for face in faces.itervalues():
            if face.neighbour is not None:
                num_neighbour = num_neighbour + 1
            f.write('{}\n'.format(elements[(face.owner.label, face.owner.instanceName)].cell_index))
        f.write(')\n')
        f.write(FOOTER + '\n')
        return num_neighbour


def write_neighbours(num_neighbours, faces):
    header = r'''/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  4.0                                   |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       labelList;
    location    "constant/polyMesh";
    object      neighbour;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
'''
    with open(os.path.join(PATH,'neighbour'), 'w') as f:
        f.write(header)
        f.write('{}\n'.format(num_neighbours))
        f.write('(\n')
        for face in faces.itervalues():
            if face.neighbour is not None:
                f.write('{}\n'.format(elements[(face.neighbour.label, face.neighbour.instanceName)].cell_index))
        f.write(')\n')
        f.write(FOOTER + '\n')


def write_boundary(faces):
    header = r'''/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  4.0                                   |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       polyBoundaryMesh;
    location    "constant/polyMesh";
    object      boundary;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
'''
    boundaries = {}
    with open(os.path.join(PATH,'boundary'), 'w') as f:
        f.write(header)
        for face in faces.itervalues():
            if face.boundary.name != 'default':
                boundaries.setdefault(face.boundary.name, {'start': face.index,
                                                           'wall': face.boundary.wall}).update({'end':face.index})
        f.write('{}\n'.format(len(boundaries)))
        f.write('(\n')
        for name, d in sorted(boundaries.iteritems(), key=lambda x: x[1]['start']):
            f.write('{}\n'.format(name))
            f.write('{\n')
            if d['wall']:
                f.write('type wall;\n')
            else:
                f.write('type patch;\n')
            f.write('nFaces {};\n'.format(d['end'] - d['start'] + 1))
            f.write('startFace {};\n'.format(d['start']))
            f.write('}\n')
 
        f.write(')\n')
        f.write(FOOTER + '\n')


def get_nodes_elements(model):
    node_index = 0
    cell_index = 0
    for instance in model.rootAssembly.allInstances.values():
        print('Instance: {}'.format(instance.name))
        # Get all of the nodes
        print('Number of nodes: {}'.format(len(instance.nodes)))
        for node in instance.nodes:
            nodes[(node.label, node.instanceName)] = Node(node_index, node)
            node_index = node_index + 1
        print('Number of elements: {}'.format(len(instance.elements)))
        # Get all of the elements
        for element in instance.elements:
            elements[(element.label, element.instanceName)] = Element(element, cell_index)
            cell_index = cell_index + 1


def extract_faces(model):
    print('Get all faces')
    _faces = OrderedDict()
    for instance in model.rootAssembly.allInstances.values():
        for element in instance.elements:
            for face in element.getElemFaces():
                f = Face(face)
                key = frozenset([nodes[(n.label, n.instanceName)].index for n in f.nodes])
                if key not in _faces:
                    _faces[key] = f
        

    print('From boundary conditions')
    for bc_name in model.boundaryConditions.keys():
        print('Extracting faces for bc {}'.format(bc_name))
        bc = model.boundaryConditions[bc_name]
        if "<type 'FluidWallConditionBC'>" == str(bc.__class__):
            wall = True
        else:
            wall = False
        boundary_condition = Boundary(bc_name, wall)
        region = bc.region
        '''
        Don't quite know the meaning of the first two flags, but the last flag is a boolean 
        indicating whether or not the set is an internal set or not.
        '''
        if region[1] == 'Assembly':
            if region[4]:
                # Internal set
                surface = model.rootAssembly.allInternalSurfaces[region[0]]
            else:
                # Regular set
                surface = model.rootAssembly.allSurfaces[region[0]]
        else:
            raise ValueError('Unhandled owner region')
        for surface_face in surface.faces:
            for face in surface_face.getElementFaces():
                f = Face(face)
                key = frozenset([nodes[(n.label, n.instanceName)].index for n in f.nodes])
                if key not in _faces:
                    # This should not happen
                    raise
                else:
                    _faces[key].boundary = boundary_condition

    print('Get internal faces (if neighbour, then internal face)')
    face_index = 0
    for key, f in _faces.iteritems():
        if f.neighbour is not None:
            f.index = face_index
            faces[key] = f
            face_index = face_index + 1

    print('Sort by boundary condition name')
    for key, f in sorted(_faces.iteritems(), key=lambda x: x[1].boundary.name):
        if f.boundary.name != 'default':
            f.index = face_index
            faces[key] = f
            face_index = face_index + 1



def main(db_name, model_name):
    openMdb(pathName=db_name)
    model = mdb.models[model_name]

    get_nodes_elements(model)

    print('Writing points ...')
    write_points(nodes)

    print('Extracting faces ...')
    extract_faces(model)

    print('Writing faces ...')
    write_faces(faces)

    print('Writing owner ...')
    num_neighbours = write_owner(faces)

    print('Writing neighbours ...')
    write_neighbours(num_neighbours, faces)

    print('Writing boundary ...')
    write_boundary(faces)

if __name__ == '__main__':
    print('{}'.format(sys.argv))
    db_name = sys.argv[8]
    model_name = sys.argv[9]
    main(db_name, model_name)