import sys
import time
import os
try:
import salome
import salome_notebook
import SMESH
from salome.smesh import smeshBuilder
except ImportError:
import inspect
print '\n\n\nUsage: salome {} args:Sofa.scn,ConnectivityFile\n\n\n'.format( inspect.getfile(inspect.currentframe()) )
sys.exit()
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
import StlToMed
# the built-in Element constructor doesn't allow including text as a non-tag item, so this function does
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(2810)
except: pass
def BasicScene(parent, background=None):
"""Define all the "standard" sofa elements
that will remain unchanged by mesh changes"""
myElement('VisualStyle', parent, displayFlags='showVisualModels')
if background: myElement('BackgroundSetting', parent, image=background)
myElement('DefaultPipeline', parent, name="DefaultCollisionPipeline",verbose="0",draw="0",depth="100")
myElement('BruteForceDetection', parent, name = "Detection")
myElement('MinProximityIntersection', parent, name="Proximity",alarmDistance="0.03",contactDistance="0.01")
myElement('DefaultContactManager', parent, name="Response", response="default")
myElement('DefaultCollisionGroupManager', parent, name="Group")
myElement('MinProximityIntersection', parent, alarmDistance="8", contactDistance="5") # This one might actually need to be changed.
myElement('EulerImplicit', parent, name="cg_odesolver", printLog="false")
myElement('CGLinearSolver', parent, iterations="25",name="linear solver",tolerance="1.0e-5",threshold="1.0e-5")
def SofaObjectNode(parent, part):
"""Creates the xml Node for an object in the Sofa display"""
name = part.name
object = myElement('Node', parent, name=name)
if part.mat == 'rigid':
myElement('MechanicalObject', object, template="Rigid3d", name="{}Object".format(name), translation="0 0 0", rotation="0 0 0")
movementComment = et.Comment('<LinearMovementConstraint template="Rigid3d" name="Displacement" keyTimes="t1 t2 t3" movements="x1 y1 z1 rx1 ry1 rz1 x2 y2 z2 rx2 ry2 rz2 x3 y3 z3 rx3 ry3 rz3"/>')
object.append(movementComment)
visu=myElement('Node', object, name="Visu")
myElement('MeshGmshLoader', visu, filename='{}.msh'.format(name), name="loader")
myElement('TriangleSetTopologyContainer', visu, name="Container", src="@loader")
if 'Probe' in part.name: materialColor = "Default Diffuse 1 0.7 0.7 0.7 1 Ambient 0 0.8 0.8 0.8 1 Specular 0 1 0 0 1 Emissive 0 1 0 0 1 Shininess 0 45" # dull grey
else: materialColor = "Default Diffuse 1 0.7 0.7 0.7 1 Ambient 1 0.8 0.8 0.8 1 Specular 0 1 0 0 1 Emissive 0 1 0 0 1 Shininess 0 45" # bone
myElement('VisualModel', visu, template="ExtVec3f", name="Visual", material=materialColor)
myElement('RigidMapping', visu, input="@..", output="@Visual")
for tie in part.Ties:
t=myElement('Node', object, name=tie)
myElement('MeshGmshLoader', t, filename='{}_@_{}.msh'.format(tie,name), name="loader")
myElement('MechanicalObject', t, src="@loader", name="Attach")
myElement('RigidMapping', t, input="@..",output="@Attach")
for i, contact in enumerate(part.Contact):
c = myElement('Node', object, name='{}Contact'.format(contact) )
myElement('MeshGmshLoader', c, filename='{}_@_{}.msh'.format(contact, name), name="loader")
myElement('Mesh', c, src="@loader")
myElement('MechanicalObject', c, src="@loader", name="Attach")
myElement('RigidMapping', c, input="@..", output="@Attach")
myElement('TriangleModel', c, name="Triangles For Collision", contactStiffness="100000", group=str(i))
myElement('LineModel', c, name="Lines For Collision", contactStiffness="100000", group=str(i))
myElement('PointModel', c, name="Points For Collision", contactStiffness="100000", group=str(i))
else: # part.mat == 'elastic'
myElement('MeshGmshLoader', object, filename='{}.msh'.format(name), name="loader", translation="0 0 0", rotation="0 0 0")
myElement('MechanicalObject', object, src="@loader", name="Volume")
myElement('include', object, href="Objects/TetrahedronSetTopology.xml", src="@loader")
myElement('UniformMass', object, totalmass="100") #hmmm
myElement('TetrahedralCorotationalFEMForceField', object, youngModulus="300", poissonRatio="0.4" ,method="large") # more logic needed!
surf = myElement('Node', object, name="Surface")
myElement('TriangleSetTopologyContainer', surf, name="Container")
myElement('TriangleSetTopologyModifier', surf, name="Modifier")
myElement('TriangleSetTopologyAlgorithms', surf, name="TopoAlgo", template="Vec3d")
myElement('TriangleSetGeometryAlgorithms', surf, name="GeomAlgo", template="Vec3d")
myElement('Tetra2TriangleTopologicalMapping', surf, input="@../Container", output="@Container")
myElement('TriangleSet', surf, contactStiffness='3000')
myElement('LineModel', surf, contactStiffness='3000')
myElement('PointModel', surf, contactStiffness='3000')
visu=myElement('Node', surf ,name="Visu")
if 'Cartilage' in part.name: materialColor = "Default Diffuse 1 0.6 0.8 1.0 1 Ambient 1 0.6 0.8 1.0 1 Specular 0 1 0 0 1 Emissive 0 1 0 0 1 Shininess 0 10" # light blue
elif 'Meniscus' in part.name: materialColor = "Default Diffuse 1 0.6 1.0 0.6 1 Ambient 1 0.6 1.0 0.6 1 Specular 0 1 0 0 1 Emissive 0 1 0 0 1 Shininess 0 10" # light green
elif 'Skin' == part.name: materialColor = "Default Diffuse 1 1.0 0.7 0.6 1 Ambient 1 1.0 0.7 0.6 1 Specular 0 1 0 0 1 Emissive 0 1 0 0 1 Shininess 0 10 "
elif 'Fat' == part.name: materialColor = "Default Diffuse 1 1.0 1.0 0.5 1 Ambient 1 1.0 1.0 0.5 1 Specular 0 1 0 0 1 Emissive 0 1 0 0 1 Shininess 0 10 "
elif 'Muscle' == part.name: materialColor = "Default Diffuse 1 1.0 0.5 0.5 1 Ambient 1 1.0 0.5 0.5 1 Specular 0 1 0 0 1 Emissive 0 1 0 0 1 Shininess 0 10 "
else: materialColor = "Default Diffuse 1 1.0 0.6 0.6 1 Ambient 1 1.0 0.6 0.6 1 Specular 0 1 0 0 1 Emissive 0 1 0 0 1 Shininess 0 10" # bubblegum
myElement('OglModel', visu, name="Visual", material=materialColor)
myElement('IdentityMapping', visu)
def WriteAsGmsh(salomeDat, dirName='Sofa'):
""""Takes a SALOME dat file of a mesh and write a gmsh format v1.0 mesh"""
outMesh = os.path.join(dirName, salomeDat.replace('.dat', '.msh') )
# outMesh = 'Sofa/' + salomeDat.replace('.dat', '.msh')
gmsh=open(outMesh, 'w')
data=open(salomeDat, 'r')
nodes, elements = [int(s) for s in data.readline().split()]
gmsh.write('$NOD\n')
gmsh.write('{}\n'.format(nodes))
#
nodeConversion={}
for i in range(1, nodes+1):
line=data.readline().split()
nodeConversion[line[0]]=str(i)
line[0]=str(i)
line = ' '.join(line) + '\n'
gmsh.write(line)
#
gmsh.write('$ENDNOD\n')
gmsh.write('$ELM\n')
gmsh.write('{}\n'.format(elements))
for i in range(1,elements+1):
line = data.readline().split()
#
# line contains Salome elem #, Salome elem type, node list
# becomes Gmsh elem #, Gmsh elem type, PG#, EE#, #ofnodes, modified node list
#
numNodes=int(line[1][2:]) # line[1] == Salome elem type --> dimension of elem, 0, # of nodes
if line[1]=='203': line[1] = '2 0 1 3' # triangle, PG#, EE#, # of nodes
elif line[1]=='304': line[1] = '4 1 1 4' # tetrahedron, PG#, EE#, #nodes
else: print line[1], ' not supported by this script\n', 'See documentation.'; sys.exit()
line[-numNodes:] = [nodeConversion[i] for i in line[-numNodes:]] # change node list to renumbered nodes
line = ' '.join(line) + '\n'
gmsh.write(line)
gmsh.write('$ENDELM\n')
gmsh.close()
data.close()
os.remove(salomeDat)
def WriteAssembly(sofa, assembly, sofaDir='Sofa', tooBig = 15000):
for part in assembly.itervalues():
SofaObjectNode(sofa, part)
#
file = '{}.dat'.format(part.name)
if part.mat == 'rigid': part.mesh.ExportDAT(file)
else:
dim = part.mesh.MeshDimension()
if dim == 3: group, = part.mesh.GetGroupByName('{}_All_Volumes'.format(part.name)) # Makes list of one
elif dim == 2: group, = part.mesh.GetGroupByName('{}_All_Faces'.format(part.name)) # Makes list of one
part.mesh.ExportDAT(file, group)
if part.mesh.NbElements() > tooBig: print '\n{} has {} elements, may be too big for real time Sofa simulation!\n'.format(part.name, part.mesh.NbElements())
WriteAsGmsh(file, sofaDir)
#
# Order matters in a Sofa scene, cannot combine loops
#
for part in assembly.itervalues():
file = '{}.dat'.format(part.name)
for contact in part.Contact:
contactFile = '{}_@_{}'.format(contact, file)
print '{}_@_{}_{}Faces'.format(part.name, contact, 'Contact')
group, = part.mesh.GetGroupByName('{}_@_{}_{}Faces'.format(part.name, contact, 'Contact'))
assembly[contact].mesh.ExportDAT(contactFile, group)
WriteAsGmsh(contactFile, sofaDir)
if part.mat == 'rigid':
for tiedName in part.Ties:
tiedFile = '{}_@_{}'.format(tiedName, file)
group, = assembly[tiedName].mesh.GetGroupByName('{}_@_{}_{}Faces'.format(tiedName, part.name, 'Ties'))
assembly[tiedName].mesh.ExportDAT(tiedFile, group)
WriteAsGmsh(tiedFile, sofaDir)
#
elems = group.GetIDs()
nodes2 = set()
for elem in elems:
nodes = assembly[tiedName].mesh.GetElemNodes(elem)
for n in nodes:
nodes2.add(n)
#
nodes2 = sorted([i - 1 for i in nodes2]) # Sofa nodes start at zero, gmsh at 1
nodes1 = range(len(nodes2)) # We use a renumbered subset of the Ties mesh as our connection
#
myElement('AttachConstraint', sofa, object1="@{}/{}/Attach".format(part.name, tiedName), object2="@{}".format(tiedName),
twoWay="False", radius="1", indices1=str(nodes1)[1:-1], indices2=str(nodes2)[1:-1])
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 StlToMed.py args:ConnectivityFile\n'
sys.exit()
#
for part in Assembly.itervalues():
file = os.path.join(os.path.dirname(part.file),'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 MakeSofa(sofaFile='Sofa.scn', connectivityFile='Connectivity.xml', backgroundImage=None):
"""Read connectivity file, get MEDs and writes Sofa scene and the necessary mesh files"""
Assembly = ReadConnectivity(connectivityFile)
#
if os.path.dirname(sofaFile): # make sure it exists
if not os.path.isdir(os.path.dirname(sofaFile)):
os.mkdir(os.path.dirname(sofaFile))
else: # use Sofa dir
sofaFile = os.path.join('Sofa', sofaFile)
if not os.path.isdir('Sofa'): os.mkdir('Sofa')
#
sofa = et.Element('Node', {'name':'root','dt':'0.02','gravity':'0 0 0'})
tree = et.ElementTree(sofa)
#
import inspect
comment = et.Comment('This document was written by\n\t{}\n\ton {}'.format(inspect.getfile(inspect.currentframe()),
time.strftime('%b %d %Y at %H:%m:%S %p')))
sofa.append(comment)
#
BasicScene(sofa, backgroundImage) # standard options for any sofa scene
WriteAssembly(sofa, Assembly, os.path.dirname(sofaFile)) # parts from connectivity
#
try: tree.write(sofaFile, xml_declaration=True, pretty_print=True)
except TypeError: tree.write(sofaFile, xml_declaration=True)
print '\nWrote', sofaFile
#
CloseSalome()
if __name__ == '__main__':
MakeSofa(*sys.argv[1:])