#=============================================================================
# Monolithic FE^2
# Nils Lange, Geralf Huetter, Bjoern Kiefer
#   Nils.Lange@imfd.tu-freiberg.de, Geralf.Huetter@imfd.tu-freiberg.de, 
#   Bjoern.Kiefer@imfd.tu-freiberg.de
# distributed under CC BY-NC-SA 4.0 license
# (https://creativecommons.org/licenses/by-nc-sa/4.0/)
# Reference: 
#   N. Lange, G. Huetter, B. Kiefer: "An efficient monolithic solution scheme for FE2 problems",
#   https://arxiv.org/abs/2101.01802
#
# Further information on the implementation, structure of the source code,
# examples and tutorials can be found in the file doc/documentation.pdf
#
#=============================================================================

from abaqus import *	#import abaqus functionalities
from abaqusConstants import *
import regionToolset
import os	#import os functionalities to get current directory path

def FE2_Inputgenerator(Modelname,Partname,jobname,rve_number,
			information=" ",upside_abq=[],downside_abq=[],
			rightside_abq=[],leftside_abq=[],gen_constraints_CAE=False,
			material2_specifier='none',material_2_picked=[],
			material_2_set='',PBC='FoamGUI'):

	if mdb.models[Modelname].parts[Partname].space.name=='TWO_D_PLANAR':
		dimens=2
	elif mdb.models[Modelname].parts[Partname].space.name=='THREE_D':
		dimens=3
	else:
		raise Exception('Only 2D planar or 3D Parts allowed!')

	
	#-> get mesh coordinates and element to node connection
	#call get_data_from_abaqus
	(COORDSGLOBAL, element_to_node, nNodes, nElem, nnodes_e,
	 ndof_e,ndof_n,ndof_macro,JTYPE)= get_data_from_abaqus(Modelname,Partname,dimens,PBC)


	#-> get the connection element to materialnumber
	#call get_material_to_element
	element_to_material\
	= get_material_to_element(nElem, material2_specifier,Modelname,Partname,
	                          material_2_picked,material_2_set)

												
	#-> get the connection between slave dof's and master dof's
	#call get_slave_to_master_connnection
	if dimens==2:
		(slaves, slaves_dict, slave_to_master_dict, boundarypoint, nslaves,
		delta_x, global_node_to_slave_node)\
		= get_slave_to_master_connnection2D(upside_abq, downside_abq,
			leftside_abq, rightside_abq, COORDSGLOBAL, nNodes, dimens,
			Modelname, Partname, gen_constraints_CAE)
	else:
		(slaves, slaves_dict, slave_to_master_dict, boundarypoint, nslaves,
		delta_x, global_node_to_slave_node)\
		=get_slave_to_master_connnection3D(Modelname,dimens,COORDSGLOBAL,nNodes,PBC)
				
				
	#->	get connection between elements and dof's
	#call get_element_to_dof
	(master_reduced_to_global,node_label_to_red_dof,element_to_global_dof,
	 element_to_reduced_dof,global_master_reduced_to_slave)\
	 = get_element_to_dof(boundarypoint, nNodes, nElem, nnodes_e, dimens,
	   element_to_node, slave_to_master_dict, slaves)
	
	
	#-> get connection element stiffness and k_matrix_values entries
	#call get_stiffnessmatrix_entries for symmetric matrix
	(rowIndex_reduced, columnIndex_reduced,values_to_global_reduced)\
	 = get_stiffnessmatrix_entries(nElem, ndof_e, element_to_reduced_dof, False)
	 
	 #-> get connection element stiffness and k_matrix_values entries
	#call get_stiffnessmatrix_entries for unsymmetric matrix (struct. symmetric)
	(rowIndex_reduced_UNSYMM, columnIndex_reduced_UNSYMM,
	 values_to_global_reduced_UNSYMM)\
	 = get_stiffnessmatrix_entries(nElem, ndof_e, element_to_reduced_dof, True)
		
	
	#-> get "equations" for enforcing the constraint
	#call get equations
	equations\
	= get_equations(dimens,slave_to_master_dict,COORDSGLOBAL,slaves)	
				
				
	#-> create Inputscript, which is read by the FE2 Program
	#call write_data_to_file
	write_data_to_file(jobname,rve_number, dimens,nnodes_e, ndof_e,ndof_n, 
		nElem, nNodes, columnIndex_reduced, columnIndex_reduced_UNSYMM,
		rowIndex_reduced, rowIndex_reduced_UNSYMM, nslaves, ndof_macro, delta_x, COORDSGLOBAL,
		master_reduced_to_global, element_to_global_dof,
		element_to_reduced_dof, values_to_global_reduced,values_to_global_reduced_UNSYMM,
		global_master_reduced_to_slave, equations, element_to_material, 
		global_node_to_slave_node,JTYPE,information)
		
def get_data_from_abaqus(Modelname, Partname, dimens,PBC):
	
###################Coordinates##########################################
	
	COORDSGLOBAL=[]	#coords of the mesh [[c11,c12],[c21,c22],[c31,c32],.
					#element and node labels start counting from
					#(as FORTRAN) but python is counting from 0!
					
	NODES=mdb.models[Modelname].parts[Partname].nodes	#get node infos
	
	for n in range(len(NODES)): #loop through nodes to get coordinates
		if dimens==2:
			COORDSGLOBAL.append([NODES[n].coordinates[0],
								NODES[n].coordinates[1]])
		if dimens==3:
			COORDSGLOBAL.append([NODES[n].coordinates[0],
								NODES[n].coordinates[1],
								NODES[n].coordinates[2]])
								
	if (dimens==3 and PBC=='FoamGUI'):		#last 3 Nodes are only reference
											#points therefore in this
											#implementation not needed
		COORDSGLOBAL.pop(len(COORDSGLOBAL)-1)
		COORDSGLOBAL.pop(len(COORDSGLOBAL)-1)
		COORDSGLOBAL.pop(len(COORDSGLOBAL)-1)
																
	nNodes=len(COORDSGLOBAL)	#number of nodes of the RVE mesh
	
###################Elements#############################################
	
	ELEMENTS=mdb.models[Modelname].parts[Partname].elements
					#get information of all Elements
	
	nElem=len(ELEMENTS)	#number of elements of the RVE Mesh
	
	nnodes_e=len(ELEMENTS[0].connectivity)
					#number of nodes of each element (3=triangular,
					#4=rectangular etc.) -> assumption: all elements
					#have the same amount of nodes!!!!!!!!!!!!!!!
					
	ndof_n=dimens #number of nodes per node (==dimens for purel mechanical problems)
	
	ndof_e=ndof_n*nnodes_e	#degrees of freedom of each Element
	
	ndof_macro=dimens*dimens	#number of macroscopic degrees of freedom (e.g in
								#purely mechanical problems its 9 entries of the
								#1.Piola Kirchhoff stress tensor resp. the
								#displacement gradient)
	
############################get elementtype JTYPE###########################

	#elementlibrary={(Abaqus Elementtype):JTYPE}
	elementlibrary={'CPE4':2001,
					'CPE8':2002,
					'CPS4':2003,
					'C3D8':2004,
					'C3D20':2005,
					'CPE4R':2006,
					'CPE8R':2007,
					'CPS4R':2008,
					'CPS8':2009,
					'CPS8R':2010,
					'C3D20R':2011,
					'CPE3':2012,
					'C3D4':2013,
					'C3D10':2014,
					'CPS3':2024}
					
	JTYPE_List=list(set(ele.type.name for ele in ELEMENTS))
	
	if len(JTYPE_List)>1:
		raise Exception('All elements in the mesh need to have the same element type!')
					
	if JTYPE_List[0] in elementlibrary:
		#get element type
		JTYPE=elementlibrary[JTYPE_List[0]]
	else:
		#display error
		raise Exception('this elementtype does not exist in the UELlib!')
	
###################get element-to node connection###########################
			
	element_to_node=[] #connectivity [[etn11,etn12,etn13],[etn21,etn22..
	
	for e in range(len(ELEMENTS)):	#get element_to_node connectivity
		A=[]						#from ELEMENTES
		for i in range(nnodes_e):
			A.append(ELEMENTS[e].connectivity[i]+1)
		element_to_node.append(A)
		
	return (COORDSGLOBAL,element_to_node,nNodes,nElem,nnodes_e,
			ndof_e,ndof_n,ndof_macro,JTYPE)
				
				
def get_material_to_element(nElem, material2_specifier,Modelname,Partname,
	                          material_2_picked,material_2_set):
	#get material number for each element; currently only 2 Materials
	#possible

	#element_to_material contains connection element<->material

	element_to_material=[1]*nElem #set all element to material connections to 1
		
	if material2_specifier=='by picking elements':
		for i in material_2_picked:
			element_to_material[i.label-1]=2
	if material2_specifier=='by specifying a set':
		A=mdb.models[Modelname].parts[Partname].sets[material_2_set].elements
		for i in A:
			element_to_material[i.label-1]=2
			
	return element_to_material
		
def get_slave_to_master_connnection2D(upside_abq,downside_abq,leftside_abq, 
									rightside_abq,COORDSGLOBAL,nNodes,
									dimens,Modelname,Partname,gen_constraints_CAE):
	
	###################get node labels of RVE edges#####################
	
	upside=[] #list of nodes with node labels of set upside
	for j in range(len(upside_abq)):
		upside.append(upside_abq[j].label)
	upside=list(set(upside)) #remove potential double entries
	
	downside=[] #list of nodes with node labels of set upside
	for j in range(len(downside_abq)):
		downside.append(downside_abq[j].label)
	downside=list(set(downside)) #remove potential double entries

	leftside=[] #list of nodes with node labels of set upside
	for j in range(len(leftside_abq)):
		leftside.append(leftside_abq[j].label)
	leftside=list(set(leftside)) #remove potential double entries
	
	rightside=[] #list of nodes with node labels of set upside
	for j in range(len(rightside_abq)):
		rightside.append(rightside_abq[j].label)
	rightside=list(set(rightside)) #remove potential double entries
	
	if len(upside)!=len(downside) or len(leftside)!=len(rightside):
		raise Exception('Meshes must be congruent!')
	
	##############get sorted node label lists of the edges##############
	
	#sort the node labels of the edge nodelable sets (upside, donwside,
	#leftside,rightside) by x1 resp. x2 in ascending order
	
	#------------------------------------------------------------------#
	tup_leftside=[] #contains tupels of (nodelabel, x2)
	for i in range(len(leftside)):
		tup_leftside.append((leftside[i],COORDSGLOBAL[leftside[i]-1][1]))
	tup_leftside.sort(key=lambda tup: tup[1]) #sort
	#------------------------------------------------------------------#
	tup_rightside=[] #contains tupels of (nodelabel, x2)
	for i in range(len(rightside)):
		tup_rightside.append((rightside[i],COORDSGLOBAL[rightside[i]-1][1]))
	tup_rightside.sort(key=lambda tup: tup[1]) #sort
	#------------------------------------------------------------------#
	tup_downside=[]	#contains tupels of (nodelabel, x1)
	for i in range(len(downside)):
		tup_downside.append((downside[i],COORDSGLOBAL[downside[i]-1][0]))
	tup_downside.sort(key=lambda tup: tup[1]) #sort	
	#------------------------------------------------------------------#
	tup_upside=[] #contains tupels of (nodelabel, x1)
	for i in range(len(upside)):
		tup_upside.append((upside[i],COORDSGLOBAL[upside[i]-1][0]))
	tup_upside.sort(key=lambda tup: tup[1]) #sort
	
	###################detect corner nodes##############################
	
	leftdown=tup_downside[0][0]							#leftdown
	leftup=tup_upside[0][0] 							#leftup
	rightdown=tup_downside[len(tup_downside)-1][0] 		#rightdown
	rightup=tup_upside[len(tup_upside)-1][0] 			#rightup
	
	#------------------------------------------------------------------#
	#if the corner nodes also appear in the lists tup_leftside and
	#tup_rightside, they are to be removed and then treaded seperatly
	if tup_leftside[0][0]==leftdown:
				
		tup_leftside.pop(0)
		tup_leftside.pop(len(tup_leftside)-1)
		tup_rightside.pop(0)
		tup_rightside.pop(len(tup_rightside)-1)
		tup_downside.pop(0)
		tup_downside.pop(len(tup_downside)-1)
		tup_upside.pop(len(tup_upside)-1)
		tup_upside.pop(0)
		
		slaves=[rightup,rightdown] #contains all slave node labels
		
		slave_to_master_dict={} #slave<->master connection
		slave_to_master_dict[rightup]=leftdown	#corner nodes
		slave_to_master_dict[rightdown]=leftup
		
		#displacements boundary condition of this point is zero ->
		#lock rigid body motion
	
		boundarypoint=leftup
	
	else:
				
		slaves=[] #contains all slave node labels
		slave_to_master_dict={} #slave<->master connection
		
		#displacements boundary condition of this point is zero ->
		#lock rigid body motion
	
		boundarypoint=tup_leftside[len(tup_leftside)-1]


	######################sort slaves in lists##########################
	
	#slaves contains all slave node labels
	
	for t_r in tup_rightside:
		slaves.append(t_r[0])
	
	for t_u in tup_upside:
		slaves.append(t_u[0])
	
	slaves.sort()	#sort the slaves list
	
	slaves_dict={} #connection slave node label <-> global node label
	for i in range(len(slaves)):
		slaves_dict[slaves[i]]=i+1
	
	################get slave to master connection######################
		
	 #each key in the dictionary is a global slave node label and has
	 #a value which is the corresponding global master node label
	
	for i in range(len(tup_leftside)):	#edge nodes
		slave_to_master_dict[tup_rightside[i][0]]=tup_leftside[i][0]
		
	for i in range(len(tup_downside)):  #edge nodes
		slave_to_master_dict[tup_upside[i][0]]=tup_downside[i][0]
			
	################additional information##############################
	
	nslaves=len(slave_to_master_dict) #number of slave nodes
	
	delta_x=[0]*dimens #length and width of the RVE
	
	delta_x[0]=COORDSGLOBAL[tup_rightside[0][0]-1][0]-\
			   COORDSGLOBAL[tup_leftside[0][0]-1][0]
	delta_x[1]=COORDSGLOBAL[tup_upside[0][0]-1][1]-\
			   COORDSGLOBAL[tup_downside[0][0]-1][1]
			   	
	###############connectivity array for FORTRAN#######################
	
	global_node_to_slave_node=[] #connection global slave node label to
								 #slave number
	
	for n in range(1,nNodes+1):
		if n in slaves_dict:
			for d in range(dimens):
				global_node_to_slave_node.append(slaves_dict[n]*dimens-dimens+1+d)
		else:
			for d in range(dimens):
				global_node_to_slave_node.append(0)
				
	#############create Equations in Abaqus CAE (for postprocessing############
	######################simulations directly in Abaqus)######################
	
	if gen_constraints_CAE:
	
		a=mdb.models[Modelname].rootAssembly
		p=mdb.models[Modelname].parts[Partname]
		
		#create Instance in the assembly
		instance_name='Part-1-1'
		create_instance=a.Instance(name=instance_name,part=p,dependent=ON)
		
		#create two reference points:
		#ref. point x: strain (resp. displ. gradient) E11,E12
		#ref. point y: strain (resp. displ. gradient) E12,E22
		refpoint_feature_x=a.ReferencePoint(point=(1,0,0)).id
		change_name=a.features.changeKey(fromName='RP-1',toName='ref_point_x_E11_E12')
		refpoint_feature_y=a.ReferencePoint(point=(0,1,0)).id
		change_name=a.features.changeKey(fromName='RP-1',toName='ref_point_y_E21_E22')
		#create sets containing the referencepoints
		r=a.referencePoints
		create_set=a.Set(name='set_ref_point_x_E11_E12',referencePoints=(r[refpoint_feature_x],))
		create_set=a.Set(name='set_ref_point_y_E21_E22',referencePoints=(r[refpoint_feature_y],))
		#loop through all pairs of slave/masters
		for i in range(len(slaves)):
			#create sets with slave/master nodes
			create_set=a.SetFromNodeLabels(name='Set-m-'+str(i),nodeLabels=(('Part-1-1',(slave_to_master_dict[slaves[i]],)),))
			create_set=a.SetFromNodeLabels(name='Set-s-'+str(i),nodeLabels=(('Part-1-1',(slaves[i],)),))
			#get nodal coordinates
			coords_m=a.sets['Set-m-'+str(i)].nodes[0].coordinates
			coords_s=a.sets['Set-s-'+str(i)].nodes[0].coordinates
			#assigne equations
			create_set=mdb.models[Modelname].Equation(name='Eqn-'+str(i)+'-1', terms=((-1.0, 'Set-s-'+str(i), 1), (1.0, 'Set-m-'+str(i) , 1),(coords_s[0]-coords_m[0],'set_ref_point_x_E11_E12', 1),(coords_s[1]-coords_m[1],'set_ref_point_x_E11_E12', 2)))
			create_set=mdb.models[Modelname].Equation(name='Eqn-'+str(i)+'-2', terms=((-1.0, 'Set-s-'+str(i), 2), (1.0, 'Set-m-'+str(i) , 2),(coords_s[0]-coords_m[0],'set_ref_point_y_E21_E22', 1),(coords_s[1]-coords_m[1],'set_ref_point_y_E21_E22', 2)))
			
		#create step
		mdb.models[Modelname].StaticStep(name='Step-1', previous='Initial')
		#create amplitudes, later to be filled with the real strain/displacement gradient values encountered by the RVE
		amplitude=mdb.models[Modelname].TabularAmplitude(name='E11',data=[(0,0),(0,0),])
		amplitude=mdb.models[Modelname].TabularAmplitude(name='E12',data=[(0,0),(0,0),])
		amplitude=mdb.models[Modelname].TabularAmplitude(name='E21',data=[(0,0),(0,0),])
		amplitude=mdb.models[Modelname].TabularAmplitude(name='E22',data=[(0,0),(0,0),])
		#create boundary conditions (strains/displacement gradient)
		mdb.models[Modelname].DisplacementBC(name='E11', createStepName='Step-1',region=regionToolset.Region(referencePoints=(r[refpoint_feature_x],)), u1=1.0, u2=UNSET, ur3=UNSET, amplitude='E11', fixed=OFF, distributionType=UNIFORM, fieldName='', localCsys=None)
		mdb.models[Modelname].DisplacementBC(name='E12', createStepName='Step-1',region=regionToolset.Region(referencePoints=(r[refpoint_feature_x],)), u1=UNSET, u2=1.0, ur3=UNSET, amplitude='E12', fixed=OFF, distributionType=UNIFORM, fieldName='', localCsys=None)
		mdb.models[Modelname].DisplacementBC(name='E21', createStepName='Step-1',region=regionToolset.Region(referencePoints=(r[refpoint_feature_y],)), u1=1.0, u2=UNSET, ur3=UNSET, amplitude='E21', fixed=OFF, distributionType=UNIFORM, fieldName='', localCsys=None)
		mdb.models[Modelname].DisplacementBC(name='E22', createStepName='Step-1',region=regionToolset.Region(referencePoints=(r[refpoint_feature_y],)), u1=UNSET, u2=1.0, ur3=UNSET, amplitude='E22', fixed=OFF, distributionType=UNIFORM, fieldName='', localCsys=None)
		#create fixed boundary condition
		mdb.models[Modelname].DisplacementBC(name='fixed_point', createStepName='Step-1',region=regionToolset.Region(nodes=a.instances[instance_name].nodes.sequenceFromLabels((boundarypoint,))), u1=0.0, u2=0.0, ur3=UNSET, amplitude=UNSET, fixed=OFF, distributionType=UNIFORM, fieldName='', localCsys=None)
		
	#################return lists#######################################
	
	return (slaves, slaves_dict, slave_to_master_dict, boundarypoint,
			nslaves,delta_x,global_node_to_slave_node)

def get_slave_to_master_connnection3D(Modelname,dimens,COORDSGLOBAL,nNodes,PBC):
	
	#In the current implementation no stand alone slave-master connection
	#is implemented. The connection must be defined in Abaqus so called
	#equations. When the slave-master connection directly lies in equations
	#its assumed, that the master is the first entry in 'terms'
	#then comes the slaves and the length of 'terms' is equal to 3; The
	#equations have the name 'Eqn-1'-'Eqn-n'. When the Micromechanics Plugin
	#was being used, then there exists sorted lists RVE_AG.._AGSORT where
	#the ..POS.. list contains the master and the ..NEG.. list contains the
	#slaves
	
	###############loop over all equations##############################
	
	slaves=[]
	slave_to_master_dict={}
	
	if PBC=='FoamGUI':
	
		for i in range(1,len(mdb.models[Modelname].constraints)+1):
			
			if len(mdb.models[Modelname].constraints['Eqn-'+str(i)].terms)==3:
				#equation sets of slave and master
				eq_set_m=mdb.models[Modelname].constraints['Eqn-'+str(i)].terms[0][1]
				eq_set_s=mdb.models[Modelname].constraints['Eqn-'+str(i)].terms[1][1]
				
				m=mdb.models[Modelname].rootAssembly.sets[eq_set_m].nodes[0].label
				s=mdb.models[Modelname].rootAssembly.sets[eq_set_s].nodes[0].label
				
				slave_to_master_dict[s]=m
				
	elif PBC=='Micromechanics Plugin':
		
		#get slave to master connection for x, y and z face
		
		direct_path=os.path.abspath(".")
		f=open(direct_path+'/'+Modelname+'.inp',"r") #open inputfile
		f_lines=f.readlines()
		
		
		master_to_slave_x={}
		master_x=[]
		slave_x=[]
		master_to_slave_y={}
		master_y=[]
		slave_y=[]
		master_to_slave_z={}
		master_z=[]
		slave_z=[]
		
		indicator='none'
		for f_line in f_lines:
			if f_line[0]=='*':
				if f_line[0:34]=='*NSET, NSET=RVE_AG_NegX_RVE_AGSORT':
					indicator='NEGX'
				elif f_line[0:34]=='*NSET, NSET=RVE_AG_PosX_RVE_AGSORT':
					indicator='POSX'
				elif f_line[0:34]=='*NSET, NSET=RVE_AG_NegY_RVE_AGSORT':
					indicator='NEGY'
				elif f_line[0:34]=='*NSET, NSET=RVE_AG_PosY_RVE_AGSORT':
					indicator='POSY'
				elif f_line[0:34]=='*NSET, NSET=RVE_AG_NegZ_RVE_AGSORT':
					indicator='NEGZ'
				elif f_line[0:34]=='*NSET, NSET=RVE_AG_PosZ_RVE_AGSORT':
					indicator='POSZ'
				else:
					indicator='none'
			elif indicator!='none':
				if indicator=='NEGX':
					master_x.extend([int(number) for number in f_line.split(',')])
				elif indicator=='POSX':
					slave_x.extend([int(number) for number in f_line.split(',')])
				elif indicator=='NEGY':
					master_y.extend([int(number) for number in f_line.split(',')])
				elif indicator=='POSY':
					slave_y.extend([int(number) for number in f_line.split(',')])
				elif indicator=='NEGZ':
					master_z.extend([int(number) for number in f_line.split(',')])
				elif indicator=='POSZ':
					slave_z.extend([int(number) for number in f_line.split(',')])
					
		f.close() #close Inputfile
		
		for i in range(len(master_x)):
			master_to_slave_x[master_x[i]]=slave_x[i]
		for i in range(len(master_y)):
			master_to_slave_y[master_y[i]]=slave_y[i]
		for i in range(len(master_z)):
			master_to_slave_z[master_z[i]]=slave_z[i]
				
		for m in master_z:
			#edges
			if (m in master_to_slave_y) and (not m in master_to_slave_x):
				slave_to_master_dict[master_to_slave_y[master_to_slave_z[m]]]=m
				slave_to_master_dict[master_to_slave_y[m]]=master_to_slave_z[m]
			elif (m in master_to_slave_x) and (not m in master_to_slave_y):
				slave_to_master_dict[master_to_slave_x[master_to_slave_z[m]]]=m
				slave_to_master_dict[master_to_slave_x[m]]=master_to_slave_z[m]
			#corners
			elif (m in master_to_slave_x) and (m in master_to_slave_y):
				#corner1
				slave_to_master_dict[master_to_slave_x[master_to_slave_y[master_to_slave_z[m]]]]=m
				#corner2
				slave_to_master_dict[master_to_slave_x[master_to_slave_y[m]]]=master_to_slave_z[m]
				#corner3
				slave_to_master_dict[master_to_slave_x[master_to_slave_z[m]]]=master_to_slave_y[m]
				#corner4
				slave_to_master_dict[master_to_slave_x[m]]=master_to_slave_y[master_to_slave_z[m]]
			#inner nodes
			else:
				slave_to_master_dict[master_to_slave_z[m]]=m
		for m in master_y:
			#edges
			if m in master_to_slave_x:
				if not master_to_slave_x[master_to_slave_y[m]] in slave_to_master_dict:
					slave_to_master_dict[master_to_slave_x[master_to_slave_y[m]]]=m
					slave_to_master_dict[master_to_slave_x[m]]=master_to_slave_y[m]
			#inner nodes
			elif not master_to_slave_y[m] in slave_to_master_dict:
				slave_to_master_dict[master_to_slave_y[m]]=m
		for m in master_x:
			if not master_to_slave_x[m] in slave_to_master_dict:
				#inner nodes
				slave_to_master_dict[master_to_slave_x[m]]=m
	
	for i in range(1,nNodes+1):
		if i in slave_to_master_dict:
			slaves.append(i)
					
	slaves_dict={} #connection slave node label <-> global node label
	for i in range(len(slaves)):
		slaves_dict[slaves[i]]=i+1
		
	nslaves=len(slave_to_master_dict) #number of slave nodes
	
	boundarypoint=slave_to_master_dict[slaves[0]] #set an arbitrary master with
												  #fixed displacement boundary
												  #condition
	
	###############connectivity array for FORTRAN#######################
	
	global_node_to_slave_node=[] #connection global slave node label to
								 #slave number
	
	for n in range(1,nNodes+1):
		if n in slaves_dict:
			for d in range(dimens):
				global_node_to_slave_node.append(slaves_dict[n]*dimens-dimens+1+d)
		else:
			for d in range(dimens):
				global_node_to_slave_node.append(0)
			   
	#########################get RVE dimensions#########################
	
	delta_x=[0.0]*dimens
	x=[]
	y=[]
	z=[]
		
	for s in slaves:
		#subtract the coordinates of slave and masters and assigne abs to x,y,z
		x.append(abs(COORDSGLOBAL[slave_to_master_dict[s]-1][0]-COORDSGLOBAL[s-1][0]))
		y.append(abs(COORDSGLOBAL[slave_to_master_dict[s]-1][1]-COORDSGLOBAL[s-1][1]))
		z.append(abs(COORDSGLOBAL[slave_to_master_dict[s]-1][2]-COORDSGLOBAL[s-1][2]))

	delta_x[0]=max(x)
	delta_x[1]=max(y)
	delta_x[2]=max(z)
			
	#################return lists#######################################


	return (slaves, slaves_dict, slave_to_master_dict, boundarypoint,
			nslaves,delta_x,global_node_to_slave_node)
			
			

def get_element_to_dof(boundarypoint,nNodes,nElem,nnodes_e,dimens,
	element_to_node,slave_to_master_dict,slaves):

	###################get the reduced dof's############################
	
	node_label_to_red_dof={}	#gives connection: master-> node label
								#to reduced dof; dependend -> node label
								#to connected master reduced dof 
	master_reduced_to_global=[]	#gives connection between master reduced
								#dof entries to master node label
	a=1 #counts from [1 to (dimens*nNodes-dimens*nSlaves-dimens)]
	
	#------------------------------------------------------------------#
	node_label_to_red_dof[boundarypoint]=0	#boundary condition
	
	for n in range(nNodes):	#independend nodes get positive sign
		if (n+1!=boundarypoint and not(n+1 in slave_to_master_dict)):
			node_label_to_red_dof[n+1]=a
			a=a+1
			for d in range(dimens):
				master_reduced_to_global.append((n+1)*dimens-dimens+1+d)
	
	for n in slaves:	#dependend nodes get negative sign
		node_label_to_red_dof[n]=\
			(-1)*node_label_to_red_dof[slave_to_master_dict[n]]
	#------------------------------------------------------------------#
		
	###################get element connections##########################
	
	element_to_global_dof=[] #for each element:element dof to global dof
	element_to_reduced_dof=[]#for each element:dof to reduced system dof
							 #with positiv number==independend
							 #negativ number==dependend
	
	for e in range(nElem):	#loop over all elements
		A=[]
		for i in range(nnodes_e):	#get element_to_global_dof
			for d in range(dimens):
				A.append(element_to_node[e][i]*dimens-dimens+1+d)
		element_to_global_dof.append(A)
		
		A=[]
		for i in range(nnodes_e):	#get element_to_reduced_dof
			for d in range(dimens):
				if node_label_to_red_dof[element_to_node[e][i]]>0:
					A.append(node_label_to_red_dof
							  [element_to_node[e][i]]*dimens-dimens+1+d)
				elif node_label_to_red_dof[element_to_node[e][i]]<0:
					A.append(node_label_to_red_dof
							  [element_to_node[e][i]]*dimens+dimens-1-d)
				else:
					A.append(0)
		element_to_reduced_dof.append(A)
		
	################slave label to master reduced dof###################
		
	global_master_reduced_to_slave=[] #connection slave node label <->
									  #master reduced dof
	for n in slaves:
		for d in range(dimens):
			if node_label_to_red_dof[n]!=0:
				global_master_reduced_to_slave.append(n*dimens-dimens+1+d)
				global_master_reduced_to_slave.append(node_label_to_red_dof[n]
				                                      *(-1)*dimens-dimens+1+d)
		
	#############################return#################################
	
	return (master_reduced_to_global, node_label_to_red_dof,
		    element_to_global_dof, element_to_reduced_dof,
		    global_master_reduced_to_slave)
	
def get_equations(dimens,slave_to_master_dict,COORDSGLOBAL,slaves):	
	#"equations" -> u_slave_i = u_master_i + Grad_U_ij * (X_master_j-X_slave_j)
	
	equations=[] #[slave label, master label, factor 11,
				 #factor 22,factor 12,factor 21] for 2D case
				 #[slave label, master label, factor 11,
				 #factor 22, factor 33, factor 12, factor 21, factor 23 
				 #factor 32, factor 13, factor 31 for 3D case
				 
	#used matrix notation for deformation gradient:
	ordering={3:{(1,1):0,(2,2):1,(3,3):2,(1,2):3,(2,1):4,(2,3):5,
		      (3,2):6,(1,3):7,(3,1):8},
		      2:{(1,1):0,(2,2):1,(1,2):2,(2,1):3}}
		
	for n in slaves: #loop over all slaves
		for d in range(dimens): #loop over all slave dof's
			A=[0.0]*(2+dimens*dimens)
			A[0]=n*dimens-dimens+1+d
			A[1]=slave_to_master_dict[n]*dimens-dimens+1+d
			for j in range(dimens):
				A[2+ordering[dimens][(d+1,j+1)]]=COORDSGLOBAL[n-1][j]-\
						COORDSGLOBAL[slave_to_master_dict[n]-1][j]
			equations.append(A)
			
	return equations
	
def get_stiffnessmatrix_entries(nElem, ndof_e, element_to_reduced_dof,UNSYMM):
	#detect the place in the global reduced stiffnessmatrix in CSR
	#Format where to put the element stiffnessmatrix entries ->
	#K_elem(row, column) -> k_matrix_values(i); UNSYMM==True then take
	#whole matrix, UNSYMM==False, only take the upper triangle
	
	##detect taken (row, column)-Entries of the sparse stiffnessmatrix##
	
	K=[] #NonZero entries in the global stiffnessmatrix (row,column)
	element_value_to_index=[] #element stiffnessmatrix entries
	for e in range(nElem): 	#iteration over all elements
		A=[]
		for i in range(ndof_e):
			for j in range(ndof_e):
				if j>=i or UNSYMM:
					c0=abs(element_to_reduced_dof[e][j])
					r0=abs(element_to_reduced_dof[e][i])
					if UNSYMM:
						c1=c0
						r1=r0
					else:
						if c0<r0:
							c1=r0
							r1=c0
						else:
							c1=c0
							r1=r0
					if c1!=0 and r1!=0:
						K.append((r1,c1))
					A.append((r1,c1))
		element_value_to_index.append(A)
		
	K=list(set(K)) #remove double entries from K
	K.sort(key=lambda tup: (tup[0],tup[1])) #sort K

	############get connection K(row, column)->values(i)################
	
	l=0 #current values entrie number
	ri=0 #current row
	c_and_r_to_value={}	#dictionary: {(row,column):i}
	rowIndex_reduced=[]	#Output for reduced System in CSR Format
	columnIndex_reduced=[] #Output for reduced System in CSR Format
	
	for k in K: #loop over all nonzero stiffnessmatrix entries
		columnIndex_reduced.append(k[1])
		l=l+1
		c_and_r_to_value[(k[0],k[1])]=l
		if k[0]!=ri:
			rowIndex_reduced.append(l)
		ri=k[0]
	rowIndex_reduced.append(l+1)
	
	######get connection elementstiffnessmatrix to k_matrix_values######
											   
	values_to_global_reduced=[]	#gives for each element the connection
								#elementstiffness entry (only upper
								#triangle) to k_matrix_values
	
	for e in range(nElem): #loop over all elements
		A=[]
		l=0
		for i in range(ndof_e):
			for j in range(ndof_e):
				if j>=i or UNSYMM:
					if (element_value_to_index[e][l][0]==0 or
						element_value_to_index[e][l][1]==0):
							A.append(-1)	#"-1" indicates that element
											#stiffness entry not needed
					else:
						A.append(c_and_r_to_value[element_value_to_index[e][l]])
					l=l+1
				else:
					A.append(-1) #lower triangle not needed for symmetric matricies
		values_to_global_reduced.append(A)
		
	########################return######################################
	
	return (rowIndex_reduced, columnIndex_reduced,values_to_global_reduced)
	
def write_data_to_file(jobname,rve_number, dimens,nnodes_e, ndof_e, ndof_n,
		nElem, nNodes, columnIndex_reduced, columnIndex_reduced_UNSYMM,
		rowIndex_reduced, rowIndex_reduced_UNSYMM, nslaves, ndof_macro, delta_x,
		COORDSGLOBAL,master_reduced_to_global, element_to_global_dof,
		element_to_reduced_dof, values_to_global_reduced,
		values_to_global_reduced_UNSYMM, global_master_reduced_to_slave,
		equations, element_to_material,global_node_to_slave_node,JTYPE,
		information):
	
	#create Inputfile for the FE2 program which is read at the beginning 
	#of the FE analysis
	
	
	direct_path=os.path.abspath(".")
	
	fobj_out = open(direct_path+"/"+jobname+".FE"+str(rve_number),"w")
	
	#output version number
	fobj_out.write('**MonolithFE2 Version 0.1\n')
	
	#general warning
	fobj_out.write('**Do not modify this inputfile, otherwise it may not work properly!\n')
	
	fobj_out.write('*Part, "')							#include a RVE description
	if not information.strip():
		raise Exception('RVE description shall not be empty!')
	fobj_out.write(information)
	fobj_out.write('"\n')
		
	fobj_out.write('*Dimension, N=')
	fobj_out.write(str(dimens))				 			#dimens
	fobj_out.write('\n')
		
	fobj_out.write('*Elementnodes, N=')
	fobj_out.write(str(nnodes_e))					 	#nnodes_e
	fobj_out.write('\n')
		
	fobj_out.write('*NodeDOF, N=')
	fobj_out.write(str(ndof_n))							#ndof_n
	fobj_out.write('\n')
		
	fobj_out.write('*NDOF_macro, N=')
	fobj_out.write(str(ndof_macro))						#ndof_macro
	fobj_out.write('\n')
		
	fobj_out.write('*Slavenodes, N=')
	fobj_out.write(str(nslaves))					 	#nslaves
	fobj_out.write('\n')
		
	fobj_out.write('*Elementtype, ')
	fobj_out.write(str(JTYPE))							#JTYPE
	fobj_out.write('\n')
		
	fobj_out.write('*Node, N=')
	fobj_out.write(str(nNodes))							#nNodes
	fobj_out.write('\n')
	for i in COORDSGLOBAL:								#COORDSGLOBAL
		for j in range(dimens):
			fobj_out.write(str(i[j]))
			if not(j==(len(i)-1)):
				fobj_out.write(', ')
			else:
				fobj_out.write('\n')
	
	fobj_out.write('*Element_to_global_DOF, N=')
	fobj_out.write(str(nElem))				 			#nElem
	fobj_out.write('\n')
	for i in element_to_global_dof:						#element_to_global_dof
		for j in range(len(i)):
			fobj_out.write(str(i[j]))
			if not(j==(len(i)-1)):
				fobj_out.write(', ')
			else:
				fobj_out.write('\n')
	
	fobj_out.write('*Element_to_reduced_DOF, N=')
	fobj_out.write(str(nElem))				 			#nElem
	fobj_out.write('\n')
	for i in element_to_reduced_dof:					#element_to_reduced_dof
		for j in range(len(i)):
			fobj_out.write(str(i[j]))
			if not(j==(len(i)-1)):
				fobj_out.write(', ')
			else:
				fobj_out.write('\n')
	
	fobj_out.write('*Equations, N=')
	fobj_out.write(str(len(equations)))					#nEquations
	fobj_out.write('\n')
	for i in equations:									#Equations
		for j in range(len(i)):
			fobj_out.write(str(i[j]))
			if not(j==len(i)-1):
				fobj_out.write(', ')
			else:
				fobj_out.write('\n')
	
	fobj_out.write('*Element_to_Material, N=')
	fobj_out.write(str(nElem))				 			#nElem
	fobj_out.write('\n')
	for i in str(element_to_material):					#element_to_material
		if not(i=='[' or i==']'):
			fobj_out.write(i)
	fobj_out.write('\n')
	
	fobj_out.write('*RVE_Dimensions, N=')
	fobj_out.write(str(dimens))				 			#dimens
	fobj_out.write('\n')
	for i in str(delta_x):								#delta_x
		if not(i=='[' or i==']'):
			fobj_out.write(i)
	fobj_out.write('\n')
	
	fobj_out.write('*ColumnIndex_reduced, N=')
	fobj_out.write(str(len(columnIndex_reduced)))		#nNonZero_reduced
	fobj_out.write('\n')
	for i in str(columnIndex_reduced):					#columnIndex_reduced
		if not(i=='[' or i==']'):
			fobj_out.write(i)
	fobj_out.write('\n')
	
	fobj_out.write('*ColumnIndex_reduced_UNSYMM, N=')
	fobj_out.write(str(len(columnIndex_reduced_UNSYMM)))#nNonZero_reduced_UNSYMM
	fobj_out.write('\n')
	for i in str(columnIndex_reduced_UNSYMM):			#columnIndex_reduced_UNSYMM
		if not(i=='[' or i==']'):
			fobj_out.write(i)
	fobj_out.write('\n')


	fobj_out.write('*RowIndex_reduced, N=')
	fobj_out.write(str(len(rowIndex_reduced)))			#nRows_reduced+1
	fobj_out.write('\n')
	for i in str(rowIndex_reduced):						#rowIndex_reduced
		if not(i=='[' or i==']'):
			fobj_out.write(i)
	fobj_out.write('\n')
	
	fobj_out.write('*RowIndex_reduced_unsymm, N=')
	fobj_out.write(str(len(rowIndex_reduced)))			#nRows_reduced+1
	fobj_out.write('\n')
	for i in str(rowIndex_reduced_UNSYMM):			#rowIndex_reduced_UNSYMM
		if not(i=='[' or i==']'):
			fobj_out.write(i)
	fobj_out.write('\n')
	
	fobj_out.write('*values_to_global_reduced, N=')
	fobj_out.write(str(nElem))				 			#nElem
	fobj_out.write('\n')
	for i in values_to_global_reduced:					#values_to_global_reduced
		for j in range(len(i)):
			fobj_out.write(str(i[j]))
			if not(j==len(i)-1):
				fobj_out.write(', ')
			else:
				fobj_out.write('\n')
	
	fobj_out.write('*values_to_global_reduced_UNSYMM, N=')
	fobj_out.write(str(nElem))				 			#nElem
	fobj_out.write('\n')
	for i in values_to_global_reduced_UNSYMM:					#values_to_global_reduced
		for j in range(len(i)):
			fobj_out.write(str(i[j]))
			if not(j==len(i)-1):
				fobj_out.write(', ')
			else:
				fobj_out.write('\n')
	
	fobj_out.write('*master_reduced_to_global, N=')
	fobj_out.write(str(len(master_reduced_to_global)))			#nRows_reduced
	fobj_out.write('\n')
	for i in str(master_reduced_to_global):				#master_reduced_to_global
		if not(i=='[' or i==']'):
			fobj_out.write(i)
	fobj_out.write('\n')
	
	fobj_out.write('*global_master_reduced_to_slave, N=')
	fobj_out.write(str(len(global_master_reduced_to_slave)/2))#number master-slave
	fobj_out.write('\n')
	j=1
	for i in global_master_reduced_to_slave:#global_master_reduced_to_slave
		fobj_out.write(str(i))
		if j==1:
			fobj_out.write(', ')
			j=2
		else:
			fobj_out.write('\n')
			j=1
	
	fobj_out.write('*global_node_to_slave_node, N=')
	fobj_out.write(str(ndof_n*nNodes))
	fobj_out.write('\n')
	for i in str(global_node_to_slave_node):			#global_node_to_slave_node
		if not(i=='[' or i==']'):
			fobj_out.write(i)
	fobj_out.write('\n')
	
	fobj_out.close()
