import sys
import time
import os
import inspect
try: from lxml import etree as et # not in my standard Salome's python but better
except ImportError: import xml.etree.cElementTree as et # functionally equavalent, but no pretty printing
try:
import salome
import salome_notebook
from salome.smesh import smeshBuilder
import SMESH
except ImportError:
print('\n\n\nUsage: salome {} args:ConnectivityFile,FeBioFile(optional),GeometryFile(optional),steps(optional) \n\n\n'.format( inspect.getfile(inspect.currentframe()) ))
sys.exit()
try:from MEDLoader import MEDLoader
except ImportError: import MEDLoader
import StlToMed
salome.salome_init()
smesh = smeshBuilder.New()
def myElement(tag, parent=None, text=None, attrib={}, **extra):
e = et.Element(tag, attrib, **extra)
if text: e.text = text
if parent is not None: # because "if parent" doesn't work as I think is should.
parent.append(e)
return e
def CloseSalome():
"""Try to close out of Salome gracefully"""
# there is nothing graceful about this
try:
sys.stdout = 'redirect to nowhere...' # so I don't see the junk it prints
from killSalomeWithPort import killMyPort
killMyPort(os.getenv('NSPORT'))
except: pass
def MaterialsXml(parent, materials, rigids):
"""Temporary cludge!
This function not fully implemented
User should manually fix material assignments for now!"""
print('Making Material')
Material = myElement('Material', parent)
#
for i, mat in enumerate(rigids):
material = myElement('material', Material, id=repr(i+1), name=mat, type='rigid body')
myElement('density', material, '1e-9')
for i, mat in enumerate(materials):
material = myElement('material', Material, id=repr(i+1+len(rigids)), name=mat, type="Mooney-Rivlin")
myElement('c1', material, '0.1')
myElement('c2', material, '0')
myElement('k', material, '100')
myElement('density', material, '1e-9')
def GeometryXml(parent, Assembly, materials, rigids):
"""Create the Geometry SubElement to parent
consists of Nodes, NodeSets, Surfaces, SurfacesPairs, and Elements.
re-looping allows for necessary placement of Geometry subsections,"""
print('Making Geometry')
geom = myElement('Geometry', parent)
#
print(' Making Nodes')
pN = 0 # pN for previousNodes
for part in Assembly.values():
nodes = myElement('Nodes', geom, name=part.name)
partCoords=MEDLoader.ReadUMeshFromFile(part.med, 0).getCoords()
for i in range(len(partCoords)): myElement('node', nodes, repr( partCoords.getTuple(i) )[1:-1], id=repr(i+1 + pN))
pN += len(partCoords)
#
print(' Making NodeSets')
pN = 0 # reset node numbering offset
for part in Assembly.values(): # iterate again for NodeSet
groups = part.mesh.GetGroups()
for group in groups: # assume triangle and tetrahedra, if bad assumption mesh.GetElementShape or mesh.GetElementGeomType might work
name = group.GetName()
if name.endswith('Nodes'):
nodeSet = myElement('NodeSet', geom, name=name)
for node in sorted(group.GetListOfID()): myElement('node', nodeSet, id=repr(node + pN))
pN += part.mesh.NbNodes()
#
print(' Making Surfaces')
pN = 0 # reset node numbering offset, but pN changed to 1 due to use of MEDCoupling
for part in Assembly.values(): # iterate again for Surfaces
dim = part.mesh.MeshDimension()
if dim == 3: medMesh = MEDLoader.ReadUMeshFromFile(part.med,-1)
elif dim == 2: medMesh = MEDLoader.ReadUMeshFromFile(part.med, 0)
groups = part.mesh.GetGroups()
for group in groups:
name = group.GetName()
print(' {}'.format(name))
if name.endswith('Faces'):
surf = myElement('Surface', geom, name=name)
for elem in sorted(group.GetListOfID()):
# index = (elem-1)*(nodeCount+1) # MEDCoupling element start at 0 not 1
# connectivity = repr( [meshConnectivity[index+c] + pN for c in range(1,nodeCount+1)] )[1:-1]
connectivity = repr( [node+pN for node in part.mesh.GetElemNodes(elem)] )[1:-1]
# if triangle bad assumption, type = medMesh.GetReprOfGeometricType(*medMesh.getAllGeoTypes()) should work
type = medMesh.getAllGeoTypes()[0]
if type == 3: type = 'tri3'
elif type == 6: type = 'tri6'
elif type == 4: type = 'quad4'
elif type == 8: type = 'quad8'
else: print('type {} not supported, check Febio documentation.'.format(medMesh.GetReprOfGeometricType()))
myElement(type, surf, connectivity, id=repr(elem))
pN += part.mesh.NbNodes()
#
# Surface Pairs MUST come after surfaces are defined!
#
print(' Making SurfacesPairs')
for master, part in Assembly.items(): # iterate again for SurfacePairs
for slave in part.Contact:
pair = myElement('SurfacePair', geom, name='{}_To_{}'.format(master, slave))
myElement('master', pair, surface='{}_@_{}_ContactFaces'.format(master, slave))
myElement('slave', pair, surface='{}_@_{}_ContactFaces'.format(slave, master))
if master in materials: # surface pairs are also needed for deformable tied surfaces
for slave in part.Ties:
if slave in materials:
pair = myElement('SurfacePair', geom, name='{}_To_{}'.format(master, slave))
myElement('master', pair, surface='{}_@_{}_TiesFaces'.format(master, slave))
myElement('slave', pair, surface='{}_@_{}_TiesFaces'.format(slave, master))
else: pass # this is handeled as a boundary condition instead
#
print(' Making Elements')
pN = 1 # reset node numbering offset, but pN changed to 1 due to use of MEDCoupling
eN = 1 #
for part in Assembly.values(): # iterate one last time for Elements
medMesh = MEDLoader.ReadUMeshFromFile(part.med, 0)
meshConnectivity = medMesh.getNodalConnectivity()
# dim = medMesh.getMeshDimension()
#
name = part.name
if name in rigids: num = str(rigids.index(name) + 1)
elif name in materials: num = str(materials.index(name) + 1 + len(rigids))
else: print(name)
#
type = medMesh.getAllGeoTypes()[0] # this assumes that meshes are not mixed types
if type == 3: vol = myElement('Elements', geom, name=name, type='tri3', mat=num)
elif type == 6: vol = myElement('Elements', geom, name=name, type='tri6', mat=num)
elif type == 4: vol = myElement('Elements', geom, name=name, type='quad4', mat=num)
elif type == 8: vol = myElement('Elements', geom, name=name, type='quad8', mat=num)
elif type == 14: vol = myElement('Elements', geom, name=name, type='tet4', mat=num)
elif type == 20: vol = myElement('Elements', geom, name=name, type='tet10', mat=num)
elif type == 16: vol = myElement('Elements', geom, name=name, type='penta6', mat=num)
elif type == 25: vol = myElement('Elements', geom, name=name, type='penta15',mat=num)
elif type == 18: vol = myElement('Elements', geom, name=name, type='hex8', mat=num)
elif type == 30: vol = myElement('Elements', geom, name=name, type='hex20', mat=num)
# elif type == 27: vol = myElement('Elements', geom, name=name, type='hex27', mat=num)
else: print('type {} not supported check Febio documentation.'.format(medMesh.GetReprOfGeometricType()))
#
nodeCount = medMesh.getNumberOfNodesInCell(0)
numElems = medMesh.getNumberOfCells()
for elem in range(numElems):
index = elem * (nodeCount+1)
connectivity = [meshConnectivity[index+c] + pN for c in range(1,nodeCount+1)]
if type == 3: pass # tri3
elif type == 6: pass # tri6
elif type == 4: pass # quad4
elif type == 8: pass # quad8
elif type == 14: #tet4
connectivity[1], connectivity[2] = connectivity[2], connectivity[1] # swap nodes n2 and n3
elif type == 20: # tet10
connectivity[1], connectivity[2] = connectivity[2], connectivity[1] # swap nodes n2 and n3
connectivity[4], connectivity[6] = connectivity[6], connectivity[4] # swap nodes n5 and n7
connectivity[8], connectivity[9] = connectivity[9], connectivity[8] # swap nodes n9 and n10
elif type == 16: # penta6
connectivity[1], connectivity[2] = connectivity[2], connectivity[1] # swap nodes n2 and n3
connectivity[4], connectivity[5] = connectivity[5], connectivity[4] # swap nodes n5 and n6
elif type == 25: # penta15
connectivity[1] , connectivity[2] = connectivity[2] , connectivity[1] # swap nodes n2 and n3
connectivity[4] , connectivity[5] = connectivity[5] , connectivity[4] # swap nodes n5 and n6
connectivity[6] , connectivity[8] = connectivity[8] , connectivity[6] # swap nodes n5 and n6
connectivity[9] , connectivity[11] = connectivity[11], connectivity[9] # swap nodes n5 and n6
connectivity[13], connectivity[14] = connectivity[14], connectivity[13] # swap nodes n5 and n6
elif type == 18: # hex8
connectivity[1], connectivity[3] = connectivity[3], connectivity[1] # swap nodes n2 and n4
connectivity[5], connectivity[7] = connectivity[7], connectivity[5] # swap nodes n6 and n8
elif type == 30: # hex20
connectivity[1], connectivity[3] = connectivity[3], connectivity[1] # swap nodes n2 and n4
connectivity[5], connectivity[7] = connectivity[7], connectivity[5] # swap nodes n6 and n8
connectivity[8], connectivity[11]= connectivity[11],connectivity[8] # swap nodes n9 and n12
connectivity[9], connectivity[10]= connectivity[10],connectivity[9] # swap nodes n10 and n11
connectivity[12],connectivity[15]= connectivity[15],connectivity[12] # swap nodes n13 and n16
connectivity[13],connectivity[14]= connectivity[14],connectivity[13] # swap nodes n14 and n15
connectivity[17],connectivity[19]= connectivity[19],connectivity[17] # swap nodes n18 and n20
# elif type == 27: # hex27, I don't have information on this type, but it is included in both MED and Febio
myElement('elem', vol, repr(connectivity)[1:-1], id=repr(elem + eN))
eN += numElems
pN += len(medMesh.getCoords())
def BoundariesXml(parent, Assembly, rigids):
"""Create the Boundary SubElement to parent
consists of attachments between Tied rigid and deformable bodies
and fixed all DoF for each rigid body
user should manually change DoF as needed."""
print('Making Boundary')
Boundary = myElement('Boundary', parent)
#
# Attach deformable bodies to rigid bodies
#
for master, part in Assembly.items():
for slave in part.Ties:
print(master, slave)
if slave in rigids:
name = '{}_With_{}'.format(master, slave)
set = '{}_@_{}_TiesNodes'.format(master,slave)
id = str(rigids.index(slave)+1)
myElement('rigid', Boundary, name=name, rb=id, node_set=set)
#
# # Make dangling (vector) boundaries
#
for master, part in Assembly.items():
for item in part.Other:
if item == 'Vector':
myElement('fix', Boundary, bc="x,y,z", node_set='{}_Vector_Nodes'.format(master,slave))
#
# Fix each rigid body in place, user can modify from here
#
DoF={'x':'fixed',
'y':'fixed',
'z':'fixed',
'Rx':'fixed',
'Ry':'fixed',
'Rz':'fixed'} # 'z':'force', "z":"prescribed"
for r in rigids:
id = str(rigids.index(r)+1)
rigid = myElement('rigid_body', Boundary, mat=id)
for d, type in DoF.items():
a = myElement(type, rigid, bc=d)
if type != 'fixed':
a.set('lc', '1')
a.text = '1.0'
#
# Comment for copy pasting other options
#
Boundary.append( et.Comment('<prescribed bc="z" lc="1">1.0</prescribed>') )
Boundary.append( et.Comment('<force bc="z" lc="1">10.0</force>') )
def ContactXml(parent, Assembly, materials):
"""Create the Contact SubElement to parent
consists of facet-tofacet sliding for each contact pair"""
print('Making Contact')
Contact = myElement('Contact', parent)
#
# Remove to prevent doubling contact i.e. a_To_b b_To_a. ContactXml is the last use of Assembly so no need to make a copy first.
for part in Assembly.values():
for contactName in part.Contact: del Assembly[contactName].Contact[part.name]
for tiedName in part.Ties: del Assembly[tiedName].Ties[part.name]
#
contactParams={'laugon':"1",'penalty':"1", 'auto_penalty':"1",'maxaug':"10"}
for master, part in Assembly.items():
for slave in part.Contact:
c = myElement('contact', Contact, type='facet-to-facet sliding', surface_pair='{}_To_{}'.format(master, slave))
for p in contactParams:
myElement(p, c, text=contactParams[p])
#
contactParams = {'penalty': "10"} # maybe tolerance is required , 'tolerance': "1.0"}
for master, part in Assembly.items():
if master in materials: # Tied contact is only for deformable to deformable surfaces
for slave in part.Ties:
if slave in materials:
c = myElement('contact', Contact, type='tied', surface_pair='{}_To_{}'.format(master, slave))
for p in contactParams:
myElement(p, c, text=contactParams[p])
def LoadDataXml(parent, steps=10):
"""Create the LoadData SubElement to parent
this assumes it will have a single smooth load id=1
and a step up function id=2"""
load = myElement('LoadData', parent)
singleStep = myElement('loadcurve', load, attrib={"id":"1", "type":"linear"})
myElement('point', singleStep, "0,0")
myElement('point', singleStep, "1,1")
manyStep = myElement('loadcurve', load, attrib={"id": "2", "type":"step"})
for i in range(steps + 1):
myElement('point', manyStep, "{:4.2f},{}".format( float(i)/steps, 1.0/steps ))
def OutputXml(parent, outputData=[], logFile=False, logData=None):
"""Create the LoadData SubElement to parent
this assumes it will have a single smooth load id=1
and a step up function id=2"""
out = myElement('Output', parent)
plot = myElement('plotfile', out, type='febio')
requiredData = ["contact gap","contact pressure","contact traction","displacement","reaction forces","stress"]
for v in requiredData: myElement('var', plot, type=v)
for v in outputData: myElement('var', plot, type=v)
plot.append(et.Comment("""<var type="prestrain stretch"/>"""))
#
if logFile:
log = myElement('logfile', out, file=logFile)
for data in logData: myElement(data, plot)
def StepXml(parent, name, steps):
step = myElement('Step', parent, name=name)
control = myElement('Control', step, type='febio')
#
myElement('time_steps', control, text=str(steps) )
myElement('step_size', control, text=str(1./steps))
#
stepper = myElement('time_stepper', control)
myElement('dtmin', stepper, text='1e-10')
myElement('dtmax', stepper, text=str(1./steps), lc="2")
myElement('max_retries', stepper, text='20')
myElement('opt_iter', stepper, text='10')
myElement('aggressiveness', stepper, text='1')
#
myElement('optimize_bw', control, text='1')
myElement('plot_level', control, text='PLOT_MUST_POINTS')
myElement('analysis', control, type='static')
myElement('min_residual', control, text='1e-10')
#
if name == 'Initialize':
comment = et.Comment("""<Constraints>\n\t\t\t<constraint type="prestrain">\n\t\t\t\t<tolerance>0.01</tolerance>\n\t\t\t\t<update>1</update>\n\t\t\t</constraint>\n\t\t</Constraints>""")
step.append(comment)
def ReadConnectivity(connectivityFile):
"""Read the connectivity file and asserts that all files it references exist."""
import ConnectivityXml as xml
print('\n') # just make some white space for easier reading of the output
print('Reading from connectivity file {}'.format(connectivityFile))
if os.path.dirname(connectivityFile) != os.getcwd() and os.path.dirname(connectivityFile):
os.chdir(os.path.dirname(connectivityFile)) # connectivity file may not give full file path information.
#
try:
Assembly = xml.ReadConnectivity(connectivityFile) # make sure
# xml.RemoveDoubleCounting(Assembly)
except IOError:
print('\n\nNo connectivity file given.\nUsage: salome MedToFebio.py args:ConnectivityFile,FeBioFile(optional),GeometryFile(optional),steps(optional)\n')
sys.exit()
#
for part in Assembly.values():
print(part.file)
if part.file.endswith('.med'): file = part.file
else: file = os.path.join(os.path.dirname(connectivityFile),'MED',os.path.basename(part.file).replace('.stl','.med'))
part.mesh,=StlToMed.OpenMesh(file) # comma necessary to unpack the list (of one) returned by opening MED files
print('')
#
return Assembly
def WriteFebioMainFile(febioFile, geoFile, Assembly, steps=10):
"""Write the main febio file
this will import the geometry section from a seperate file due to its size."""
febio = et.Element('febio_spec', {'version': "2.5"})
tree = et.ElementTree(febio)
#
comment = et.Comment('This document was written by\n {}\n on {}'.format(inspect.getfile(inspect.currentframe()),
time.strftime('%b %d %Y at %H:%m:%S %p')))
febio.append(comment)
#
myElement('Module' , febio, type='solid')
#
matList = []; rigidList = []
for name, part in Assembly.items():
if part.mat.lower() == 'rigid': rigidList.append(name)
else: matList.append(name)
#
MaterialsXml(febio, matList, rigidList)
myElement('Geometry', febio, attrib={'from': geoFile}) # as a dictionary because from is a keyword
BoundariesXml(febio, Assembly, rigidList)
ContactXml(febio, Assembly, matList)
#
LoadDataXml(febio, steps=steps)
OutputXml(febio)
StepXml(febio, 'Initialize', steps)
StepXml(febio, 'LoadingStep', steps)
#
if not os.path.dirname(febioFile):
febioFile = os.path.join(os.getcwd(), 'Febio', febioFile)
if not os.path.isdir('Febio'): os.mkdir('Febio')
print('Writing file {}'.format(febioFile))
try: tree.write(febioFile, xml_declaration=True, pretty_print=True)
except TypeError:
f = open(febioFile,'w')
f.write('<?xml version="1.0" encoding="ISO-8859-1"?>\n')
f.write(et.tostring(febio).decode().replace('>','>\n'))
f.close()
def WriteFebioGeometryFile(geoFile, Assembly):
febio = et.Element('febio_spec', {'version': "2.5"})
tree = et.ElementTree(febio)
#
comment = et.Comment('This document was written by\n {}\n on {}'.format(inspect.getfile(inspect.currentframe()),
time.strftime('%b %d %Y at %H:%m:%S %p')))
febio.append(comment)
#
matList = []; rigidList = []
for name, part in Assembly.items():
if part.mat.lower() == 'rigid': rigidList.append(name)
else: matList.append(name)
#
GeometryXml(febio, Assembly, matList, rigidList)
#
if not os.path.dirname(geoFile):
geoFile = os.path.join(os.getcwd(), 'Febio', geoFile)
if not os.path.isdir('Febio'): os.mkdir('Febio')
print('Writing file {}'.format(geoFile))
try: tree.write(geoFile, xml_declaration=True, pretty_print=True)
except TypeError:
f = open(geoFile,'w')
f.write('<?xml version="1.0" encoding="ISO-8859-1"?>\n')
f.write(et.tostring(febio).decode().replace('>','>\n'))
f.close()
def MakeFebio(connectivityFile='Connectivity.xml', febioFile='FeBio.feb', geoFile='Geometry.feb', steps=10):
"""Read connectivity file, get MEDs and write Febio input file"""
Assembly = ReadConnectivity(connectivityFile)
#
WriteFebioGeometryFile(geoFile, Assembly)
WriteFebioMainFile(febioFile, geoFile, Assembly, steps)
#
# raw_input('...') # soetimes I want to stop before closing when testing
CloseSalome()
if __name__ == '__main__':
MakeFebio(*sys.argv[1:])