<html><head><meta name="color-scheme" content="light dark"></head><body><pre style="word-wrap: break-word; white-space: pre-wrap;">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:args:ConnectivityFile,SofaFile(optional),backgroundImage(optional)\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('RequiredPlugin', parent, pluginName='SofaMiscCollision')
	myElement('RequiredPlugin', parent, pluginName='SofaOpenglVisual')
	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="0.1", contactDistance="0.5")	# This one might actually need to be changed.
	myElement('EulerImplicitSolver', 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, groupNumber=1):
	"""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")
		fixedComment = et.Comment('&lt;PartialFixedConstraint template="Rigid3d" name="partialFixedConstraint0"  indices="0"  fixedDirections="0 0 0 0 0 0" /&gt;')
		object.append(fixedComment)
		movementComment = et.Comment('&lt;LinearMovementConstraint template="Rigid3d" name="Displacement" keyTimes="t1 t2 t3"\n movements="x1 y1 z1 rx1 ry1 rz1\n x2 y2 z2 rx2 ry2 rz2\n x3 y3 z3 rx3 ry3 rz3"/&gt;')
		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('OglModel', visu, template="Vec3d", 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='{}_@_{}_Ties.msh'.format(tie,name), name="loader")
			myElement('MechanicalObject', t, src="@loader", name="Attach")
			myElement('RigidMapping', t, input="@..",output="@Attach")
		for contact in part.Contact:
			c = myElement('Node', object, name='{}Contact'.format(contact) )
			myElement('MeshGmshLoader', c, filename='{}_@_{}_Contact.msh'.format(name,contact), name="loader")
			myElement('MeshTopology', c, src="@loader")
			myElement('MechanicalObject', c, src="@loader", name="Contact")
			myElement('RigidMapping', c, input="@..", output="@Contact")
			myElement('TriangleCollisionModel', c, name="Triangles For Collision", contactStiffness="100000", group=str(groupNumber))
			myElement('LineCollisionModel', c, name="Lines For Collision", contactStiffness="100000", group=str(groupNumber))
			myElement('PointCollisionModel', c, name="Points For Collision", contactStiffness="100000", group=str(groupNumber))
	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('TriangleCollisionSet', surf, contactStiffness='3000')
		# myElement('LineCollisionModel',   surf, contactStiffness='3000')
		# myElement('PointCollisionModel',  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 "	# tanish
		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 "	# yellowish
		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 "	# bubblegum
		elif 'Bone'   ==	part.name:	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
		else:							materialColor = "Default Diffuse 1 1.0 0.0 1.0 1 Ambient 1 1.0 0.0 1.0 1 Specular 0 1 0 0 1 Emissive 0 1 0 0 1 Shininess 0 10"	# orange
		myElement('OglModel', visu, name="Visual", material=materialColor)
		myElement('IdentityMapping', visu)
		for i, contact in enumerate(part.Contact):
			c = myElement('Node', object, name='{}Contact'.format(contact))
			myElement('MeshGmshLoader', c, filename='{}_@_{}_Contact.msh'.format(name, contact), name="loader")
			myElement('MeshTopology', c, src="@loader")
			myElement('MechanicalObject', c, src="@loader", name="Contact")
			myElement('BarycentricMapping', c, input="@..", output="@Contact")
			myElement('TriangleCollisionModel', c, name="Triangles For Collision", contactStiffness="3000", group=str(groupNumber))
			myElement('LineCollisionModel', c, name="Lines For Collision", contactStiffness="3000", group=str(groupNumber))
			myElement('PointCollisionModel', c, name="Points For Collision", contactStiffness="3000", group=str(groupNumber))


def WriteAsGmsh(salomeDat, dirName='Sofa'):
	""""Takes a SALOME dat file of a mesh and write a gmsh format v1.0 mesh"""
	path, salomeFile = os.path.split(salomeDat)
	outMesh = os.path.join(path, dirName, salomeFile.replace('.dat', '.msh') )
	with open(outMesh, 'w') as gmsh, open(salomeDat, 'r') as data:
		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 --&gt; dimension of elem, 2 digit # of nodes
			if   line[1]=='203':	line[1] = '2 0 1 3'  # triangle,    PG#, EE#, # of nodes
			elif line[1]=='206':							 # quadratic triangle
				numNodes = 3
				line = line[0:numNodes+2]
				line[1] = '2 0 1 3'
			elif line[1]=='304':	line[1] = '4 1 1 4'  # tetrahedron, PG#, EE#, #nodes
			elif line[1]=='310':						 # quadratic tetrahedron
				numNodes = 4
				line = line[0:numNodes+2]
				line[1] = '4 1 1 4'
			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')

def WriteAssembly(sofa, assembly, sofaDir='Sofa', tooBig = 15000):
	for num, part in enumerate(assembly.values()):
		SofaObjectNode(sofa, part, num)
		#
		file = part.med.replace('.med','.dat')
		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() &gt; 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 name, part in assembly.items():
		print(name, 'Contact')
		for contact in part.Contact:
			print('\t', contact)
			contactFile = part.med.replace('.med', '{}_@_{}_{}.dat'.format(name , contact,'Contact') )
			print('		{}_@_{}_{}Faces'.format(name, contact, 'Contact'))
			group, = part.mesh.GetGroupByName('{}_@_{}_{}Faces'.format(name, contact, 'Contact'))
			assembly[contact].mesh.ExportDAT(contactFile, group)
			WriteAsGmsh(contactFile, sofaDir)
			os.remove(contactFile)
		if part.mat == 'rigid':
			print(name, 'Ties')
			for tiedName in part.Ties:
				print('\t', tiedName)
				switch=False
				if part.mat != 'rigid':
					name,tiedName = tiedName,name	# switch names
					switch=True
					print('switching!')
				if switch:  tiedFile = part.med.replace('.med', '{}_@_{}_{}.dat'.format(tiedName, name, 'Ties') )
				else:       tiedFile = part.med.replace('.med',     '{}_@_{}_{}.dat'.format(tiedName, name, 'Ties') )
				print('	{}_@_{}_{}Faces'.format(tiedName, name, 'Ties'))
				group, = assembly[tiedName].mesh.GetGroupByName('{}_@_{}_{}Faces'.format(tiedName, name, 'Ties'))
				assembly[tiedName].mesh.ExportDAT(tiedFile, group)
				WriteAsGmsh(tiedFile, sofaDir)
				os.remove(tiedFile)
				#
				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(name, tiedName), object2="@{}".format(tiedName),
							twoWay="False", radius="1",
							# indices1=str(nodes1)[1:-1],  " ".join(str(n) for n in nodes1)
							# indices2=str(nodes2)[1:-1],
							# constraintFactor=str(range(len(nodes1))) )
							indices1=" ".join(str(n) for n in nodes1),
							indices2=" ".join(str(n) for n in nodes2),
							constraintFactor= " ".join(str(1) for n in range(len(nodes1)))  )
				if switch:	name,tiedName = tiedName,name	# switch back

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 MedToSofa.py args:ConnectivityFile,SofaFile(optional),backgroundImage(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 MakeSofa(connectivityFile='Connectivity.xml', sofaFile='Sofa.scn', backgroundImage=None):
	"""Read connectivity file, get MEDs and writes Sofa scene and the necessary mesh files"""
	Assembly = ReadConnectivity(connectivityFile)
	#
	matList = [];	rigidList = []
	for name, part in Assembly.items():
		if part.mat.lower() == 'rigid':	rigidList.append(name)
		else:							matList.append(name)
	#
	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:
		f = open(sofaFile, 'w')
		f.write('&lt;?xml version="1.0" encoding="ISO-8859-1"?&gt;\n')
		f.write(et.tostring(sofa).decode().replace('&gt;', '&gt;\n'))
		f.close()
	print('\nWrote', sofaFile)
	#
	CloseSalome()

if __name__ == '__main__':
	MakeSofa(*sys.argv[1:])

</pre></body></html>