# -*- coding: utf-8 -*-
import sys
import os
import inspect
import time
try:
	import salome
	import salome_notebook
	import SMESH
	from salome.smesh import smeshBuilder
except ImportError:
	import inspect
	print('\n\n\nUsage: salome {} args:ConnectivityFile\n\n\n'.format( inspect.getfile(inspect.currentframe()) )  )
	sys.exit()

try:					from MEDLoader import MEDLoader
except ImportError:		import MEDLoader
import numpy as np

salome.salome_init()
smesh = smeshBuilder.New()

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 CutGroups(mesh, main, tool):
	"""takes a mesh and two (lists of) groups, main and tool
		removes elements from main that are also in tool"""
	name = main[0].GetName()
	new = mesh.CutListOfGroups(main, tool, name)
	if len(new.GetListOfID()) > 0:
		mesh.RemoveGroup(main[0])
		# mesh.RemoveGroup(oldNodes)
	else:
		print(name)
		print(    'WARNING!		No elements remain in contact surface!\n')
		log.write(name)
		log.write('WARNING!		No elements remain in contact surface!\n')
		mesh.RemoveGroup(new)
		# mesh.RemoveGroup(newNodes)

def FindCenter(mesh, group=None):
	"""Takes a mesh and list of nodes on that mesh
		get node positions and calculates the mean
		returns the geometric center, i.e. barycenter"""
	if group:	nodes = group.GetNodeIDs()	# only nodes within group
	else:		nodes = mesh.GetNodesId()	# all nodes in mesh
	positions = np.ones( (len(nodes),3) )
	for i, node in enumerate(nodes):
		positions[i] = mesh.GetNodeXYZ(node)
	center = np.mean(positions, axis=0)
	return center

def FindOrientedNormals(mesh, nodes, vector, threshold=0):
	"""nodes is a list of nodes that will be used to find the elements and then normals
		vector is the vector to compare element normal against
		returns a list of faces ids use to make group"""
	#  ensure vector is normalized
	v = vector
	mag = (v[0]**2+v[1]**2+v[2]**2)**0.5
	vector = np.array([v[0]/mag, v[1]/mag, v[2]/mag])
	# uses nodes to make a node group
	tempNodes = mesh.MakeGroupByIds('TempNodes', SMESH.NODE, nodes)
	# uses nodes group to make a faces group
	tempFaces = mesh.CreateDimGroup(tempNodes, SMESH.FACE, "Temp", SMESH.MAJORITY, 0)
	# extact faces from face group
	faces = tempFaces.GetListOfID()
	normals = np.ones((len(faces), 3))
	for i, face in enumerate(faces):
		normals[i] = mesh.GetFaceNormal(face, normalized=True)
	dot = np.dot(normals, vector)
	#
	for i in range(len(dot)-1, -1, -1):  # iterate backwards so the higher indices don't change mid loop!
		if dot[i] <= threshold:
			faces.pop(i)
	# remove temporary groups from mesh
	mesh.RemoveGroup(tempNodes)
	mesh.RemoveGroup(tempFaces)
	#
	return faces	# a list of faces ids use to make group

def FindVectorFaces(part):
	"""Similar to FindOrientedNormals, but expanded for use with MakeVectorSurface
		vector is the vector to compare element normal against
		nodes is a list of nodes numbers that will be used to find the elements and then normals
		nodeArray is a list of node coordinates that will be used to determine if the are above or below the barycenter
		returns a list of faces ids use to make group"""
	mesh = part.mesh
	#  ensure vector is normalized
	v = part.Other['Vector']
	mag = (v[0]**2+v[1]**2+v[2]**2)**0.5
	vector = np.array([v[0]/mag, v[1]/mag, v[2]/mag])
	#
	#
	try:	AllNodeArray = MEDLoader.ReadUMeshFromFile(part.med, 0).getCoords().toNumPyArray()
	except AttributeError:
		coords = MEDLoader.ReadUMeshFromFile(part.med, 0).getCoords()
		AllNodeArray = np.zeros((len(coords), 3))
		for i, node in enumerate(coords):
			AllNodeArray[i] = node
	#
	allNodes  = mesh.GetElementsByType(SMESH.NODE)	# all nodes
	#
	filter = smesh.GetFilter(SMESH.FACE, SMESH.FT_FreeFaces)
	surface = mesh.GroupOnFilter(SMESH.NODE, "Surface", filter)
	surfNodes = surface.GetListOfID()
	mesh.RemoveGroup(surface)
	#
	center = FindCenter(mesh)
	diff = AllNodeArray-center
	diffdot = np.dot(diff,vector)	# positive values on the normal side of CoM
	#
	for i in range(len(diffdot)-1, -1, -1):  	# this iterates backwards so the higher indices don't change mid loop!
		if diffdot[i] <= 0:
			allNodes.pop(i)
	#
	nodes = set(surfNodes)
	nodes = nodes.intersection(allNodes)
	#
	# use nodes to make a node group
	tempNodes = mesh.MakeGroupByIds('TempNodes', SMESH.NODE, list(nodes))
	# use nodes group to make a faces group
	tempFaces = mesh.CreateDimGroup(tempNodes, SMESH.FACE, "Temp", SMESH.AT_LEAST_ONE, 0)
	# extact faces from face group
	faces = tempFaces.GetListOfID()
	normals = np.ones((len(faces), 3))
	for i, face in enumerate(faces):
		normals[i] = mesh.GetFaceNormal(face, normalized=True)
	#
	angle = np.arccos( np.dot(normals, vector) )
	for i in range(len(angle)-1, -1, -1):  # this iterates backwards so the higher indices don't change mid loop!
		if angle[i]>0.35:	# about 20 degrees
			faces.pop(i)
	# remove temporary groups from mesh
	mesh.RemoveGroup(tempNodes)
	mesh.RemoveGroup(tempFaces)
	#
	return faces	# a list of faces ids use to make group

def FindProximity(part1, part2, multiplier):
	"""mesh1 and mesh2 are mesh data,
		multiplier is a mesh characteristic length modifier
		returns list of nodes in proximity"""
	nodes1, array1 = NodesAsArray(part1)
	nodes2, array2 = NodesAsArray(part2)	# for point based surface finding nodes2 should be 1, and array2 should be part2 ie the point
	#
	if multiplier:
		minSize, maxSize= part1.mesh.GetMinMax(SMESH.FT_Length2D)	# approximates a 'characteristic length' of the mesh
		meshSize = (maxSize+minSize)/2.0 * multiplier
	else:
		box = part1.mesh.BoundingBox()
		meshSize = ((box[3] - box[0]) ** 2 + (box[4] - box[1]) ** 2 + (box[5] - box[2]) ** 2) ** 0.5
	#
	proximityNodes = []
	for n in nodes1:
		node = array1[n-1]
		minDiff = np.min(np.linalg.norm(array2-node, axis=1))
		if minDiff < meshSize:
			proximityNodes.append(n)
	if not proximityNodes:		# the given multiplier is too small, help find a better one!
		box = part1.mesh.BoundingBox()
		print(    '	No nodes in proxmitiy for {} to {}'.format(part1.name, part2.name),)
		log.write('	No nodes in proxmitiy for {} to {}'.format(part1.name, part2.name))
		boxSize = ((box[3] - box[0]) ** 2 + (box[4] - box[1]) ** 2 + (box[5] - box[2]) ** 2) ** 0.5
		for newSize in np.logspace( np.log10(meshSize), np.log10(boxSize) , 10):	# increase multiplier logarithmically
			# print(type(newSize))
			for n in nodes1:
				node = array1[n - 1]
				minDiff = np.min(np.linalg.norm(array2 - node, axis=1))
				if minDiff < newSize:
					proximityNodes.append(n)
			if proximityNodes:
				print(    ', estimated new multiplier: ', newSize/(meshSize/multiplier), 'old multiplier', multiplier)
				log.write(', estimated new multiplier: {}\n'.format(newSize/(meshSize/multiplier)))
				break	# stop when some nodes are found
	if not proximityNodes:		# objects are really far apart and proximity is NOT a good criteria
		print(    ', within total size of {}!'.format(part1.name))
		log.write(', within total size of {}!\n'.format(part1.name))
	return proximityNodes

def MakeGroups(Assembly):
	print('Making groups:');		log.write('Making groups:\n')
	for name, part in Assembly.items():
		print(name);		log.write('{}\n'.format(name))
		for tiedName, tie in part.Ties.items():
			if 'local_normal' in tie:		MakeLocalGroup    (part, Assembly[tiedName])
			elif 'local'      in tie:		MakeLocalGroup   (part, Assembly[tiedName])
			elif 'new'        in tie:		MakeNewGroup     (part, Assembly[tiedName])
			elif 'contains'   in tie:		MakeContainsGroup(part, Assembly[tiedName])
			elif 'normals'    in tie:		MakeNormalsGroup (part, Assembly[tiedName])
			else:							MakeProximityGroup(part, Assembly[tiedName])
		tiedGroups = [g for g in part.mesh.GetGroups() if 'TiesFaces' in g.GetName()]
		tiedNodeGroups = [g for g in part.mesh.GetGroups() if 'TiesNodes' in g.GetName()]
		for contactName, contact in part.Contact.items():
			if 'local_normals' in contact:		MakeLocalGroup    (part, Assembly[contactName], 'Contact')
			elif 'local'       in contact:		MakeLocalGroup    (part, Assembly[contactName], 'Contact')
			elif 'new'         in contact:		MakeNewGroup      (part, Assembly[contactName], 'Contact')
			elif 'proximity'   in contact:		MakeProximityGroup(part, Assembly[contactName], 'Contact')
			elif 'normals'     in contact:		MakeNormalsGroup  (part, Assembly[contactName], 'Contact')
			elif 'contains'    in contact:		MakeContainsGroup (part, Assembly[contactName], 'Contact')
			else:  # if type all, or anything else
				print(    'Contact between {} and {} is of type All.\n'.format(name, contactName))
				log.write('Contact between {} and {} is of type All.\n'.format(name, contactName))
				part.mesh.MakeGroupByIds('{}_@_{}_ContactFaces'.format(name, contactName), SMESH.FACE, part.mesh.GetElementsByType(SMESH.FACE))
			# remove any tied elements from contact surface, even if contact is ALL
			contactGroup, = part.mesh.GetGroupByName('{}_@_{}_ContactFaces'.format(name, contactName))
			print('Original contact surface has {} faces'.format(len(contactGroup.GetListOfID())) )
			#
			if tiedGroups:
				noCommonNodes = True
				if noCommonNodes:
					print('no common nodes')
					nodesOfContact = part.mesh.CreateDimGroup(contactGroup, SMESH.NODE, 'NodesOfContact', SMESH.ALL_NODES, 0)	# contactGroup may need to be a list of one
					unionTieNodes = part.mesh.UnionListOfGroups(tiedNodeGroups, 'UnionOfTieNodes')
					nodesOfContactTieEdge = part.mesh.IntersectGroups(unionTieNodes, nodesOfContact,	'NodesOfContactTieEdge')
					if len(nodesOfContactTieEdge.GetListOfID()) > 0:
						print('Tie and contact overlap {} faces'.format(len(nodesOfContactTieEdge.GetListOfID())) )
						nodesOfCutTrimmed = part.mesh.CutListOfGroups([nodesOfContact], [nodesOfContactTieEdge], '{}_@_{}_ContactNodes'.format(name, contactName))
						if len(nodesOfCutTrimmed.GetListOfID()) >= 3:
							print('reduced contact still has {} faces'.format( len(nodesOfCutTrimmed.GetListOfID()) )  )
							part.mesh.RemoveGroup(contactGroup)
							part.mesh.CreateDimGroup(nodesOfContact, SMESH.FACE,'{}_@_{}_ContactFaces'.format(name, contactName),SMESH.ALL_NODES, 0)
							part.mesh.RemoveGroup(nodesOfContact)
						else:
							part.mesh.RemoveGroup(nodesOfCutTrimmed)
							print(    'WARNING!		No elements remain in contact surface!\n')
							log.write(name)
							log.write('WARNING!		No elements remain in contact surface!\n')
					part.mesh.RemoveGroup(unionTieNodes)
					part.mesh.RemoveGroup(nodesOfContactTieEdge)
				else:
					CutGroups(part.mesh, contactGroup, tiedGroups)
		for otherName, other in part.Other.items():
			print('otherName:', otherName)
			if   otherName == 'Vector':		MakeVectorSurface(part)
			elif otherName == 'Remaining':	MakeRemainingSurface(part)

def MakeAllGroups(mesh, dim='3'):
	"""Make mesh Group of type "All" for specified dimensions
		mesh is mesh data, dim a the dimension of element wanted
		for example dim = '01' would make all groups for nodes and edges
		no return value but mesh is augmented with groups"""
	name = mesh.GetName()
	if   '0' in dim:	mesh.MakeGroupByIds('{}_All_Nodes'.format(name),   SMESH.NODE,   mesh.GetElementsByType(SMESH.NODE))
	elif '1' in dim:	mesh.MakeGroupByIds('{}_All_Edges'.format(name),   SMESH.EDGE,   mesh.GetElementsByType(SMESH.EDGE))
	elif '2' in dim:	mesh.MakeGroupByIds('{}_All_Faces'.format(name),   SMESH.FACE,   mesh.GetElementsByType(SMESH.FACE))
	elif '3' in dim:	mesh.MakeGroupByIds('{}_All_Volumes'.format(name), SMESH.VOLUME, mesh.GetElementsByType(SMESH.VOLUME))
def MakeNormalsGroup(part1, part2, modifier='Ties'):
	"""part1 and part2 are  are SimulationPart defined in ConnectivityXml.py
		modifer is a type name that will be added to the group name for identification
		creates mesh groups on both meshes based on element normal vectors and limited in extent by proximity
		no returm value, but part1.mesh and part2.mesh will be augmented with groups"""
	name1, name2 = part1.name, part2.name
	mesh1, mesh2 = part1.mesh, part2.mesh
	#
	center1 = FindCenter(mesh1)
	center2 = FindCenter(mesh2)
	#
	vector1 = center2 - center1
	#
	multiplier1 = getattr(part1, modifier)[name2]['normals']
	# except KeyError:	multiplier1 = 0.0
	prox1 = FindProximity(part1, part2, multiplier1)
	faces1 = FindOrientedNormals(mesh1, prox1, vector1)
	#
	surface1 = mesh1.MakeGroupByIds('{}_@_{}_{}Faces'.format(name1, name2, modifier), SMESH.FACE, faces1)
	nodes1 =   mesh1.CreateDimGroup( surface1, SMESH.NODE, '{}_@_{}_{}Nodes'.format(name1, name2, modifier), SMESH.ALL_NODES, 0)
	return surface1
def MakeProximityGroup(part1, part2, modifier='Ties'):
	"""part1 and part2 are SimulationPart defined in ConnectivityXml.py
		modifer is a type name that will be added to the group name for identification
		creates mesh groups on part1 based of inter mesh proximity
		no returm value, but part1.mesh will be augmented with group"""
	name1, name2 = part1.name, part2.name
	mesh1 = part1.mesh
	#
	multiplier1 = getattr(part1, modifier)[name2]['proximity']
	prox1 = FindProximity(part1, part2, multiplier1)
	surf1 = mesh1.MakeGroupByIds('{}_@_{}_{}Nodes'.format(name1, name2, modifier), SMESH.NODE, prox1 )
	#
	#	SMESH.AT_LEAST_ONE	- include if one or more node is common
	#	SMESH.MAJORITY		- include if half of nodes or more are common.
	#	SMESH.ALL_NODES		- include if all nodes are common
	mesh1.CreateDimGroup( surf1, SMESH.FACE, '{}_@_{}_{}Faces'.format(name1, name2, modifier), SMESH.MAJORITY, 0)
def MakeContainsGroup(part1, part2, modifier='Ties', threshold=0.0):
	"""Normal Vector based contact for type contains"""
	name1, name2 = part1.name, part2.name
	mesh1, mesh2 = part1.mesh, part2.mesh
	#
	multiplier1, orient1 = getattr(part1, modifier)[name2]['contains']  # inner or outer
	center1 = FindCenter(mesh2)
	#
	#Multis_febio
	prox1 = FindProximity(part1, part2, multiplier1)	# nodes in proximity
	#
	tempNodes = mesh1.MakeGroupByIds('TempNodes', SMESH.NODE, prox1)	# uses prox1 to make a node group
	tempFaces = mesh1.CreateDimGroup(tempNodes, SMESH.FACE, "Temp", SMESH.AT_LEAST_ONE, 0)	# uses tempNodes to make a face group
	# extact faces from face group
	faces1 = tempFaces.GetListOfID()
	#
	contains = []
	for i, face in enumerate(faces1):
		normal = mesh1.GetFaceNormal(face, normalized=True)
		if   orient1 == 'inner':		estimatedNormal = center1-mesh1.BaryCenter(face)
		elif orient1 == 'outer':		estimatedNormal = mesh1.BaryCenter(face)-center1
		# estimatedNormal = estimatedNormal/np.linalg.norm(estimatedNormal)
		dot = np.dot(normal, estimatedNormal)
		# Could find angle with np.arccos(dot)
		if dot > threshold:		contains.append(face)
	#
	surf1 = mesh1.MakeGroupByIds('{}_@_{}_{}Faces'.format(name1, name2, modifier), SMESH.FACE, contains)
	mesh1.CreateDimGroup(surf1, SMESH.NODE, '{}_@_{}_{}Nodes'.format(name1, name2, modifier) )
	#
	# remove temporary groups from mesh
	mesh1.RemoveGroup(tempNodes)
	mesh1.RemoveGroup(tempFaces)
def MakeLocalGroup(part1, part2, modifier='Ties', threshold=0.0):
	"""part1 and part2 are  are SimulationPart defined in ConnectivityXml.py
		modifer is a type name that will be added to the group name for identification
		creates mesh groups on both meshes based on local element normal vectors
		no returm value, but part1.mesh and part2.mesh will be augmented with groups"""
	name1, name2 = part1.name, part2.name
	mesh1, mesh2 = part1.mesh, part2.mesh
	#
	print('MakeLocalGroup()', name1, name2)
	#
	nodes1, array1 = NodesAsArray(part1)
	nodes2, array2 = NodesAsArray(part2)
	#
	elemList = {}
	for n in nodes1:
		node = array1[n - 1]
		elems = mesh1.GetNodeInverseElements(n)
		for i in range(len(elems) - 1, -1, -1):  # iterate backwards so the higher indices don't change mid loop!
			if part1.mesh.GetElementType(elems[i]) == SMESH.VOLUME:
				elems.pop(i)	# remove volumes, only operate on faces
			# elif elems[i] in elemList:
			# 	elems.pop(i)	# removes elements already included in group
		normals = {elem:part1.mesh.GetFaceNormal(elem, normalized=True) for elem in elems}
		average = np.mean(np.array([n for n in normals.values()]))
		print(average)
		if normals:
			norms = np.array([n for n in normals.values()])
			#
			minDiff = np.min(np.linalg.norm(array2 - node, axis=1))
			minIndex, = np.where(minDiff == np.nanmin(minDiff))	# closest node in other mesh
			#
			oppoElems = part2.mesh.GetNodeInverseElements(nodes2[minIndex])	# list of elements including closest node
			oppoNormals = []
			for j in range(len(oppoElems) - 1, -1, -1):  # iterate backwards so the higher indices don't change mid loop!
				if part1.mesh.GetElementType(oppoElems[j]) == SMESH.VOLUME:
					oppoElems.pop(j)  # remove volumes, only operate on faces
			for op in oppoElems:
				oppoNormals.append(part2.mesh.GetFaceNormal(op, normalized=True))
			oppoNormals = np.array(oppoNormals)
			avg = np.mean(oppoNormals, axis=0)  # this should be the average normal of opposing side
			#
			dot = np.dot(norms, avg)
			for i in range(len(dot)):
				if dot[i] >= threshold:
					elemList[elems[i]] = ''  # searching dictionary keys is faster than searching a list for inclusion.
	#
	surface1 = mesh1.MakeGroupByIds('LOCAL{}_@_{}_{}Faces'.format(name1, name2, modifier), SMESH.FACE, elemList.keys())
	nodes1 =   mesh1.CreateDimGroup( surface1, SMESH.NODE, 'LOCAL{}_@_{}_{}Nodes'.format(name1, name2, modifier), SMESH.ALL_NODES, 0)

def MakeNewGroup(part1, part2, modifier='Ties'):
	"""part1 and part2 are  are SimulationPart defined in ConnectivityXml.py
		modifer is a type name that will be added to the group name for identification
		creates mesh groups on both meshes based on local element normal vectors
		no returm value, but part1.mesh and part2.mesh will be augmented with groups"""
	name1, name2 = part1.name, part2.name
	mesh1, mesh2 = part1.mesh, part2.mesh
	#
	center2 = FindCenter(mesh2)
	#
	addToGroup = 0
	try:						threshold = getattr(part1, modifier)[name2]['local']
	except KeyError:
		try:					threshold = getattr(part1, modifier)[name2]['local']
		except KeyError:		threshold = 0.0
	print('threshold', threshold)
	print('MakeNewGroup()', name1, name2)
	#
	nodes1, array1 = NodesAsArray(part1)
	#
	elemList = {}
	for n in nodes1:
		node = array1[n - 1]
		localVector = node-center2		# vector from node to center of other object
		localVector = NormalizeVector(localVector)
		elems = mesh1.GetNodeInverseElements(n)
		for i in range(len(elems) - 1, -1, -1):  # iterate backwards so the higher indices don't change mid loop!
			if part1.mesh.GetElementType(elems[i]) == SMESH.VOLUME:
				elems.pop(i)	# remove volumes, only operate on faces
			elif elems[i] in elemList:
				elems.pop(i)	# removes elements already included in group
		normals = {elem:part1.mesh.GetFaceNormal(elem, normalized=True) for elem in elems}
		if normals:
			norms = np.array([n for n in normals.values()])
			dot = np.dot(norms, localVector)
			#
			for i in range(len(dot)):
				if dot[i] <= threshold:
					elemList[elems[i]] = ''  # searching dictionary keys is faster than searching a list for inclusion.
					addToGroup+=1
	#
	surface1 = mesh1.MakeGroupByIds('NEW{}_@_{}_{}Faces'.format(name1, name2, modifier), SMESH.FACE, elemList.keys())
	nodes1 =   mesh1.CreateDimGroup( surface1, SMESH.NODE, 'NEW{}_@_{}_{}Nodes'.format(name1, name2, modifier), SMESH.ALL_NODES, 0)
	#
	print(name1, name2)
	print('added', addToGroup)

# def NormalizeVector(vector):
# 	v = vector
# 	mag = (v[0] ** 2 + v[1] ** 2 + v[2] ** 2) ** 0.5
# 	vector = np.array([v[0] / mag, v[1] / mag, v[2] / mag])
# 	return vector

def MakeVectorSurface(part):
	"""Make a surface based off a known predefined vector
		This should be used for the occasional
		'I don't know how to define this' surface"""
	name, mesh = part.name, part.mesh
	#
	faces = FindVectorFaces(part)
	surf = mesh.MakeGroupByIds('{}_Vector_Faces'.format(name), SMESH.FACE, faces)
	#
	mesh.CreateDimGroup(surf, SMESH.NODE, '{}_Vector_Nodes'.format(name))
def MakeRemainingSurface(part):
	"""Creates a contact surface that is all surfaces not already on a group
		no return value but mesh is augmented"""
	mesh = part.mesh
	groups = mesh.GetGroups(SMESH.FACE)
	all=[]
	tool=[]
	for g in groups:
		if  g.GetName() == "{}_All_Faces".format(mesh.GetName()):	all.append(g)
		else:														tool.append(g)
	name = mesh.GetName()
	if tool:
		mesh.CutListOfGroups(all, tool, '{}_Remaining_Faces'.format(name))
	else:
		filter = smesh.GetFilter(SMESH.FACE, SMESH.FT_FreeFaces)
		mesh.GroupOnFilter(SMESH.FACE, '{}_Remaining_Faces'.format(name), filter)
def MakePointSurface(part, point):
	"""Make a surface based off a known predefined point
		part is SimulationPart defined in ConnectivityXml.py
		creates mesh groups on part based off mesh proximity to point
		no returm value, but part.mesh will be augmented with group"""
	#
	# multiplier1 = getattr(part1, modifier)[name2]['proximity']
	#
	prox1 = FindProximity(part, point, multiplier=1.0)		# issue here as point is not a part!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	surf1 = part.mesh.MakeGroupByIds('{}_Point_Nodes'.format(part.name), SMESH.NODE, prox1)
	#
	#	SMESH.AT_LEAST_ONE	- include if one or more node is common
	#	SMESH.MAJORITY		- include if half of nodes or more are common.
	#	SMESH.ALL_NODES		- include if all nodes are common
	part.mesh.CreateDimGroup(surf1, SMESH.FACE, '{}_Point_Faces'.format(part.name), SMESH.AT_LEAST_ONE, 0)
def MakePlaneSurface(part):
	"""Make a surface based off a plane
		I am not sure of a use case for this method
		so it isn't implemented"""
	pass

def MakeMesh(meshFile, material='elastic', meshName='', ifSurface=True, ifAll=True, ifQuad=None):
	"""Takes an stl mesh file
		open file,
		mesh in 3D,
		set all groups that require no other information.
		returns mesh data"""
	mesh = OpenMesh(meshFile)
	if meshFile.endswith('.stl'):
		if material != 'rigid':
			meshMin, meshMax = mesh.GetMinMax(SMESH.FT_Length2D)
			print('	optimizing.');	log.write('	optimizing.\n')
			#
			#	3D Meshing and parameters
			#
			NETGEN_3D = mesh.Tetrahedron(geom=None)
			NETGEN_3D_Parameters = NETGEN_3D.Parameters()
			NETGEN_3D_Parameters.SetMaxSize(meshMax)
			NETGEN_3D_Parameters.SetOptimize(1)
			NETGEN_3D_Parameters.SetMinSize(meshMin)
			NETGEN_3D_Parameters.SetFineness(4)
			#
			if not mesh.Compute():		# this fixes an error I had with one of my test cases
				print('\n\n\nChanging NETGEN parameters and remeshing\n\n\n')
				NETGEN_3D_Parameters.SetFineness(3)		# I don't know why changing the fineness matters
				mesh.Compute()
			if ifQuad is not None:				mesh.ConvertToQuadratic(1)
	elif meshFile.endswith('.unv'):
		print('\n########################################################################\n')
		dim = mesh.MeshDimension()
		print(dim)
		pass    # I don't think I want to mesh if it still 2D...hmmm
	elif meshFile.endswith('.med'):
		mesh=mesh[0]
		groups = mesh.GetGroups(SMESH.ALL)
		print('	removing existing groups.');	log.write('	removing existing groups.\n')
		for g in groups:
			info = smesh.GetMeshInfo(g)
			keys = info.keys();	keys.sort()
			volumeGroup = False
			for key in keys:
				if     ('Quad'  in str(key) and info[key] != 0)\
					or ('Tetra' in str(key) and info[key] != 0)\
					or ('Penta' in str(key) and info[key] != 0):	volumeGroup = True
			if not volumeGroup:										mesh.RemoveGroup(g)
	#
	if meshName:	mesh.SetName(meshName)
	else:			mesh.SetName(mesh.GetName().split('.')[0])	  # remove '.stl' from file name
	#
	dim = mesh.MeshDimension()
	if ifAll:					MakeAllGroups(mesh, str(dim) )
	if ifSurface and dim != 2:	MakeAllGroups(mesh, '2')
	return mesh

def NodesAsArray(part):
	"""Returns list of node IDs and numpy array of all nodes in the mesh
		part is SimulationPart defined in ConnectivityXml.py
			the array can be indexed by (node-1)"""
	from ConnectivityXml import SimulationPart
	if isinstance(part, SimulationPart):
		nodes = NodesInGroup(part.mesh)
		try:	AllNodeArray = MEDLoader.ReadUMeshFromFile(part.med, 0).getCoords().toNumPyArray()
		except AttributeError:
			coords = MEDLoader.ReadUMeshFromFile(part.med, 0).getCoords()
			AllNodeArray = np.zeros((len(coords),3))
			for i,node in enumerate(coords):
				try:	AllNodeArray[i] = node
				except ValueError:
					node = [n for n in node].append(0.0)
					AllNodeArray[i] = node
	elif isinstance(part, np.ndarray):
		nodes = part.reshape(-1, 3).shape[0]
		AllNodeArray = part
	else:	print('NodesAsArray() does not accept type', type(part))
	
	# print(part.name)
	# print('nuber of nodes:',len(nodes), '\nNumbr in array:',len(AllNodeArray))
	# print(AllNodeArray)
	# print(nodes)
	# for i in range(len(nodes)+10):
	# 	if i not in nodes:
	# 		print(i)
	
	return nodes, AllNodeArray

def NodesInGroup(mesh, groupName=None):
	"""mesh is mesh data
		returns list of node IDs and numpy array with dim nodes * 3"""
	if groupName:
		group, = mesh.GetGroupByName(groupName)
		nodes = group.GetListOfID()
	else:
		filter = smesh.GetFilter(SMESH.FACE, SMESH.FT_FreeFaces)
		tempSurface = mesh.GroupOnFilter(SMESH.FACE, "ActualSurface", filter)
		tempNodes = mesh.CreateDimGroup(tempSurface, SMESH.NODE, "TempNodes", SMESH.AT_LEAST_ONE, 0)
		nodes = tempNodes.GetListOfID()
		#
		mesh.RemoveGroup(tempSurface)
		mesh.RemoveGroup(tempNodes)
	return nodes

def OpenMesh(meshFile):
	"""Opens mesh file, creates mesh data"""
	print('Opening {}'.format(meshFile))
	try:	log.write('Opening {}\n'.format(meshFile))
	except NameError:	pass
	if   meshFile.endswith('.stl'):     mesh = smesh.CreateMeshesFromSTL(meshFile)
	elif meshFile.endswith('.med'):     mesh, status = smesh.CreateMeshesFromMED(meshFile)
	elif meshFile.endswith('.unv'):
		print('elif endswith unv')
		mesh = smesh.CreateMeshesFromUNV(meshFile)
		print('ending elif unv')
	else:
		print('Mesh type {} not handled by OpenMesh'.format(meshFile.split['.'][-1]))
		try:	log.write('Mesh type {} not handled by OpenMesh\n'.format(meshFile.split['.'][-1]))
		except NameError:	pass
		mesh = None
	return mesh

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(os.path.realpath(connectivityFile)))
	print('Reading from connectivity file {}'.format(connectivityFile))
	log.write('Reading from connectivity file {}\n'.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)
	except IOError:
		print('\n\nNo connectivity file given.\nUsage: salome StlToMed.py args:ConnectivityFile\n')
		sys.exit()
	#
	for part in Assembly.values():  # Assert that files exist
		if not os.path.isfile(part.file):
			print('\n\nScript cannot find file {}\n	associated with {}'.format(part.file, part.name))
			print('Please check connectivity file.\n')
			log.write('Script cannot find file {}\n	associated with {}\n'.format(part.file, part.name))
			log.write('Please check connectivity file.\n')
			#CloseSalome()
			sys.exit()
	#
	# now becasue everything is in place, make sure the MED directory exists
	medDir = os.path.join(os.path.dirname(connectivityFile), 'MED')
	if not os.path.isdir(medDir):    os.mkdir(medDir)
	#
	return Assembly

def Logfile(logName='Connectivity.xml'):
	try:				log = open(logName.replace('xml', 'log'), 'w')
	except IndexError:	log = open('Connectivity.log', 'w')
	log.write('Logfile created by {}\n'.format(inspect.getfile(inspect.currentframe())))
	log.write('	on {}\n\n'.format(time.strftime('%b %d %Y at %H:%m:%S %p')))
	return log

def _runTest(connectivityFile='Connectivity.xml', close='False'):
	"""Usually similar to Make Meds
	 with stuff i'm testing and more printed output"""
	print(sys.argv)
	print('\n\n---TESTING---\n\nMade with _runTest() function NOT running as expected!\n')
	log.write('---TESTING---\n\nMade with _runTest() function NOT running as expected!\n')
	Assembly = ReadConnectivity(connectivityFile)
	#
	medDir = os.path.join(os.path.dirname(connectivityFile), 'MED')
	if not os.path.isdir(medDir):    os.mkdir(medDir)
	#
	for part in Assembly.values():
		part.mesh=MakeMesh(part.file, part.mat, part.name, ifQuad=part.quad)
		if part.file.endswith('.stl') or part.file.endswith('.unv'):
			if os.path.isfile(part.med):    os.remove(part.med)  # seems to be necessary on Windows
			part.mesh.ExportMED(part.med)	# write it out so I can access nodes with MEDLoader, because its much faster
	#
	MakeGroups(Assembly)
	#
	for part in Assembly.values():
		if os.path.isfile(part.med):    os.remove(part.med)  # seems to be necessary on Windows
		part.mesh.ExportMED(part.med)	# write med with all groups added
	print('All done exporting.')
	#
	if eval(close):		CloseSalome()  # eval is necessary because command line arguments are passed as strings
	else: 			print('NOT closing Salome!!!')
	print('\n\n---TESTING---\n\nMade with _runTest() function NOT running as expected!\n')
	log.write('---TESTING---\n\nMade with _runTest() function NOT running as expected!\n')


def MakeMeds(connectivityFile='Connectivity.xml', close='True'):
	"""Read connectivity file, then stls and makes MED files and groups"""
	Assembly = ReadConnectivity(connectivityFile)
	#
	medDir = os.path.join(os.path.dirname(connectivityFile), 'MED')
	if not os.path.isdir(medDir):    os.mkdir(medDir)
	#
	for part in Assembly.values():
		part.mesh=MakeMesh(part.file, part.mat, part.name, ifQuad=part.quad)
		if part.file.endswith('.stl') or part.file.endswith('.unv'):
			if os.path.isfile(part.med):	os.remove(part.med)	#seems to be necessary on Windows
			part.mesh.ExportMED(part.med)	# write it out so I can access nodes with MEDLoader, because its much faster
	MakeGroups(Assembly)
	#
	for part in Assembly.values():
		if os.path.isfile(part.med):	os.remove(part.med)	#seems to be necessary on Windows
		part.mesh.ExportMED(part.med)	# write med with all groups added
	print('All done exporting.');	log.write('All done exporting.\n')
	#
	print('Script finished.')
	if eval(close):		CloseSalome()		# eval is necessary because command line arguments are passed as strings
	else: 				print('NOT closing Salome!!!')


if __name__ == '__main__':
	sys.path.append(os.path.dirname(sys.argv[0]))
	try:				log = Logfile(sys.argv[1])
	except IndexError:	log = Logfile()
	MakeMeds(*sys.argv[1:])
	# _runTest(*sys.argv[1:])