<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:Abaqus.inp,ConnectivityFile  \n\n\n'.format( inspect.getfile(inspect.currentframe()) )
	sys.exit()
from MEDLoader import MEDLoader

from collections import Counter
import itertools
chain = itertools.chain.from_iterable

import StlToMed

salome.salome_init()
theStudy = salome.myStudy

notebook = salome_notebook.NoteBook(theStudy)
smesh = smeshBuilder.New(theStudy)

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 WritePart(part, abaqus, formatting=False):
	"""Output an Abaqus Part"""
	abaqus.write('*Part, Name = {}\n'.format(part.name))
	mesh = part.mesh
	dimension = mesh.MeshDimension()
	# Nodes
	abaqus.write('	*Node\n')
	if formatting:	writeString ='	{:6d}, {: 12.8e}, {: 12.8e}, {: 12.8e}\n'
	else:			writeString ='	{}, {}, {}, {}\n'
	nodes = MEDLoader.ReadUMeshFromFile(part.med, 0).getCoords()
	for node in xrange(len(nodes)):
		abaqus.write(writeString.format(node+1, *nodes.getTuple(node)) )
	#
	if part.mat == 'rigid':
		abaqus.write(writeString.format(node+2, *StlToMed.FindCenter(mesh)))
		abaqus.write('	*NSet, NSet=RigidNode\n\t{}\n'.format(node+2))
	# Node Sets
	groups = mesh.GetGroups()
	for group in groups:
		name = group.GetName()
		if name.endswith('Nodes'):
			abaqus.write('	*NSet, NSet={}\n'.format(name))
			WriteListOfItems(sorted(group.GetListOfID()), abaqus,formatting)
	# Elements and Surfaces
	medMesh = MEDLoader.ReadUMeshFromFile(part.med, 0)
	if medMesh.getMeshDimension() == 2:
		nodeCount = 3
		# Elements
		abaqus.write('	*Element, Type=S3\n')
		meshConnectivity = medMesh.getNodalConnectivity()
		numElems = medMesh.getNumberOfCells()
		if formatting:		writeString = '	{:6d}, {:6d}, {:6d}, {:6d}\n'
		else:				writeString = '	{}, {}, {}, {}\n'
		for elem in xrange(numElems):
			index = elem * (nodeCount+1)
			connectivity = [meshConnectivity[index + c] + 1 for c in range(1, nodeCount + 1)]	# sliceing does't work quite right as this is a Salome structure
			abaqus.write(writeString.format(elem+1, *connectivity))
		# Elements Sets and Surfaces(built from Element Sets)
		for group in groups:
			name = group.GetName()
			if name.endswith('Surface') or name.endswith('Faces'):
				abaqus.write('	*ElSet, ElSet=_{}, internal\n'.format(name))
				WriteListOfItems(sorted(group.GetListOfID()), abaqus, formatting)
				#
				abaqus.write('	*Surface, type=Element, name={}\n'.format(name))
				abaqus.write('	_{}\n'.format(name))
			elif name.endswith('All'):				pass
			elif name.endswith('Nodes'):			pass
			else:
				print 'Unknown Surface type', name
				sys.exit()
		abaqus.write('	*Shell Section, elset=_{}Surface, material={}\n'.format(part.name, part.name))
	elif medMesh.getMeshDimension() == 3:
		nodeCount = 4
		# Elements
		abaqus.write('	*Element, TYPE=C3D4\n')
		meshConnectivity = medMesh.getNodalConnectivity()
		numElems = medMesh.getNumberOfCells()
		if formatting:	writeString = '	{:6d}, {:6d}, {:6d}, {:6d}, {:6d}\n'
		else:			writeString = '	{}, {}, {}, {}, {}\n'
		# variable length is slightly slower than fixed  #	writeString = ('	{:6d}'+', {:6d}'*nodeCount+'\n')
		for elem in xrange(numElems):
			index = elem * (nodeCount + 1)
			connectivity = [meshConnectivity[index + c] + 1 for c in range(1, nodeCount + 1)]
			connectivity[1], connectivity[2] = connectivity[2], connectivity[1]  # swap nodes n2 and n3
			abaqus.write(writeString.format(elem+1, *connectivity))
		#
		# Elements Sets and Surfaces(built from Element Sets)
		for group in groups:
			name = group.GetName()
			numFaces = mesh.NbFaces()
			if   name.endswith('Surface') or name.endswith('Faces'):
				faceSets = {i: [] for i in range(nodeCount)}
				for i, elem in enumerate(sorted(group.GetListOfID())):
					triNodes = set(mesh.GetElemNodes(elem))  # a list of 3 nodes in the face, a set so order doesn't matter
					#
					tetraFromNodes=[]
					for node in triNodes:
						tetraFromNodes.append( mesh.GetNodeInverseElements(node) )	# make a list of element in lists
					joinedTetraList = chain(tetraFromNodes)							# flatten the lists
					count = Counter(joinedTetraList)		# count the number of times the elements show up in the list
					count[elem] = 0	# remove the count of the face element that always shows up.  I don't want to find it accidentally.
					tetra = count.most_common(1)[0][0]		# this will be the tetrahedral element
					#
					for j in range(nodeCount):
						faceNodes = set(mesh.GetElemFaceNodes(tetra, j))
						if faceNodes == triNodes:				# match node is tetra face and triangle group face
							faceSets[j].append(tetra-numFaces)	# add tetra to faceSet associated with the face that is in the surface
							#										-numFaces because Abaqus doesn't count faces as seperate elements and Salome does
							break			# only one face of the tetrahedrin can be the matching face
				#
				for i, fSet in faceSets.iteritems():  # Build internal (hidden) element sets used to build surface
					abaqus.write('	*ElSet, ElSet={}_{}, internal\n'.format(name, i + 1))
					WriteListOfItems(sorted(fSet), abaqus, formatting)
				#
				abaqus.write('	*Surface, type=Element, name={}\n'.format(name))
				for i in faceSets:		abaqus.write('	{}_{}, S{}\n'.format(name, i + 1, i + 1))
			elif name.endswith('Volumes'):
				abaqus.write('	*ElSet, ElSet={}\n'.format(name))
				WriteListOfItems(sorted(group.GetListOfID()), abaqus, formatting)
			elif name.endswith('Nodes'):	pass
			else:
				print 'Unknown Surface type', name
				sys.exit()
		abaqus.write('	*Solid Section, elset={}_ALL_VOLUMES, material={}\n'.format(part.name, part.name))
	else:
		print 'Bad mesh Dimension:', part.mesh.MeshDimension()
		sys.exit()
	#
	abaqus.write('*End Part\n'.format(part.name))

def WritePartTimimg(part, abaqus, formatting=False):
	"""Output an Abaqus Part"""
	print ''
	begin = time.clock()
	abaqus.write('*Part, Name = {}\n'.format(part.name))
	dimension = part.mesh.MeshDimension()
	# Nodes
	start = time.clock()
	abaqus.write('	*Node\n')
	if formatting:	writeString ='	{:6d}, {: 12.8e}, {: 12.8e}, {: 12.8e}\n'
	else:			writeString ='	{}, {}, {}, {}\n'
	nodes = MEDLoader.ReadUMeshFromFile(part.med, 0).getCoords()
	for node in xrange(len(nodes)):
		abaqus.write(writeString.format(node+1, *nodes.getTuple(node)) )
	#
	if part.mat == 'rigid':
		abaqus.write(writeString.format(node+2, *StlToMed.FindCenter(part.mesh)))
		abaqus.write('	*NSet, NSet=RigidNode\n\t{}\n'.format(node+2))
	end = time.clock()
	print 'Nodes took', end-start
	# Node Sets
	groups = part.mesh.GetGroups()
	for group in groups:
		name = group.GetName()
		if name.endswith('Nodes'):
			abaqus.write('	*NSet, NSet={}\n'.format(name))
			WriteListOfItems(sorted(group.GetListOfID()), abaqus,formatting)
	# Elements and Surfaces
	start = time.clock()
	faceSetTime = 0
	medMesh = MEDLoader.ReadUMeshFromFile(part.med, 0)
	mesh = part.mesh
	if medMesh.getMeshDimension() == 2:
		nodeCount = 3
		# Elements
		abaqus.write('	*Element, Type=S3\n')
		meshConnectivity = medMesh.getNodalConnectivity()
		numElems = medMesh.getNumberOfCells()
		if formatting:		writeString = '	{:6d}, {:6d}, {:6d}, {:6d}\n'
		else:				writeString = '	{}, {}, {}, {}\n'
		for elem in xrange(numElems):
			index = elem * (nodeCount+1)
			connectivity = [meshConnectivity[index + c] + 1 for c in range(1, nodeCount + 1)]	# sliceing does't work quite right as this is a Salome structure
			abaqus.write(writeString.format(elem+1, *connectivity))
		# Elements Sets and Surfaces(built from Element Sets)
		for group in groups:
			name = group.GetName()
			if name.endswith('Surface') or name.endswith('Faces'):
				abaqus.write('	*ElSet, ElSet=_{}, internal\n'.format(name))
				WriteListOfItems(sorted(group.GetListOfID()), abaqus, formatting)
				#
				abaqus.write('	*Surface, type=Element, name={}\n'.format(name))
				abaqus.write('	_{}\n'.format(name))
			elif name.endswith('All'):				pass
			elif name.endswith('Nodes'):			pass
			else:
				print 'Unknown Surface type', name
				sys.exit()
		abaqus.write('	*Shell Section, elset=_{}Surface, material={}\n'.format(part.name, part.name))
	elif medMesh.getMeshDimension() == 3:
		nodeCount = 4
		# Elements
		go = time.clock()
		abaqus.write('	*Element, TYPE=C3D4\n')
		meshConnectivity = medMesh.getNodalConnectivity()
		numElems = medMesh.getNumberOfCells()
		if formatting:	writeString = '	{:6d}, {:6d}, {:6d}, {:6d}, {:6d}\n'
		else:			writeString = '	{}, {}, {}, {}, {}\n'
		# variable length is slightly slower than fixed  #	writeString = ('	{:6d}'+', {:6d}'*nodeCount+'\n')
		for elem in xrange(numElems):
			index = elem * (nodeCount + 1)
			connectivity = [meshConnectivity[index + c] + 1 for c in range(1, nodeCount + 1)]
			connectivity[1], connectivity[2] = connectivity[2], connectivity[1]  # swap nodes n2 and n3
			abaqus.write(writeString.format(elem+1, *connectivity))
		stop = time.clock()
		print 'Elements took', stop - go
		#
		# Elements Sets and Surfaces(built from Element Sets)
		go = time.clock()
		triNodesTime = [];		tetraFromNodesTime = [];			facesTime = []
		for group in groups:
			name = group.GetName()
			numFaces = mesh.NbFaces()
			if   name.endswith('Surface') or name.endswith('Faces'):
				a = time.clock()
				faceSets = {i: [] for i in range(nodeCount)}
				for i, elem in enumerate(sorted(group.GetListOfID())):
					here = time.clock()
					triNodes = set(mesh.GetElemNodes(elem))  # a list of 3 nodes in the face, a set so order doesn't metter
					there = time.clock()
					triNodesTime.append(there-here)


					tetraFromNodesList=[]
					here = time.clock()
					for node in triNodes:
						tetraFromNodesList.append( mesh.GetNodeInverseElements(node) )	# make a list of element in lists
					joinedTetraList = chain(tetraFromNodesList)							# flatten the lists
					count = Counter(joinedTetraList)		# count the number of times the elements show up in the list
					count[elem] = 0	# remove the count of the face element that always shows up.  I don't want to find it accidentally.
					tetra = count.most_common(1)[0][0]		# this will be the tetrahedral element

					there = time.clock()
					tetraFromNodesTime.append(there-here)

					here = time.clock()
					for j in range(nodeCount):
						faceNodes = set(mesh.GetElemFaceNodes(tetra, j))
						if faceNodes == triNodes:				# match node is tetra face and triangle group face
							faceSets[j].append(tetra-numFaces)	# add tetra to faceSet associated with the face that is in the surface
							#										-numFaces because Abaqus doesn't count faces as seperate elements and Salome does
							break			# only one face of the tetrahedrin can be the matching face
					there = time.clock()
					facesTime.append(there-here)

				b = time.clock()
				print 'Making FaceSet took {}	{}'.format(b-a, name)
				faceSetTime += b-a
				#
				for i, fSet in faceSets.iteritems():  # Build internal (hidden) element sets used to build surface
					abaqus.write('	*ElSet, ElSet={}_{}, internal\n'.format(name, i + 1))
					WriteListOfItems(sorted(fSet), abaqus, formatting)
				#
				abaqus.write('	*Surface, type=Element, name={}\n'.format(name))
				for i in faceSets:		abaqus.write('	{}_{}, S{}\n'.format(name, i + 1, i + 1))
			elif name.endswith('Volumes'):
				abaqus.write('	*ElSet, ElSet={}\n'.format(name))
				WriteListOfItems(sorted(group.GetListOfID()), abaqus, formatting)
			elif name.endswith('Nodes'):	pass
			else:
				print 'Unknown Surface type', name
				sys.exit()
		stop = time.clock()
		print 'ElementSet and Surfaces took', stop - go
		print 'Face Set Time ', faceSetTime

		print  'triNodes took', sum(triNodesTime)
		print  'tetraFromNodes took', sum(tetraFromNodesTime)
		print  'faces took', sum(facesTime)


		abaqus.write('	*Solid Section, elset={}_ALL_VOLUMES, material={}\n'.format(part.name, part.name))
	else:
		print 'Bad mesh Dimension:', part.mesh.MeshDimension()
		sys.exit()
	end = time.clock()
	print 'Element and Sets and Surfaces took', end-start
	#
	abaqus.write('*End Part\n'.format(part.name))
	finish = time.clock()
	print 'Total Writing Part took', finish-begin, '\n'

def WriteListOfItems(listOfItems, abaqus, formatting=False):
	"""Writes the set of nodes or elements...reusing code is good
		one line is a good balance between efficiency and memory use"""
	if formatting:	writeString = '	{:6d},	{:6d},	{:6d},	{:6d},	{:6d},	{:6d},	{:6d},	{:6d},	{:6d},	{:6d},\n'
	else:			writeString = '	{},	{},	{},	{},	{},	{},	{},	{},	{},	{},\n'	# no formatting is much faster, but harder to read output
	#
	n = 10		# number of elements on a line.  I prefer 10 to the 16 allowable in Abaqus.
	for i in xrange( len(listOfItems)/n ):
		abaqus.write(writeString.format(*listOfItems[i*n:(i+1)*n]))
	# And the last partial line
	i=len(listOfItems)/n
	if formatting:	abaqus.write(('	{:6d},'*len(listOfItems[i*n:(i+1)*n]) + '\n').format(*listOfItems[i*n:(i+1)*n]))
	else:			abaqus.write(('	{},'*len(listOfItems[i*n:(i+1)*n]) + '\n').format(*listOfItems[i*n:(i+1)*n]))

def WriteAssembly(assembly, abaqus, name='AutoAssembly', formatting=False):
	"""assembly is from xml reading
		abaqus is a file already opened for writing"""
	for part in assembly.itervalues():
		print 'Writing Part {}'.format(part.name)
		WritePart(part, abaqus, formatting)
	print 'Creating Assembly'
	abaqus.write('*Assembly, Name = {}\n'.format(name))
	for part in assembly.itervalues():
		abaqus.write('	*Instance, Name={}, Part={}\n'.format(part.name, part.name))
		# possible part positioning data could go here, but that not relevant for the current project
		abaqus.write('	*End Instance\n')
	# Make Rigid bodies
	for name, part in assembly.iteritems():
		if part.mat == 'rigid':
			abaqus.write('	*Rigid Body, ref node={}.RigidNode, elset={}._{}Surface\n'.format(name, name, name))
	# Set ties
	for master, part in assembly.iteritems():
		switched = False
		for slave in part.Ties:
			minSize, maxSize = assembly[master].mesh.GetMinMax(SMESH.FT_Length2D)  # approximate 'characteristic length' of the mesh
			try:				multiplier = part.Ties[slave]['proximity']
			except KeyError:	multiplier=1.0
			positionTol = (maxSize + minSize) * multiplier
			#
			if assembly[slave].mat == 'rigid':	# rigid body need to be the master surface
				master, slave = slave, master
				switched = True
			name = '{}_Tied_{}'.format(master, slave)
			abaqus.write('	*Tie, Name={}, Adjust=no, Position Tolerance={}\n'.format(name, positionTol))
			abaqus.write('	{}.{}_@_{}_TiesFaces, {}.{}_@_{}_TiesFaces\n'.format(slave,slave,master,	master,master,slave))
			if switched:						# return master to original value if changed
				master = slave					#	no need to change slave as loop is ending
	#
	abaqus.write('*End Assembly\n')

def WriteContact(assembly,abaqus, name='AutoAssembly'):
	print 'Creating Contact conditions'
	#
	for master, part in assembly.iteritems():
		switched = False
		for slave in part.Contact:
			if assembly[slave].mat == 'rigid':
				master, slave = slave, master
			contact = '{}_Contact_{}'.format(master, slave)
			abaqus.write('*Surface Interaction, Name={}\n'.format(contact))
			abaqus.write('*Surface Behavior, Augmented Lagrange\n')
			abaqus.write('*Contact\n')
			abaqus.write('*Contact Inclusions\n')
			abaqus.write('	{}.{}.{}_@_{}_ContactFaces,' .format(name,slave,slave,master))
			abaqus.write('	{}.{}.{}_@_{}_ContactFaces\n'.format(name,master,master,slave))
			if switched: 				# return master to original value if changed
				master = slave			#	no need to change slave as loop is ending

def WriteMaterials(assembly,abaqus):
	"""Only writes a generic elastic material for each part
		The user should alter the material definitions as needed."""
	print 'Creating basic Material Definitions'
	for part in assembly.itervalues():
		abaqus.write('*Material, Name={}\n'.format(part.name))
		abaqus.write('	*Elastic\n')
		abaqus.write('	100, 0.3\n')

def WriteBoundaries(assembly,abaqus):
	"""Creates a Step and sets boundary conditions for rigid objects"""
	print 'Creating basic rigid body Boundary conditions.'
	abaqus.write('*Step, name = Step - 1, nlgeom = YES\n')
	abaqus.write('	*Static\n')
	abaqus.write('	0.1, 1., 1e-05, 0.1\n')		# Initial time increment, Total Time period, Min increment, Max increment
	DoF = [i+1 for i in range(6)]
	for part in assembly.itervalues():
		if part.mat == 'rigid':
			abaqus.write('	*Boundary\n'.format(part.name))
			for d in DoF:
				abaqus.write('		{}.RigidNode, {}, {}\n'.format(part.name, d, d))
	abaqus.write('*End Step\n')

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 MakeAbaqus(abaqusFile='Abaqus.inp', connectivityFile='Connectivity.xml'):
	"""Read connectivity file, get MEDs and write Abaqus input file"""
	Assembly = ReadConnectivity(connectivityFile)
	#
	if not os.path.dirname(abaqusFile):
		abaqusFile = os.path.join('Abaqus', abaqusFile)
		if not os.path.exists('Abaqus'):    os.makedirs('Abaqus')
	abaqus = open(abaqusFile, 'w')
	#
	WriteAssembly(Assembly, abaqus, 'OpenKnee')
	WriteContact(Assembly, abaqus, 'OpenKnee')
	WriteMaterials(Assembly, abaqus)
	WriteBoundaries(Assembly, abaqus)
	#
	abaqus.close()
	print 'Writing to {} file'.format(abaqusFile)
	#
	CloseSalome()

def runTest(abaqusFile='MEDCoupling.inp', connectivityFile='Connectivity.xml'):
	print '\n\n\n  WARNING - Testing MEDCouplng speed up\n\n\n'
	import ConnectivityXml as xml
	print 'Reading from connectivity file {}'.format(connectivityFile)
	Assembly = xml.ReadConnectivity(connectivityFile)
	xml.RemoveDoubleCounting(Assembly)
	#
	if not os.path.dirname(abaqusFile):
		abaqusFile = os.path.join('Abaqus', abaqusFile)
		if not os.path.exists('Abaqus'):    os.makedirs('Abaqus')
	abaqus = open(abaqusFile, 'w')
	#
	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
	#
	for part in Assembly.itervalues():
		print '\nWriting Part {}'.format(part.name)
		WritePartTimimg(part, abaqus)
	#
	abaqus.close()
	#
	print 'Writing to {} file'.format(abaqusFile)
	print '\n\n\n  WARNING - Testing MEDCouplng speed up\n\n\n'
	# try to close out of salome
	try:
		from killSalomeWithPort import killMyPort
		killMyPort(2810)
	except:	pass

if __name__ == '__main__':
	# runTest(*sys.argv[1:])
	MakeAbaqus(*sys.argv[1:])

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