#=============================================================================
# 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 odbAccess
from viewerModules import *
from abaqus import session
from scipy.cluster.vq import vq, kmeans
import odbAccess
import regionToolset
import os   #import os functionalities to get current directory path
import numpy as np #import numpy for numerical tools
import sys
import subprocess
import platform
from Hyperreduction import *
from set_FE2_Analysisparameters import set_FE2_Analysisparameters

class meshparameters: #define the same object as in FORTRAN MODULE type_meshparameters
    def __init__(para,dimens,Modelname,Partname): #initialize the object
        
        para.Modelname=Modelname #name of the Abaqus Model
        para.Partname=Partname #name of the Abaqus Part in that Model
        para.dimens=dimens #number of dimensions (2D: 2, 3D:3)
        if para.dimens==2:          #number of macroscopic degrees of freedom (e.g in
            para.ndof_macro=3       #purely mechanical problems is 6 independent entries 
        elif para.dimens==3:        #of the cauchy stress resp. biot stress (finite deformations)
            para.ndof_macro=6       #-> needs to be changed in case of generalized continua
        para.ndof_n=0 #number of dof of a node
        para.JTYPE=0 #element type of UEL lib
        para.NGP_local=0 #number of integration points of the micro elements
        para.equations=[] #equations for enforcing the periodic boundary conditions (constraint)
        para.master_reduced_to_global=[] #connection master reduced dof <-> global master node label
        para.rowIndex_reduced=[]#CSR definition array for reduced stiffnessmatrix
        para.columnIndex_reduced=[] #CSR definition array for reduced stiffnessmatrix
        para.values_to_global_reduced=[] #element stiffnessmatrix to k_matrix_values
        #for unsymmetric stiffness matrix
        para.rowIndex_reduced_UNSYMM=[]
        para.columnIndex_reduced_UNSYMM=[]
        para.values_to_global_reduced_UNSYMM=[]
        para.COORDSGLOBAL=[] #coordinates of the Mesh
        para.element_to_node=[] #connection element to node labels
        para.element_to_global_dof=[] #element connection :element dof <-> global dof label
        para.element_to_reduced_dof=[] #element connection :element dof <-> reduced system dof
        para.element_to_material=[]#conntection element <-> material label
        para.global_node_to_slave_node=[] #connection: node label <-> slave label
        para.RVE_Volume=0 #volume of the RVE (including holes, pores etc.)
        para.PROPS=[] #array containing a set of materialparameters in each row
        para.slaves=[] #contains all slave node labels
        para.slaves_dict={} #connection slave node label <-> global node label
        para.slave_to_master_dict={} #slave<->master connection
        para.boundarypoint=0 #displacements boundary condition zero -> lock rigid body motion
                             #used matrix notation for strain/stretch
        para.ordering={3:{(1,1):(0,1.0),(2,2):(1,1.0),(3,3):(2,1.0),(1,2):(3,0.5),
                          (2,1):(3,0.5),(1,3):(4,0.5),(3,1):(4,0.5),(2,3):(5,0.5),(3,2):(5,0.5),
                          0:(1,1),1:(2,2),2:(3,3),3:(1,2),4:(1,3),5:(2,3)},
                       2:{(1,1):(0,1.0),(2,2):(1,1.0),(1,2):(2,0.5),(2,1):(2,0.5),
                          0:(1,1),1:(2,2),2:(1,2)}}
        para.Node_ABQ_to_Node={} #Abaqus meshes could have nodes not connected to any element -> remove them
                                 #and name all nodes from 1 to n
        
    def get_data_from_abaqus(para):
        
        ###################Elements#############################################
        
        #get information of all Elements
        ELEMENTS=mdb.models[para.Modelname].parts[para.Partname].elements
        
        ############################get elementtype JTYPE###########################
        
        #-> assumption: all elements have the same amount of nodes!
        
        #elementlibrary={(Abaqus Elementtype):JTYPE,NGP_local,ndof_n}
        elementlibrary={'CPE4':(2001,4,2),
                        'CPE8':(2002,9,2),
                        'CPS4':(2003,4,2),
                        'C3D8':(2004,8,3),
                        'C3D20':(2005,27,3),
                        'CPE4R':(2006,1,2),
                        'CPE8R':(2007,4,2),
                        'CPS4R':(2008,1,2),
                        'CPS8':(2009,9,2),
                        'CPS8R':(2010,4,2),
                        'C3D20R':(2011,8,3),
                        'CPE3':(2012,1,2),
                        'C3D4':(2013,1,3),
                        'C3D10':(2014,4,3),
                        'CPS3':(2024,1,2),
                        'CPE6':(2029,3,2),
                        'CPS6':(2030,3,2)}
                        
        JTYPE_List=list(set(ele.type.name for ele in ELEMENTS))
        
        if len(JTYPE_List)>1:
            WARNUNG=getWarningReply('All elements in the mesh need to have the same element type. Please remesh the RVE!',(CANCEL,))
            raise Exception ('Execution stopped because of some problem')
                        
        if JTYPE_List[0] in elementlibrary:
            #get element type
            para.JTYPE=elementlibrary[JTYPE_List[0]][0]
            para.NGP_local=elementlibrary[JTYPE_List[0]][1]
            para.ndof_n=elementlibrary[JTYPE_List[0]][2]
        else:
            #display error
            WARNUNG=getWarningReply('This elementtype is currently not supported by the UELlib.\n'+
                            'Supported Elements: CPE4, CPE8, CPS4, C3D8, C3D20, CPE4R,\n'+
                            'CPE8R, CPS4R, CPS8, CPS8R, C3D20R, CPE3, C3D4, C3D10, CPS3,\n'+
                            'CPE6, CPS6',(CANCEL,))
            raise Exception ('Execution stopped because of some problem')
        
        ###################get element-to node connection###########################
        
        #assure that only the nodes are taken, which appear in the element to node
        #connection and that they are numbered from 1...nNodes
        
        labelmax=0
        Node_ABQ_to_Node={}
        for e in range(len(ELEMENTS)): #at first get all nodes which appear in the elements
            for i in range(len(ELEMENTS[0].connectivity)):
                #be aware: python and the node label counts from 1, but only in the MeshElement it counts from 0!
                Node_ABQ_to_Node[ELEMENTS[e].connectivity[i]]=0
                if ELEMENTS[e].connectivity[i]>labelmax:
                    labelmax=ELEMENTS[e].connectivity[i]
        
        label=1
        for i in range(labelmax+1): #give all ABQ node labels a new label starting from 1
            if i in Node_ABQ_to_Node:
                Node_ABQ_to_Node[i]=label
                label+=1
        
        for e in range(len(ELEMENTS)):  #get element_to_node connectivity from ELEMENTES
            para.element_to_node.append([Node_ABQ_to_Node[ELEMENTS[e].connectivity[i]] for i in range(len(ELEMENTS[e].connectivity))])
        
        ########################## material properties ################################
        
        ASSIGNMENT=mdb.models[para.Modelname].parts[para.Partname].sectionAssignments
        
        para.element_to_material=[0]*len(para.element_to_node) #Element_label ->assignment-> Material_label
        for a in range(len(ASSIGNMENT)):
            A=ASSIGNMENT[a].getSet().elements
            for i in A:
                para.element_to_material[i.label-1]=a+1
            S=mdb.models[para.Modelname].sections[ASSIGNMENT[a].sectionName]
            try:
                para.PROPS.append(mdb.models[para.Modelname].materials[S.material].userMaterial.mechanicalConstants) #Material_label ->assignment-> [list of properties]
            except:
                WARNUNG=getWarningReply('All elements must be assigned to a region having a material with a user material and not empty constants',(CANCEL,))
                raise Exception ('Execution stopped because of some problem')
            
        for e in para.element_to_material:
            if e==0:
                WARNUNG=getWarningReply('At least one element does not have a material definition!',(CANCEL,))
                raise Exception ('Execution stopped because of some problem')
        
        ###################Coordinates##########################################
        
        NODES=mdb.models[para.Modelname].parts[para.Partname].nodes #get node infos
        
        for n in range(len(NODES)): #loop through nodes to get coordinates
            if n in Node_ABQ_to_Node:
                para.COORDSGLOBAL.append([NODES[n].coordinates[i] for i in range(para.dimens)])
                para.Node_ABQ_to_Node[NODES[n].label]=Node_ABQ_to_Node[n]
                        
    def get_slave_to_master_connnection2D(para,upside_abq,downside_abq,leftside_abq, 
                                        rightside_abq):
        
        ###################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(para.Node_ABQ_to_Node[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(para.Node_ABQ_to_Node[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(para.Node_ABQ_to_Node[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(para.Node_ABQ_to_Node[rightside_abq[j].label])
        rightside=list(set(rightside)) #remove potential double entries
        
        if len(upside)!=len(downside) or len(leftside)!=len(rightside):
            WARNUNG=getWarningReply('Meshes must be congruent!',(CANCEL,))
            raise Exception ('Execution stopped because of some problem')
        
        ##############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],para.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],para.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],para.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],para.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)
            
            para.slaves=[rightup,rightdown,leftdown]
            
            para.slave_to_master_dict={} #slave<->master connection
            para.slave_to_master_dict[rightup]=leftup   #corner nodes
            para.slave_to_master_dict[rightdown]=leftup
            para.slave_to_master_dict[leftdown]=leftup
            
            para.boundarypoint=leftup
        
        else:
            
            para.boundarypoint=tup_leftside[len(tup_leftside)-1][0]
    
    
        ######################sort slaves in lists##########################
        
        #slaves contains all slave node labels
        
        for t_r in tup_rightside:
            para.slaves.append(t_r[0])
        
        for t_u in tup_upside:
            para.slaves.append(t_u[0])
        
        para.slaves.sort()  #sort the slaves list
        
        for i in range(len(para.slaves)):
            para.slaves_dict[para.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
            para.slave_to_master_dict[tup_rightside[i][0]]=tup_leftside[i][0]
            
        for i in range(len(tup_downside)):  #edge nodes
            para.slave_to_master_dict[tup_upside[i][0]]=tup_downside[i][0]
            
        ################additional information##############################
            
        delta_x=[0]*para.dimens #length and width of the RVE
        
        delta_x[0]=para.COORDSGLOBAL[tup_rightside[0][0]-1][0]-\
                para.COORDSGLOBAL[tup_leftside[0][0]-1][0]
        delta_x[1]=para.COORDSGLOBAL[tup_upside[0][0]-1][1]-\
                para.COORDSGLOBAL[tup_downside[0][0]-1][1]
        
        para.RVE_Volume=delta_x[0]*delta_x[1]
        
        
    def get_slave_to_master_connnection3D(para,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##############################
        
        if PBC=='3D - Foam GUI':
        
            for i in range(1,len(mdb.models[para.Modelname].constraints)+1):
                
                if len(mdb.models[para.Modelname].constraints['Eqn-'+str(i)].terms)==3:
                    #equation sets of slave and master
                    eq_set_m=mdb.models[para.Modelname].constraints['Eqn-'+str(i)].terms[0][1]
                    eq_set_s=mdb.models[para.Modelname].constraints['Eqn-'+str(i)].terms[1][1]
                    
                    m=para.Node_ABQ_to_Node[mdb.models[para.Modelname].rootAssembly.sets[eq_set_m].nodes[0].label]
                    s=para.Node_ABQ_to_Node[mdb.models[para.Modelname].rootAssembly.sets[eq_set_s].nodes[0].label]
                    
                    para.slave_to_master_dict[s]=m
        
        
        elif PBC=='3D - Micromechanics Plugin':
            
            #get slave to master connection for x, y and z face
            
            direct_path=os.path.abspath(".")
            f=open(direct_path+'/'+para.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]
            
            masters=[]
            
            for m in master_z: #first detect corners
                if (m in master_to_slave_x) and (m in master_to_slave_y):
                    para.slave_to_master_dict[master_to_slave_x[m]]=m
                    para.slave_to_master_dict[master_to_slave_y[m]]=m
                    para.slave_to_master_dict[master_to_slave_z[m]]=m
                    para.slave_to_master_dict[master_to_slave_x[master_to_slave_y[m]]]=m
                    para.slave_to_master_dict[master_to_slave_x[master_to_slave_z[m]]]=m
                    para.slave_to_master_dict[master_to_slave_y[master_to_slave_z[m]]]=m
                    para.slave_to_master_dict[master_to_slave_x[master_to_slave_y[master_to_slave_z[m]]]]=m
                    masters.append(m)
            for m in master_z: #next detect edges
                if (m in master_to_slave_y) and (not (m in master_to_slave_x)) and (not (m in masters)) and (not (m in para.slave_to_master_dict)):
                    para.slave_to_master_dict[master_to_slave_z[m]]=m
                    para.slave_to_master_dict[master_to_slave_y[master_to_slave_z[m]]]=m
                    para.slave_to_master_dict[master_to_slave_y[m]]=m
                    masters.append(m)
                elif (m in master_to_slave_x) and (not (m in master_to_slave_y)) and (not (m in masters)) and (not (m in para.slave_to_master_dict)):
                    para.slave_to_master_dict[master_to_slave_z[m]]=m
                    para.slave_to_master_dict[master_to_slave_x[master_to_slave_z[m]]]=m
                    para.slave_to_master_dict[master_to_slave_x[m]]=m
                    masters.append(m)
            for m in master_x: #detect edges
                if (m in master_to_slave_y) and (not (m in master_to_slave_z))and (not (m in masters)) and (not (m in para.slave_to_master_dict)):
                    para.slave_to_master_dict[master_to_slave_x[m]]=m
                    para.slave_to_master_dict[master_to_slave_y[m]]=m
                    para.slave_to_master_dict[master_to_slave_x[master_to_slave_y[m]]]=m
                    masters.append(m)
            for m in master_z: #finally get the inner nodes
                if (not (m in master_to_slave_x)) and (not (m in master_to_slave_y)) and (not (m in masters)) and (not (m in para.slave_to_master_dict)):
                    para.slave_to_master_dict[master_to_slave_z[m]]=m
                    masters.append(m)
            for m in master_y: #get the inner nodes
                if (not (m in master_to_slave_x)) and (not (m in master_to_slave_z)) and (not (m in masters)) and (not (m in para.slave_to_master_dict)):
                    para.slave_to_master_dict[master_to_slave_y[m]]=m
                    masters.append(m)
            for m in master_x: #get the inner nodes
                if (not (m in master_to_slave_y)) and (not (m in master_to_slave_z)) and (not (m in masters)) and (not (m in para.slave_to_master_dict)):
                    para.slave_to_master_dict[master_to_slave_x[m]]=m
                    masters.append(m)
        
        for i in range(1,len(para.COORDSGLOBAL)+1):
            if i in para.slave_to_master_dict:
                para.slaves.append(i)
                        
        for i in range(len(para.slaves)):
            para.slaves_dict[para.slaves[i]]=i+1
        
        #set an arbitrary point to prevent rigid body motion
        para.boundarypoint=para.slave_to_master_dict[para.slaves[0]]
        
        x=[]
        y=[]
        z=[]
        
        #########################get RVE dimensions#########################
        
        for s in para.slaves:
            #subtract the coordinates of slave and masters and assigne abs to x,y,z
            x.append(abs(para.COORDSGLOBAL[para.slave_to_master_dict[s]-1][0]-para.COORDSGLOBAL[s-1][0]))
            y.append(abs(para.COORDSGLOBAL[para.slave_to_master_dict[s]-1][1]-para.COORDSGLOBAL[s-1][1]))
            z.append(abs(para.COORDSGLOBAL[para.slave_to_master_dict[s]-1][2]-para.COORDSGLOBAL[s-1][2]))
        
        para.RVE_Volume=max(x)*max(y)*max(z)
        
    def get_element_to_dof(para):
    
        ###################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 
        a=1 #counts from [1 to (ndof_n*nNodes-ndof_n*nSlaves-ndof_n)]
        
        #------------------------------------------------------------------#
        node_label_to_red_dof[para.boundarypoint]=0 #boundary condition
        
        for n in range(1,len(para.COORDSGLOBAL)+1): #independend nodes get positive sign
            if (n!=para.boundarypoint and not(n in para.slave_to_master_dict)):
                node_label_to_red_dof[n]=a
                a=a+1
                for d in range(1,para.ndof_n+1):
                    para.master_reduced_to_global.append((n-1)*para.ndof_n+d)
                
        for n in para.slaves:   #dependend nodes get negative sign
            node_label_to_red_dof[n]=(-1)*node_label_to_red_dof[para.slave_to_master_dict[n]]
        #------------------------------------------------------------------#
        
        ###################get element connections##########################
        
        for e in range(len(para.element_to_node)):  #loop over all elements
            A=[]; B=[]
            for i in range(len(para.element_to_node[e])):
                ###########get element_to_global_dof######################
                for d in range(1,para.ndof_n+1):
                    A.append((para.element_to_node[e][i]-1)*para.ndof_n+d)
                ##########get element_to_reduced_dof######################
                    if node_label_to_red_dof[para.element_to_node[e][i]]>0:
                        B.append((node_label_to_red_dof[para.element_to_node[e][i]]-1)*para.ndof_n+d)
                    elif node_label_to_red_dof[para.element_to_node[e][i]]<0:
                        B.append((node_label_to_red_dof[para.element_to_node[e][i]]+1)*para.ndof_n-d)
                    else:
                        B.append(0)
            para.element_to_global_dof.append(A)
            para.element_to_reduced_dof.append(B)
        
        ###############connection global dof to slave dof####################
        
        for n in range(1,len(para.COORDSGLOBAL)+1):
            if n in para.slaves_dict:
                for d in range(1,para.dimens+1):
                    para.global_node_to_slave_node.append((para.slaves_dict[n]-1)*para.dimens+d)
            else:
                para.global_node_to_slave_node.extend([0]*para.dimens)
        
    
    def get_stiffnessmatrix_entries(para):
        #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##
        
        for UNSYMM in [True,False]:
            K=[] #NonZero entries in the global stiffnessmatrix (row,column)
            element_value_to_index=[] #element stiffnessmatrix entries
            for e in range(len(para.element_to_reduced_dof)):   #iteration over all elements
                A=[]
                for i in range(len(para.element_to_reduced_dof[e])):
                    for j in range(len(para.element_to_reduced_dof[e])):
                        if j>=i or UNSYMM:
                            c0=abs(para.element_to_reduced_dof[e][j])
                            r0=abs(para.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(len(para.element_to_reduced_dof)): #loop over all elements
                A=[]
                l=0
                for i in range(len(para.element_to_reduced_dof[e])):
                    for j in range(len(para.element_to_reduced_dof[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)
            
            if UNSYMM==True:
                para.rowIndex_reduced_UNSYMM=rowIndex_reduced
                para.columnIndex_reduced_UNSYMM=columnIndex_reduced
                para.values_to_global_reduced_UNSYMM=values_to_global_reduced
            else:
                para.rowIndex_reduced=rowIndex_reduced
                para.columnIndex_reduced=columnIndex_reduced
                para.values_to_global_reduced=values_to_global_reduced
            
    def get_equations(para):
        #"equations" -> u_slave_i = u_master_i + Grad_U_ij * (X_master_j-X_slave_j)
        
        #equations contain:  [slave label, master label, factor 11,
                            #factor 22,factor 12] for 2D case
                            #[slave label, master label, factor 11,
                            #factor 22, factor 33, factor 12, factor 13 
                            #factor 23] for 3D case
            
        for n in para.slaves: #loop over all slaves
            for i in range(1,para.dimens+1): #loop over all slave dof's
                A=[0.0]*(2+para.ndof_macro)
                A[0]=(n-1)*para.dimens+i
                A[1]=(para.slave_to_master_dict[n]-1)*para.dimens+i
                for j in range(1,para.dimens+1):
                    A[2+para.ordering[para.dimens][(i,j)][0]]=(para.COORDSGLOBAL[n-1][j-1]-\
                    para.COORDSGLOBAL[para.slave_to_master_dict[n]-1][j-1])*para.ordering[para.dimens][(i,j)][1]
                para.equations.append(A)
    
    def create_equations_CAE(para,PBC):
        #create Equations in Abaqus CAE (for postprocessing simulations directly in Abaqus)
        
        a=mdb.models[para.Modelname].rootAssembly
        p=mdb.models[para.Modelname].parts[para.Partname]
        
        #create Instance in the assembly
        if len(a.allInstances.values())<1:
            instance_name=para.Partname+'-1'
            create_instance=a.Instance(name=instance_name,part=p,dependent=ON)
        elif len(a.allInstances.values())==1:
            instance_name=a.allInstances.values()[0].name
        else:
            WARNUNG=getWarningReply('Only one instance allowed.',(CANCEL,))
            raise Exception ('Execution stopped because of some problem')
        
        #create three reference points:  strain (resp. displ. gradient) E11,E22,E12 etc.
        refpoint_features=[]
        for i in range(para.ndof_macro):
            refpoint_features.append(a.ReferencePoint(point=(0,0,0)).id)
            change_name=a.features.changeKey(fromName='RP-1',toName='ref_point_E'+str(para.ordering[para.dimens][i][0])+str(para.ordering[para.dimens][i][1]))
        #create sets containing the referencepoints
        r=a.referencePoints
        for i in range(para.ndof_macro):
            a.Set(name='set_ref_point_E'+str(para.ordering[para.dimens][i][0])+str(para.ordering[para.dimens][i][1]),referencePoints=(r[refpoint_features[i]],))
        #loop through all pairs of slave/masters
        for i in range(len(para.slaves)):
            #create sets with slave/master nodes
            a.SetFromNodeLabels(name='Set-m-'+str(i),nodeLabels=((instance_name,(para.slave_to_master_dict[para.slaves[i]],)),))
            a.SetFromNodeLabels(name='Set-s-'+str(i),nodeLabels=((instance_name,(para.slaves[i],)),))
            #assigne equations
            for j in range(para.dimens):
                Terms=[(-1.0, 'Set-s-'+str(i), j+1),(1.0, 'Set-m-'+str(i) , j+1)]
                for k in range(para.ndof_macro):
                    Terms.append((para.equations[i*para.dimens+j][2+k],'set_ref_point_E'+str(para.ordering[para.dimens][k][0])+str(para.ordering[para.dimens][k][1]),1))
                mdb.models[para.Modelname].Equation(name='Eqn-'+str(i)+'-'+str(j), terms=tuple(Terms))
        #create step
        stepname='Step-1'
        mdb.models[para.Modelname].StaticStep(name=stepname, previous='Initial')
        for k in range(para.ndof_macro):
            strain_name='E'+str(para.ordering[para.dimens][k][0])+str(para.ordering[para.dimens][k][1])
            #create amplitudes, later to be filled with the real strain/displacement gradient values encountered by the RVE
            amplitude=mdb.models[para.Modelname].TabularAmplitude(name=strain_name,data=[(0,0),(0,0),])
            #create boundary conditions (strains/displacement gradient)
            mdb.models[para.Modelname].DisplacementBC(name=strain_name, createStepName=stepname,region=regionToolset.Region(referencePoints=(r[refpoint_features[k]],)), u1=1.0,amplitude=strain_name, fixed=OFF, distributionType=UNIFORM, fieldName='', localCsys=None)
        #create fixed boundary condition
        n=1
        while True:
            if n in para.Node_ABQ_to_Node:
                if para.Node_ABQ_to_Node[n]==para.boundarypoint:
                    boundarypoint=n #get the actual label that Abaqus uses (not the MonolithFE2 label, that could be different)
                    break
            n=n+1
        mdb.models[para.Modelname].DisplacementBC(name='fixed_point', createStepName=stepname,region=regionToolset.Region(nodes=a.instances[instance_name].nodes.sequenceFromLabels((boundarypoint,))), u1=0.0, u2=0.0, u3=0.0, amplitude=UNSET, fixed=OFF, distributionType=UNIFORM, fieldName='', localCsys=None)
    
    def write_data_to_file(para,jobname,rve_number,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 2.0\n')
        
        #general warning
        fobj_out.write('**Do not modify this inputfile, otherwise it may not work properly!\n')
        
        fobj_out.write('*Part, "')
        if not information.strip():
            WARNUNG=getWarningReply('RVE description shall not be empty.',(CANCEL,))
            raise Exception ('Execution stopped because of some problem')
        fobj_out.write(information)
        fobj_out.write('"\n')
            
        fobj_out.write('*Dimension, N=')
        fobj_out.write(str(para.dimens))
        fobj_out.write('\n')
        
        fobj_out.write('*RVE_Volume, V=')
        fobj_out.write(str(para.RVE_Volume))
        fobj_out.write('\n')
        
        fobj_out.write('*Elementtype, ')
        fobj_out.write(str(para.JTYPE))
        fobj_out.write('\n')
        
        fobj_out.write('*NGP_local, N=')
        fobj_out.write(str(para.NGP_local))
        fobj_out.write('\n')
        
        fobj_out.write('*Elementnodes, N=')
        fobj_out.write(str(len(para.element_to_node[0])))
        fobj_out.write('\n')
            
        fobj_out.write('*NodeDOF, N=')
        fobj_out.write(str(para.ndof_n))
        fobj_out.write('\n')
            
        fobj_out.write('*NDOF_macro, N=')
        fobj_out.write(str(para.ndof_macro))
        fobj_out.write('\n')
        
        fobj_out.write('*Node, N=')
        fobj_out.write(str(len(para.COORDSGLOBAL)))
        fobj_out.write('\n')
        for i in para.COORDSGLOBAL:
            for j in range(para.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(len(para.element_to_global_dof)))
        fobj_out.write('\n')
        for i in para.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(len(para.element_to_reduced_dof)))
        fobj_out.write('\n')
        for i in para.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('*Element_to_Material, N=')
        fobj_out.write(str(len(para.element_to_material)))
        fobj_out.write('\n')
        for i in str(para.element_to_material):
            if not(i=='[' or i==']'):
                fobj_out.write(i)
        fobj_out.write('\n')
        
        fobj_out.write('*User_Material, N=')
        fobj_out.write(str(len(para.PROPS)))
        fobj_out.write('\n')
        for j in range(len(para.PROPS)):
            fobj_out.write('*Material_'+str(j+1)+', N=')
            fobj_out.write(str(len(para.PROPS[j])))
            fobj_out.write('\n')
            for i in str(para.PROPS[j]):
                if not(i=='[' or i==']'):
                    fobj_out.write(i)
            fobj_out.write('\n')
    
        fobj_out.write('*Equations, N=')
        fobj_out.write(str(len(para.equations)))
        fobj_out.write('\n')
        for i in para.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('*ColumnIndex_reduced, N=')
        fobj_out.write(str(len(para.columnIndex_reduced)))
        fobj_out.write('\n')
        for i in str(para.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(para.columnIndex_reduced_UNSYMM)))
        fobj_out.write('\n')
        for i in str(para.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(para.rowIndex_reduced)))
        fobj_out.write('\n')
        for i in str(para.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(para.rowIndex_reduced_UNSYMM)))
        fobj_out.write('\n')
        for i in str(para.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(len(para.values_to_global_reduced)))
        fobj_out.write('\n')
        for i in para.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(len(para.values_to_global_reduced_UNSYMM)))
        fobj_out.write('\n')
        for i in para.values_to_global_reduced_UNSYMM:
            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(para.master_reduced_to_global)))
        fobj_out.write('\n')
        for i in str(para.master_reduced_to_global):
            if not(i=='[' or i==']'):
                fobj_out.write(i)
        fobj_out.write('\n')
        
        fobj_out.write('*global_node_to_slave_node, N=')
        fobj_out.write(str(len(para.global_node_to_slave_node)))
        fobj_out.write('\n')
        for i in str(para.global_node_to_slave_node):
            if not(i=='[' or i==']'):
                fobj_out.write(i)
        fobj_out.write('\n')
        
        fobj_out.close()
    
def generate_postprocessing_data(macro_odbName,element,integration_point,source,para):
    
    data_error=0
    
    for i in range(para.ndof_macro):
        strain_name='E'+str(para.ordering[para.dimens][i][0])+str(para.ordering[para.dimens][i][1])
        #extract the strains from the odb
        if source=='SDV':
            try:
                data=session.xyDataListFromField(odb=session.odbs[macro_odbName], outputPosition=INTEGRATION_POINT, 
                variable=(('SDV'+str(i+1), INTEGRATION_POINT), ), elementLabels=((element.instanceName,element.label),))
            except:
                WARNUNG=getWarningReply('No SDV odb data created by MonolithFE2 was found!',(CANCEL,))
                raise Exception ('Execution stopped because of some problem')
        else:
            try:
                data=session.xyDataListFromField(odb=session.odbs[macro_odbName], outputPosition=INTEGRATION_POINT, 
                variable=(('E', INTEGRATION_POINT, ((COMPONENT, strain_name), )), ), elementLabels=((element.instanceName,element.label),))
            except:
                pass
        #create the Amplitudes
        if source=='SDV':
            try: #try to create amplitude data for the selected integration point
                data=session.xyDataObjects['SDV'+str(i+1)+' PI: '+str(element.instanceName)+' E: '+str(element.label)+' IP: '+str(integration_point)]
                amplitude=mdb.models[para.Modelname].TabularAmplitude(name=strain_name,data=[tup for tup in data])
                print 'Amplitude for '+strain_name+' created.'
            except: #raise error if the integration point does not exist
                delete_xydata_SDV()
                WARNUNG=getWarningReply('Integration point does not exist!',(CANCEL,))
                raise Exception ('Execution stopped because of some problem')
                
        else:
            try: #try to create amplitude data for the selected integration point
                data=session.xyDataObjects['E:'+strain_name+' PI: '+str(element.instanceName)+' E: '+str(element.label)+' IP: '+str(integration_point)]
                amplitude=mdb.models[para.Modelname].TabularAmplitude(name=strain_name,data=[tup for tup in data])
                print 'Amplitude for '+strain_name+' created.'
            except:
                amplitude=mdb.models[para.Modelname].TabularAmplitude(name=strain_name,data=[(0.0,0.0),(1000000.0,0.0)])
                data_error=data_error+1
    
    if source!='SDV' and data_error==para.ndof_macro:
        WARNUNG=getWarningReply('No strain odb data was found!',(CANCEL,))
        raise Exception ('Execution stopped because of some problem')
    
    delete_xydata_SDV()
    
def delete_xydata_SDV():
    #delete all the created xy-Data
    getridkeys=session.xyDataObjects.keys()
    for i in range(len(getridkeys)):
        del session.xyDataObjects[getridkeys[i]] 

def unspecific_training_directions(n_variations,dimens,amplitude,Type,nlgeom):
    #this subroutine outputs monotonic training directions and varies each (stress or strain) entry n_variations times 
    
    variations=np.linspace(-1.0,1.0,n_variations)
    values=[]
    if dimens==2:
        for e11 in variations:
            for e22 in variations:
                for e12 in variations:
                    N=np.array([[e11,e12],[e12,e22]])
                    N_norm=np.linalg.norm(N)
                    if N_norm>0.00001:
                        if nlgeom=='YES' and Type=='strain':
                            values.append([N[0,0],N[1,0],0.0,N[0,1],N[1,1],0.0,0.0,0.0,0.0]/N_norm)
                        else:
                            values.append([N[0,0],N[1,1],0.0,N[0,1]]/N_norm)
    else:
        for e11 in variations:
            for e22 in variations:
                for e33 in variations:
                    for e12 in variations:
                        for e23 in variations:
                            for e13 in variations:
                                N=np.array([[e11,e12,e13],[e12,e22,e23],[e13,e23,e33]])
                                N_norm=np.linalg.norm(N)
                                if N_norm>0.00001:
                                    if nlgeom=='YES' and Type=='strain':
                                        values.append([N[0,0],N[1,0],N[2,0],N[0,1],N[1,1],N[2,1],N[0,2],N[1,2],N[2,2]]/N_norm)
                                    else:
                                        values.append([N[0,0],N[1,1],N[2,2],N[0,1],N[0,2],N[1,2]]/N_norm)
    
    values=np.array(values)*amplitude
    
    return values

def cluster_training_directions(values,prescribe_number_clusters,n_clusters,cluster_rel_error_drop):
    #this subroutine takes the stresses or strains (in columns) for every macroscopic GP
    #and performs k-means clustering from k=1...k=N Clusters until the deviation from the
    #actual data does not drop anymore by taking in more clusters
    
    if prescribe_number_clusters=='prescribed number of clusters:':
        cluster_centers,_=kmeans(values,n_clusters,iter=10) #compute cluster centers with fixed k
    else:
        errors=[] #collect the sum of squared errors in an array
        
        k=0
        while True:
            
            k=k+1 #increase number of clusters
            
            cluster_centers,_=kmeans(values,k,iter=10) #compute cluster centers
            
            codebook,distorsions=vq(values,cluster_centers) #assing data points
            
            errors.append(np.sum(distorsions**2)) #get the error
            
            if k>3: #if decline drops below 5% of decline rate stop clustering
                error_decline=(errors[k-1]-errors[k-3])/2.0
                if abs(error_decline)<abs(cluster_rel_error_drop*reference_error_decline):
                    break
            elif k==3: #set initial error decline as reference
                reference_error_decline=(errors[2]-errors[0])/2.0
            
    return cluster_centers

def get_training_data_from_odb(identifier,odb_for_clustering,timestep_numbers_clustering,GroupBox_clusteringmode,
                               prescribe_number_clusters,n_clusters,cluster_rel_error_drop):
    #this routine goes through all macro GPs and takes the stresses or strains from all desired timesteps and
    #saves the result in values, subsequently the data is clustered
    
    ODB=session.odbs[odb_for_clustering]
    
    if GroupBox_clusteringmode=='cluster all time steps': #culster all timesteps
        Frames=range(len(session.odbs[ODB.name].steps[ODB.steps.keys()[0]].frames))
        normalization=False
    elif GroupBox_clusteringmode=='cluster only speficific timesteps':
        Frames=[int(a) for a in timestep_numbers_clustering.split(',')]
        normalization=False
    elif GroupBox_clusteringmode=='cluster normalized values of the last timestep':
        Frames=[len(session.odbs[ODB.name].steps[ODB.steps.keys()[0]].frames)-1]
        normalization=True
        max_norm=0.0
    timesteps=[]
    
    NGP=len(session.odbs[ODB.name].steps[ODB.steps.keys()[0]].frames[0].fieldOutputs[identifier].values)
    values=[[] for i in range(NGP)] #get empty values array (stresses/strains from all GPs)
    
    for f in Frames:
        timesteps.append(session.odbs[ODB.name].steps[ODB.steps.keys()[0]].frames[f].frameValue)
        for i in range(NGP):
            if identifier=='LE': #get the right stretch from the logarithmic strain
                LE=session.odbs[ODB.name].steps[ODB.steps.keys()[0]].frames[f].fieldOutputs[identifier].values[i].data
                if len(LE)==4:
                    le=np.array([[   LE[0], 0.5*LE[3], 0.0],
                                [0.5*LE[3],     LE[1], 0.0],
                                [0.0,             0.0, LE[2]]])
                else:
                    le=np.array([   [LE[0], 0.5*LE[3], 0.5*LE[4]],
                                [0.5*LE[3],     LE[1], 0.5*LE[5]],
                                [0.5*LE[4], 0.5*LE[5],     LE[2]]])
                
                #decompose into eigenvalues and eigenvectors
                val,vec=np.linalg.eig(le)
                
                #get the right stretch tensor
                U=np.exp(val[0])*np.matmul(np.reshape(vec[:,0],(3,1)),np.reshape(vec[:,0],(1,3)))+\
                  np.exp(val[1])*np.matmul(np.reshape(vec[:,1],(3,1)),np.reshape(vec[:,1],(1,3)))+\
                  np.exp(val[2])*np.matmul(np.reshape(vec[:,2],(3,1)),np.reshape(vec[:,2],(1,3)))
                
                values[i].extend([U[0,0]-1.0,U[1,0],U[2,0],U[0,1],U[1,1]-1.0,U[2,1],U[0,2],U[1,2],U[2,2]-1.0])
                
                if normalization:
                    norm=np.linalg.norm(values[i])
                
            else: #get (cauchy) stress or (infinitesimal) strain
                values[i].extend(session.odbs[ODB.name].steps[ODB.steps.keys()[0]].frames[f].fieldOutputs[identifier].values[i].data)
                if normalization:
                    norm=np.linalg.norm(np.concatenate((values[i][:3],values[i][3:len(values[i])]/np.sqrt(2))))
            
            if normalization: #normalize the values (only last timestep!)
                if norm>max_norm: #get maximum of norms
                    max_norm=norm
                if norm>0.00000000001:
                    values[i]=values[i]/norm
                else:
                    values[i]=values[i]*0.0
    
    values=np.array(values)
    
    if normalization:
        values=values*max_norm
    
    cluster_centers=cluster_training_directions(values,prescribe_number_clusters,n_clusters,cluster_rel_error_drop)
    
    return cluster_centers,timesteps

def generate_training_inputfile(training_directions,training_mode,odb_for_clustering,
                                timestep_numbers_clustering,GroupBox_clusteringmode,
                                strain_amp,stress_amp,dimens,nlgeom,jobname,Method,
                                simulation_time,n_variations,dtime_data_dump,ncpus,
                                program_directory,start_simulations,
                                prescribe_number_clusters,n_clusters,
                                cluster_rel_error_drop):
    #this subroutine first gathers the training data according to the desired
    #training mode specified by the user and saves it to a inputfile for
    #UMAT_Driver
    
    strains=[]
    stresses=[]
    
    if training_directions=='clustering':
        if training_mode=='stress and strain' or training_mode=='strain':
            if nlgeom=='NO':
                strains,timesteps=get_training_data_from_odb('E',odb_for_clustering,timestep_numbers_clustering,
                                                            GroupBox_clusteringmode,prescribe_number_clusters,n_clusters,
                                                            cluster_rel_error_drop)
            else:
                strains,timesteps=get_training_data_from_odb('LE',odb_for_clustering,timestep_numbers_clustering,
                                                             GroupBox_clusteringmode,prescribe_number_clusters,n_clusters,
                                                             cluster_rel_error_drop)
        if training_mode=='stress and strain' or training_mode=='stress':
            stresses,timesteps=get_training_data_from_odb('S',odb_for_clustering,timestep_numbers_clustering,
                                                          GroupBox_clusteringmode,prescribe_number_clusters,n_clusters,
                                                            cluster_rel_error_drop)
        simulation_time=timesteps[-1]
    else:
        if training_mode=='stress and strain' or training_mode=='strain':
            strains=unspecific_training_directions(n_variations,dimens,strain_amp,'strain',nlgeom)
        if training_mode=='stress and strain' or training_mode=='stress':
            stresses=unspecific_training_directions(n_variations,dimens,stress_amp,'stress',nlgeom)
        timesteps=[1.0] #assume monotonic loading
    
    with open(jobname+'.inp',"w") as f:
        ######################## basic informations ##################################
        f.write('*Material\nMonolithFE2\n1\n1\n')
        f.write('*Depvar\n1\n')
        if dimens==2:
            f.write('*NDI\n3\n*NSHR\n1\n')
        else:
            f.write('*NDI\n3\n*NSHR\n3\n')
        #now define the steps
        for i in range(len(strains)):
            f.write('*Step\n*name\nTrainingstep-'+str(i+1)+'\n')
            f.write('*Nlgeom\n'+nlgeom+'\n')
            f.write('*Static\n0.01,'+str(timesteps[-1])+',0.00001,0.1\n')
            if nlgeom=='NO':
                f.write('*STRAN\n')
            else:
                f.write('*DFGRD\n')
            for j in range(len(timesteps)):
                f.write(str(timesteps[j])+'\n')
                for k in range(len(strains[i])/len(timesteps)):
                    if k>0:
                        f.write(',')
                    f.write(str(strains[i,len(strains[i])/len(timesteps)*j+k]))
                f.write('\n')
            f.write('*end_Step\n')
        for i in range(len(stresses)):
            f.write('*Step\n*name\nTrainingstep-'+str(len(strains)+i+1)+'\n')
            f.write('*Nlgeom\n'+nlgeom+'\n')
            f.write('*Static\n0.01,'+str(timesteps[-1])+',0.00001,0.1\n')
            f.write('*STRESS\n')
            for j in range(len(timesteps)):
                f.write(str(timesteps[j])+'\n')
                for k in range(len(stresses[i])/len(timesteps)):
                    if k>0:
                        f.write(',')
                    f.write(str(stresses[i,len(stresses[i])/len(timesteps)*j+k]))
                f.write('\n')
            f.write('*end_Step\n')
    
    if Method=="ROM":
        reduction='full'
    else:
        reduction='reduced'
    
    #set the Analysisparameters so that the training data is generated
    set_FE2_Analysisparameters('staggered',0.00001,12,'notsaved','indefinite','symm',reduction,
                               training_data='ROM',simulation_time=simulation_time,dtime_data_dump=dtime_data_dump)
    #start the training simulations
    if start_simulations:
        if platform.system()=='Linux':
            subprocess.Popen(program_directory+" cpus="+str(ncpus)+" job="+jobname, shell=True)
        else:
            WARNUNG=getWarningReply('On Windows currently directly starting UMAT_Driver is not possible.\n'+
                            'Please manually start:\n'
                            +program_directory+" cpus="+str(ncpus)+" job="+jobname,('OK',))
            raise Exception ('Execution stopped because of some problem')

def MonolithFE2Module(mode='',Modelname='',Partname='',jobname='',rve_number='',
            information='',PBC='',leftside_abq='',downside_abq='',rightside_abq='',
            upside_abq='',gen_constraints='',scheme='',convergence_ratio='',
            max_iters='',save_soe='',indefinite_matrix='',symmetric_matrix='',
            solving_process='',macro_odbName='',macro_modelName='',element='',
            integration_point='',source='',micro_modelName='',micro_partName='',
            ncpus='',Method='',training_directions='',training_mode='',nlgeom='',
            dtime_data_dump='',program_directory='',start_simulations='',
            odb_for_clustering='',part_for_clustering='',GroupBox_clusteringmode='',
            timestep_numbers_clustering='',n_variations='',strain_amp='',
            stress_amp='',simulation_time='',n_modes='',NGP='',path_to_inputfile='',
            prescribe_number_clusters='',n_clusters='',cluster_rel_error_drop=''):
            #='' makes all arguments optional
            #actual interface for the Abaqus Plugin MonolithFE2

    if not (mode=='Set Analysisparameters to control the Simulation' or
            mode=='Evaluate the training simulations'):
        try:
            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:
                WARNUNG=getWarningReply('Only 2D planar or 3D RVE Parts allowed!',(CANCEL,))
                raise Exception ('Execution stopped because of some problem')
                
        except:
            WARNUNG=getWarningReply('Select the micro model in the General Dialogue!',(CANCEL,))
            raise Exception ('Execution stopped because of some problem')
        
        para=meshparameters(dimens,Modelname,Partname) #create a meshparameters object
    
    if mode=='Generate a Inputscript for MonolithFE2':
        
        #-> get mesh coordinates and element to node connection
        para.get_data_from_abaqus()
    
        #-> get the connection between slave dof's and master dof's
        #call get_slave_to_master_connnection
        if para.dimens==2:
            para.get_slave_to_master_connnection2D(upside_abq, downside_abq,
                leftside_abq,rightside_abq)
        else:
            para.get_slave_to_master_connnection3D(PBC)
        
        #-> get connection between elements and dof's
        para.get_element_to_dof()
        
        #-> get connection element stiffness and k_matrix_values entries
        para.get_stiffnessmatrix_entries()
        
        #-> get "equations" for enforcing the constraint
        para.get_equations()
        
        #-> create "equations" in Abaqus CAE for Postprocessing
        if gen_constraints:
            para.create_equations_CAE(PBC)
        
        #-> create Inputscript, which is read by the FE2 Program
        para.write_data_to_file(jobname,rve_number,information)
    
    elif mode=='Set Analysisparameters to control the Simulation':
        #-> create the Analysisparameters file
        set_FE2_Analysisparameters(scheme,convergence_ratio,max_iters,save_soe,
                                indefinite_matrix,symmetric_matrix,solving_process)
    
    elif mode=='Generate Training data for ROM simulations':
        generate_training_inputfile(training_directions,training_mode,odb_for_clustering,
                             timestep_numbers_clustering,GroupBox_clusteringmode,
                             strain_amp,stress_amp,dimens,nlgeom,jobname,Method,
                             simulation_time,n_variations,dtime_data_dump,ncpus,
                             program_directory,start_simulations,prescribe_number_clusters,
                             n_clusters,cluster_rel_error_drop)
    elif mode=='Evaluate the training simulations':
        #-> evaluate the training data to get the ROM modes and hyperintegration points
        evaluate_training_data(path_to_inputfile,n_modes,NGP,ncpus,Method,from_stored_SVD=True)
    
    elif mode=='Extract postprocessing data for a resimulation':
        #write load history data of macro element to amplitudes
        generate_postprocessing_data(macro_odbName,element,integration_point,
                                    source,para)
