#!/usr/bin/env python

#import xml.etree.ElementTree as ET
import sys
import string

class FebioElement:
	def __init__(self, material, name):
		self.mat = material
		self.name = name
		self.type = None
		self.connect = {}		# contains matched pairs of id# and element connectivity
		self.nodes=set()

	def __repr__(self):
		return "<FebioElement type:{} mat:{} connect:{} elements>".format(self.type, self.mat, len(self.connect))


gmshMapping={\
'1':'line',\
'2':'tri3',\
'3':'quad4',\
'4':'tet4',\
'5':'hex8',\
'8':'line3',\
'9':'tri6',\
'11':'tet10',\
'12':'hex27',\
'15':'point',\
'16':'quad8',\
'17':'hex20',\
# gmsh and febio second order tetrahedra switch the ordering of nodes 9 and 10
}
#for g in gmshMapping: print g, gmshMapping[g]

nodeData={}
elemData={}
physGroups={}	# physGroups{PG#] = "name in mesh"

# process msh data from sys.arg[1]
#
#	Could loop over all sys.args for multiple meshes
# 	 But leads to a problem with node and element numbering !
#
meshFile = open(sys.argv[1], 'r')
print 'Parsing Mesh {}'.format(sys.argv[1])

line = meshFile.readline()
while line:
	if line.startswith('$Nodes') or line.startswith('$NOD'):
		print 'Reading Nodes...'
		nodes = int( meshFile.readline() )
		for i in xrange(nodes):
			line=meshFile.readline().split()	#readline and split by whitespace
			nodeData[int(line[0])]= tuple(float(i) for i in line[-3:])
	elif line.startswith('$Elements'):
		print 'Reading Elements...'
		elements = int( meshFile.readline() )
		for i in xrange(elements):
			line=meshFile.readline().split()	#readline and split by whitespace
			mapping = gmshMapping[line[1]]
			if   mapping.startswith('line'):	print 'Lines are currently ignored.'
			elif mapping.startswith('quad'):	print 'Quadrangles are currently ignored.'
			elif mapping.startswith('hex'):	print 'Hexahedrons are currently ignored.'
			elif mapping.startswith('point'):	print 'Points are currently ignored.'
			elif mapping.startswith('t'):	# tet and tri elements	 # this code could be used for any element type being converted
				nodeNum = int(mapping.strip(string.ascii_letters))	# number of nodes in the element
				elemData[line[3]].connect[int(line[0])] = [int(j) for j in line[-nodeNum:]]
				elemData[line[3]].type=mapping
				if mapping.startswith('tri'):	# triangles only but would be used with all surface meshes
					elemData[line[3]].nodes.update( int(j) for j in line[-nodeNum:] ) 
				elif mapping.startswith('tet10'):
					c = elemData[line[3]].connect[int(line[0])]
					c[-2],c[-1] = c[-1],c[-2]		# swap last two elements in connectivity for tet10
			else:				print 'Elemenet type {} unaccounted for'.format(line[1])
	elif line.startswith('$ELM'):
		print 'Reading ELeMents...'
		elements = int( meshFile.readline() )
		for i in xrange(elements):
			line=meshFile.readline().split()	#readline and split by whitespace
			mapping = gmshMapping[line[1]]
			if   mapping.startswith('line'):	print 'Lines are currently ignored.'
			elif mapping.startswith('quad'):	print 'Quadrangles are currently ignored.'
			elif mapping.startswith('hex'):	print 'Hexahedrons are currently ignored.'
			elif mapping.startswith('point'):	print 'Points are currently ignored.'
			elif mapping.startswith('t'):	# tet and tri elements	 # this code could be used for any element type being converted
				nodeNum = int(mapping.strip(string.ascii_letters))	# number of nodes in the element
				try:	elemData[line[2]].connect[int(line[0])] = [int(j) for j in line[-nodeNum:]]
				except KeyError:
					s=''
					if nodeNum == 3:	s='_surf_'+line[2]
					elemData[line[2]] = FebioElement(line[1], sys.argv[1].split('.')[0]+s)
					elemData[line[2]].connect[int(line[0])] = [int(j) for j in line[-nodeNum:]]
				elemData[line[2]].type=mapping
				if mapping.startswith('tri'):	# triangles only but would be used with all surface meshes
					elemData[line[2]].nodes.update( int(j) for j in line[-nodeNum:] ) 
				elif mapping.startswith('tet10'):
					c = elemData[line[2]].connect[int(line[0])]
					c[-2],c[-1] = c[-1],c[-2]		# swap last two elements in connectivity for tet10
			else:				print 'Elemenet type {} unaccounted for'.format(line[1])
	elif line.startswith('$PhysicalNames'):
		names = int( meshFile.readline() )
		for i in xrange(names):
			line=meshFile.readline().split()
			# line[0]	dimensionality of group as a string
			# line[1]	physical group # as a string
			# line[2]	name as a string already in quotes
			physGroups[line[1]] = line[2]
			#   Populate FebioElementData here
			elemData[line[1]] = FebioElement(line[1], line[2])
	line=meshFile.readline()


for g in physGroups: print g, physGroups[g]

outFile = sys.argv[1].split('.')[0]+'.feb'
output = open(outFile, 'w')
print 'Converting...'

output.write('<?xml version="1.0" encoding="ISO-8859-1"?>\n')
import time
output.write('<!-- This document was written by\n	{}\n	on {}-->\n\n\n'.format(sys.argv[0], time.strftime('%b %d %Y at %H:%m:%S %p') ))
output.write('<febio_spec version="2.0">\n')

# Material
output.write('	<Material>\n')
#
output.write('		<material id="1" name="Rigid" type="rigid body">\n')
#~ output.write('			<density>2.8e-09</density>\n')
output.write('			<density>7e-09</density>\n')
output.write('		</material>\n')
#
#~ output.write('		<material id="2" name="Muscle" type="Ogden">\n')
#~ output.write('			<density>1e-09</density>\n')
#~ output.write('			<k>1.128</k>\n')
#~ output.write('			<c1>0.05737</c1>\n')
#~ output.write('			<m1>28.82</m1>\n')
#~ output.write('		</material>\n')
#~ #
#~ output.write('		<material id="3" name="Fat" type="Ogden">\n')
#~ output.write('			<density>1e-09</density>\n')
#~ output.write('			<k>0.3251</k>\n')
#~ output.write('			<c1>0.01653</c1>\n')
#~ output.write('			<m1>-0.77</m1>\n')
#~ output.write('		</material>\n')
#~ #
output.write('		<material id="2" name="Skin" type="Ogden">\n')
output.write('			<density>1e-09</density>\n')
output.write('			<k>0.7979</k>\n')
output.write('			<c1>0.04057</c1>\n')
output.write('			<m1>22.71</m1>\n')
output.write('		</material>\n')
#
output.write('	</Material>\n')

#Geometry
output.write('	<Geometry>\n')
output.write('		<Nodes>\n')	# Nodes
#~ for n in nodeData:	output.write('			<node id="{}"> {: 8.5f}, {: 8.5f}, {: 8.5f} </node>\n'.format(n, *nodeData[n]) )
for n in nodeData:	output.write('			<node id="{}"> {: 10.7e}, {: 10.7e}, {: 10.7e} </node>\n'.format(n, *nodeData[n]) )
output.write('		</Nodes>\n')
for e in sorted(elemData.itervalues()):		# Elements and Surfaces
	if not e.type:raw_input('None type found')
	if e.type.startswith('tri'):	#  or e.type.startswith('quad'):
		if '_surf' in e.name:
			output.write('		<Surface name={}>\n'.format(e.name))
			for id, c in sorted(e.connect.iteritems()):	output.write('			<{} id="{}"> {} </{}>\n'.format( e.type, id, str(c)[1:-1], e.type ))
			output.write('		</Surface>\n')
		#
		# # 	Not logic, Just a shortcut!
		#
		if 'bone' in e.name or 'Symmetry' in e.name or 'Fixed' in e.name:
			output.write('		<NodeSet name={}>\n'.format(e.name))
			for id in sorted(e.nodes):	output.write('			<node id="{}"/>\n'.format(id) )
			output.write('		</NodeSet>\n')
	else:
		output.write('		<Elements type="{}" mat="{}" name={}>\n'.format(e.type, e.mat, e.name) )
		# FeBio does not require different element to have unique elem ids
		for id, c in sorted(e.connect.iteritems()):	output.write('			<elem id="{}"> {} </elem>\n'.format(id,str(c)[1:-1]) )
		output.write('		</Elements>\n')
# ElementData?
output.write('	</Geometry>\n')

#Boundary
#	Repeat fix as necessary (options: x, y, z displacement and u, v, w rotation on x, y, z)
output.write('	<Boundary>\n')
for g in physGroups:
	if 'bone' in physGroups[g]: 
		output.write('		<fix bc="xyz" set={}/>\n'.format(physGroups[g]) )
	elif 'Symmetry' in physGroups[g]:
		output.write('		<fix bc="{}" set={}/>\n'.format(physGroups[g][1], physGroups[g]) )
	elif 'Fixed' in physGroups[g]:
		output.write('		<fix bc={} set={}/>\n'.format(physGroups[g].replace('Fixed',''), physGroups[g]) )
output.write('	</Boundary>\n')


#Initial Constraints
output.write('	<Constraints>\n')
output.write('		<rigid_body mat="1">\n')	# REQUIRES some logic, also in Step (below)
output.write('			<fixed bc="x"/>\n')
output.write('			<fixed bc="y"/>\n')
#~ output.write('			<fixed bc="z"/>\n')
output.write('			<fixed bc="Rx"/>\n')
output.write('			<fixed bc="Ry"/>\n')
output.write('			<fixed bc="Rz"/>\n')
output.write('		</rigid_body>\n')
output.write('	</Constraints>\n')

#Load Data
output.write('	<LoadData>\n')
output.write('		<loadcurve id="1" type="smooth">\n')
output.write('			<point>0,0</point>\n')
output.write('			<point>1,1</point>\n')
output.write('		</loadcurve>\n')
#
output.write('		<loadcurve id="2" type="step">\n')
steps=50
for p in xrange(steps+1):
	output.write('			<point>{:4.2f},{}</point>\n'.format(float(p)/steps,1.0/steps) )
#~ exit()
output.write('		</loadcurve>\n')
#
output.write('	</LoadData>\n')

#Output requirements
output.write('	<Output>\n')
output.write('		<plotfile type="febio">\n')
output.write('			<var type="contact gap"/>\n')
output.write('			<var type="contact pressure"/>\n')
output.write('			<var type="contact traction"/>\n')
output.write('			<var type="displacement"/>\n')
output.write('			<var type="reaction forces"/>\n')
output.write('			<var type="stress"/>\n')
output.write('		</plotfile>\n')
output.write('	</Output>\n')


# Step
output.write('	<Step name="Step01">\n')
#
output.write('		<Module type="solid"/>\n')
#
output.write('		<Control>\n')
output.write('			<time_steps>{}</time_steps>\n'.format(steps))		# required tag
output.write('			<step_size>{}</step_size>\n'.format(1.0/steps))	# required tag
output.write('			<time_stepper>\n')
output.write('				<dtmin>1e-10</dtmin>\n')
output.write('				<dtmax lc="2">0.05</dtmax>\n')		# Load curve necessary for MustPoint printing
output.write('				<max_retries>25</max_retries>\n')
output.write('				<opt_iter>10</opt_iter>\n')
output.write('				<aggressiveness>1</aggressiveness>\n')
output.write('			</time_stepper>\n')
output.write('			<optimize_bw>1</optimize_bw>\n')
output.write('			<plot_level>PLOT_MUST_POINTS</plot_level>\n')
output.write('			<analysis type="static"/>\n')
output.write('		</Control>\n')
#
output.write('		<Contact>\n')
output.write('			<contact type="facet-to-facet sliding" name="Indentor Skin">\n')
output.write('				<laugon>1</laugon>\n')
output.write('				<penalty>10</penalty>\n')
output.write('				<two_pass>0</two_pass>\n')
output.write('				<maxaug>10</maxaug>\n')
output.write('				<surface type="master" set="Indentor_surf"/>\n')	# How do I choose my surfaces?
output.write('				<surface type="slave"  set="skin_surf"/>\n')
output.write('			</contact>\n')
output.write('		</Contact>\n')
#
output.write('		<Constraints>\n')
output.write('			<rigid_body mat="1">\n')	# REQUIRES SOME LOGIC!!!
#~ output.write('				<prescribed bc="x" lc="1">-20</prescribed>\n')
#~ output.write('				<prescribed bc="y" lc="1">-100</prescribed>\n')
output.write('				<prescribed bc="z" lc="1">-5</prescribed>\n')
output.write('			</rigid_body>\n')
output.write('		</Constraints>\n')
#
output.write('	</Step>\n')


output.write('</febio_spec>\n')
output.close()

print 'FEBio XML written to', outFile