#=============================================================================
# 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",
#   DOI: https://doi.org/10.1016/j.cma.2021.113886
#   N. Lange, G. Huetter, B. Kiefer: "A monolithic hyper ROM FE2 method with
#                                     clustered training at finite deformations"
#   DOI: https://doi.org/10.1016/j.cma.2023.116522
#
# 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 itertools import product
from Hyperreduction import *

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.equations=[] #equations for enforcing the periodic boundary conditions (constraint)
        para.Element_definition=[] #n:[Abaqus element type name, real properties, number state variables per IP]
        para.COORDSGLOBAL=[] #coordinates of the Mesh
        para.RVE_Volume=0 #volume of the RVE (including holes, pores etc.)
        para.period_vectors=[] #vectors in the direction where the periodic RVE repeats
        para.reaction_force_dof=[] #connection macro value to dummy micro value -> with reaction force obtained
        para.additional_dof=[] #connection macro value to dummy micro value -> without reaction force obtained
        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)}}
        
    def get_data_from_abaqus(para):
        
        ################# define all UEL elementtypes ##########################
        
        #elementlibrary={(Abaqus):  JTYPE,NGP,NTENS,element_dof}
        elementlibrary={'CPE4':     {'JTYPE':2001,  'NGP':4,  'NTENS':4,    'NDOFN':2,  'element_dof':[[i,j] for i in range(1,4+1)  for j in range(1,2+1)]},
                        'CPE8':     {'JTYPE':2002,  'NGP':9,  'NTENS':4,    'NDOFN':2,  'element_dof':[[i,j] for i in range(1,8+1)  for j in range(1,2+1)]},
                        'CPS4':     {'JTYPE':2003,  'NGP':4,  'NTENS':3,    'NDOFN':2,  'element_dof':[[i,j] for i in range(1,4+1)  for j in range(1,2+1)]},
                        'C3D8':     {'JTYPE':2004,  'NGP':8,  'NTENS':6,    'NDOFN':3,  'element_dof':[[i,j] for i in range(1,8+1)  for j in range(1,3+1)]},
                        'C3D20':    {'JTYPE':2005,  'NGP':27, 'NTENS':6,    'NDOFN':3,  'element_dof':[[i,j] for i in range(1,20+1) for j in range(1,3+1)]},
                        'CPE4R':    {'JTYPE':2006,  'NGP':1,  'NTENS':4,    'NDOFN':2,  'element_dof':[[i,j] for i in range(1,4+1)  for j in range(1,2+1)]},
                        'CPE8R':    {'JTYPE':2007,  'NGP':4,  'NTENS':4,    'NDOFN':2,  'element_dof':[[i,j] for i in range(1,8+1)  for j in range(1,2+1)]},
                        'CPS4R':    {'JTYPE':2008,  'NGP':1,  'NTENS':3,    'NDOFN':2,  'element_dof':[[i,j] for i in range(1,4+1)  for j in range(1,2+1)]},
                        'CPS8':     {'JTYPE':2009,  'NGP':9,  'NTENS':3,    'NDOFN':2,  'element_dof':[[i,j] for i in range(1,8+1)  for j in range(1,2+1)]},
                        'CPS8R':    {'JTYPE':2010,  'NGP':4,  'NTENS':3,    'NDOFN':2,  'element_dof':[[i,j] for i in range(1,8+1)  for j in range(1,2+1)]},
                        'C3D20R':   {'JTYPE':2011,  'NGP':8,  'NTENS':6,    'NDOFN':3,  'element_dof':[[i,j] for i in range(1,20+1) for j in range(1,3+1)]},
                        'CPE3':     {'JTYPE':2012,  'NGP':1,  'NTENS':4,    'NDOFN':2,  'element_dof':[[i,j] for i in range(1,3+1)  for j in range(1,2+1)]},
                        'C3D4':     {'JTYPE':2013,  'NGP':1,  'NTENS':6,    'NDOFN':3,  'element_dof':[[i,j] for i in range(1,4+1)  for j in range(1,3+1)]},
                        'C3D10':    {'JTYPE':2014,  'NGP':4,  'NTENS':6,    'NDOFN':3,  'element_dof':[[i,j] for i in range(1,10+1) for j in range(1,3+1)]},
                        'CPS3':     {'JTYPE':2024,  'NGP':1,  'NTENS':3,    'NDOFN':2,  'element_dof':[[i,j] for i in range(1,3+1)  for j in range(1,2+1)]},
                        'CPE6':     {'JTYPE':2029,  'NGP':3,  'NTENS':4,    'NDOFN':2,  'element_dof':[[i,j] for i in range(1,6+1)  for j in range(1,2+1)]},
                        'CPS6':     {'JTYPE':2030,  'NGP':3,  'NTENS':3,    'NDOFN':2,  'element_dof':[[i,j] for i in range(1,6+1)  for j in range(1,2+1)]}}
        
        ################## create instance in CAE if not present ###############
        
        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')
        
        ################################## Elements ###########################
        
        #get information of all Elements
        ELEMENTS=mdb.models[para.Modelname].parts[para.Partname].elements
        
        ########################## Assignments ################################
        
        ASSIGNMENTS=mdb.models[para.Modelname].parts[para.Partname].sectionAssignments
        
        ################ get all necessary element definitions ################
                
        Element_to_definition={} #e: Element_definition_label_n
        for a in range(len(ASSIGNMENTS)):
            #section name
            S=mdb.models[para.Modelname].sections[ASSIGNMENTS[a].sectionName]
            try: #acessing the element definition
                props=mdb.models[para.Modelname].materials[S.material].userMaterial.mechanicalConstants
                n_depvar=mdb.models[para.Modelname].materials[S.material].depvar.n
            except:
                WARNUNG=getWarningReply('All elements must be assigned to a region having a material with a user material, not empty constants and depvar defined.',(CANCEL,))
                raise Exception ('Execution stopped because of some problem')
            for e in ASSIGNMENTS[a].getSet().elements: #go through all elements in the assignment
                if not ELEMENTS[e.label-1].type.name in elementlibrary:
                    #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')
                else:
                    NSTATEV=(n_depvar+elementlibrary[ELEMENTS[e.label-1].type.name]['NTENS'])*elementlibrary[ELEMENTS[e.label-1].type.name]['NGP']
                    definition={'JTYPE':elementlibrary[ELEMENTS[e.label-1].type.name]['JTYPE'],'PROPS':props,'NDOFN':elementlibrary[ELEMENTS[e.label-1].type.name]['NDOFN'],'NSTATEV':NSTATEV,'NTENS':elementlibrary[ELEMENTS[e.label-1].type.name]['NTENS'],'element_dof':elementlibrary[ELEMENTS[e.label-1].type.name]['element_dof'],'element_to_node':[]}
                    if not definition in para.Element_definition:
                        para.Element_definition.append(definition)
                    Element_to_definition[e.label-1]=para.Element_definition.index(definition) #get the connection element to element definition
        
        ################## get element-to node connection ######################
        
        #check if all elements have a element assignment
        for e in range(len(ELEMENTS)):
            if not e in Element_to_definition:
                #display error
                WARNUNG=getWarningReply('Some elements are not assigned to a section, or the elements are not labeled from 1 to n in steps of 1.',(CANCEL,))
                raise Exception ('Execution stopped because of some problem')
            else:
                para.Element_definition[Element_to_definition[e]]['element_to_node'].append([ELEMENTS[e].connectivity[i]+1 for i in range(len(ELEMENTS[e].connectivity))])
        
        ############################## 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
            para.COORDSGLOBAL.append([NODES[n].coordinates[i] for i in range(para.dimens)])
        
    def get_equations_PBC(para,set_boundary):
        #get the periodic boundary conditions by finding the node pairs using the set_boundary which contains all nodes at the boundary
        
        #test if the boundary set is in the part
        if not set_boundary in mdb.models[para.Modelname].parts[para.Partname].sets.keys():
            #display error
            WARNUNG=getWarningReply('The RVE boundary set was not found in the part.',(CANCEL,))
            raise Exception ('Execution stopped because of some problem')
        else:
            boundary_nodes=mdb.models[para.Modelname].parts[para.Partname].sets[set_boundary].nodes
        
        #create instance in CAE if not present
        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 these reference points:  strain (resp. stretch) E11,E22,E12 etc., and macro displacement U
        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='E'+str(para.ordering[para.dimens][i][0])+str(para.ordering[para.dimens][i][1]))
        r=a.referencePoints
        for i in range(para.ndof_macro):
            a.Set(name='E'+str(para.ordering[para.dimens][i][0])+str(para.ordering[para.dimens][i][1]),referencePoints=(r[refpoint_features[i]],))
        for i in range(para.dimens):
            refpoint_features.append(a.ReferencePoint(point=(0,0,0)).id)
            change_name=a.features.changeKey(fromName='RP-1',toName='U'+str(i+1))
        r=a.referencePoints
        for i in range(para.dimens):
            a.Set(name='U'+str(i+1),referencePoints=(r[refpoint_features[para.ndof_macro+i]],))
        
        #get the "length" of the RVE (just for tolerance purpose)
        X_min=boundary_nodes[0].coordinates[0]; X_max=boundary_nodes[0].coordinates[0]
        for N in boundary_nodes:
            if N.coordinates[0]>X_max:
                X_max=N.coordinates[0]
            if N.coordinates[0]<X_min:
                X_min=N.coordinates[0]
        tol=0.00001*(X_max-X_min)
        
        #loop through all nodes to get the slave master connection
        slaves=[] #list with all slaves
        master_to_slave={} #dictionary master_i:slave_1,slave_2,slave_3,...]
        coords={}
        for n in range(len(boundary_nodes)):
            if not boundary_nodes[n].label in slaves:
                coords[boundary_nodes[n].label]=boundary_nodes[n].coordinates
                for comb in product([-1.0,0.0,1.0],repeat=len(para.period_vectors)):
                    if not all([c==0.0 for c in comb]):
                        vector=np.zeros(para.dimens)
                        for i in range(len(para.period_vectors)):
                            vector=vector+para.period_vectors[i]*comb[i]
                        Vector=[0.0]*3
                        for i in range(para.dimens):
                            Vector[i]=boundary_nodes[n].coordinates[i]+vector[i]
                        slave=boundary_nodes.getClosest(tuple(Vector),numToFind=1,searchTolerance=tol)
                        if slave!=None:
                            if not slave.label in slaves:
                                slaves.append(slave.label)
                                coords[slave.label]=slave.coordinates
                                if boundary_nodes[n].label in master_to_slave:
                                    master_to_slave[boundary_nodes[n].label].append(slave.label)
                                else:
                                    master_to_slave[boundary_nodes[n].label]=[slave.label]
                if not boundary_nodes[n].label in master_to_slave:
                    WARNUNG=getWarningReply('The mesh is not periodic within tolerance.',(CANCEL,))
                    raise Exception ('Execution stopped because of some problem')
        
        #take one node where the rigid body motion is prevented and set it to zero
        n=0 #take the master with the most slaves -> efficiency
        for master in master_to_slave:
            if len(master_to_slave[master])>n:
                n=len(master_to_slave[master])
                zero_displacement_master=master
        
        #equation for zero displacement
        a.SetFromNodeLabels(name='zero_displacement_node',nodeLabels=((instance_name,(zero_displacement_master,)),))
        for i in range(1,para.dimens+1):
            mdb.models[para.Modelname].Equation(name='rigid_body_motion'+str(i),terms=((-1.0,'zero_displacement_node',i),(1.0,'U'+str(i),1)))
        
        #loop through all pairs of slave/masters
        for master in master_to_slave:
            #create sets with master node
            if not master==zero_displacement_master:
                master_setname='Set-m-'+str(master)
                a.SetFromNodeLabels(name=master_setname,nodeLabels=((instance_name,(master,)),))
            for slave in master_to_slave[master]:
                slave_setname='Set-s-'+str(slave)
                a.SetFromNodeLabels(name=slave_setname,nodeLabels=((instance_name,(slave,)),))
                for i in range(1,para.dimens+1):
                    #assigne equations u_i^+ = u_i^- + H_{ij} {x^+ - x^-}
                    if master==zero_displacement_master:
                        Terms=[(-1.0,slave_setname,i)]
                    else:
                        Terms=[(-1.0,slave_setname,i),(1.0,master_setname,i)]
                    for j in range(1,para.dimens+1):
                        delta_x=coords[slave][j-1]-coords[master][j-1]
                        if abs(delta_x)>tol:
                            Terms.append((delta_x*para.ordering[para.dimens][i,j][1],'E'+str(para.ordering[para.dimens][para.ordering[para.dimens][i,j][0]][0])+str(para.ordering[para.dimens][para.ordering[para.dimens][i,j][0]][1]),1))
                    mdb.models[para.Modelname].Equation(name='Eqn-'+'m'+str(master)+'-s'+str(slave)+'-'+str(i), 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
        for k in range(para.dimens):
            mdb.models[para.Modelname].DisplacementBC(name='rigid_body_motion_'+str(k), createStepName=stepname,region=regionToolset.Region(referencePoints=(r[refpoint_features[para.ndof_macro+k]],)), u1=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
        
        ################### read the equations from Abaqus CAE #################
        
        a=mdb.models[para.Modelname].rootAssembly
        p=mdb.models[para.Modelname].parts[para.Partname]
        
        #give connection between micro dummy node and macro value
        dic_macro_values_to_micro_node_label={}
        para.COORDSGLOBAL.append([0.0]*para.dimens) #[U1,U2,(U3)]
        for i in range(1,para.dimens+1):
            dic_macro_values_to_micro_node_label['U'+str(i)]=[len(para.COORDSGLOBAL),i]
            para.additional_dof.append([len(para.COORDSGLOBAL),i])
        node=len(para.COORDSGLOBAL)+1
        dof=1
        para.COORDSGLOBAL.append([0.0]*para.dimens) #[E11,E22,(E33)]
        para.COORDSGLOBAL.append([0.0]*para.dimens) #[E11,(E22,E33)]
        for j in range(para.ndof_macro):
            dic_macro_values_to_micro_node_label['E'+str(para.ordering[para.dimens][j][0])+str(para.ordering[para.dimens][j][1])]=[node,dof]
            para.reaction_force_dof.append([node,dof])
            if dof==para.dimens:
                dof=1
                node=node+1
            else:
                dof=dof+1
        
        deleted_dofs=[]
        for key in list(mdb.models[para.Modelname].constraints.keys()):
            para.equations.append([])
            #check the equation
            deleted_dofs.append(mdb.models[para.Modelname].constraints[key].terms[0][1])
            if any([dd in deleted_dofs for dd in [mdb.models[para.Modelname].constraints[key].terms[i][1] for i in range(1,len(mdb.models[para.Modelname].constraints[key].terms))]]):
                WARNUNG=getWarningReply('Error in Equation '+key+', each set shall contain one node resp. one macro Reference Point.',(CANCEL,))
                raise Exception ('Execution stopped because of some problem')
            #put the Equation into the format of MonolithFE2
            for term in mdb.models[para.Modelname].constraints[key].terms:
                if term[1] in dic_macro_values_to_micro_node_label:
                    para.equations[-1].append([dic_macro_values_to_micro_node_label[term[1]],term[0]])
                else:
                    if len(mdb.models[para.Modelname].rootAssembly.sets[term[1]].nodes)==1:
                        if mdb.models[para.Modelname].rootAssembly.sets[term[1]].nodes[0].label<=len(para.COORDSGLOBAL):
                            para.equations[-1].append([mdb.models[para.Modelname].rootAssembly.sets[term[1]].nodes[0].label,term[2],term[0]])
                        else:
                            WARNUNG=getWarningReply('Error in Equation '+key+', a bad node label appeared in set'+term[1]+'.',(CANCEL,))
                            raise Exception ('Execution stopped because of some problem')
                    else:
                        WARNUNG=getWarningReply('Error in Equation '+key+', each set shall contain one node resp. one macro Reference Point.',(CANCEL,))
                        raise Exception ('Execution stopped because of some problem')
        
        ####################### write file #####################################
        
        direct_path=os.path.abspath(".")
        
        fobj_out = open(direct_path+"/"+jobname+".FE"+str(rve_number),"w")
        
        #output version number
        fobj_out.write('**MonolithFE2 Version 3.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('*RVE_Volume, V='+para.RVE_Volume+'\n')
        
        fobj_out.write('*Coupling, Ndof=')
        fobj_out.write(str(para.Element_definition[0]['NDOFN']))
        fobj_out.write('\n')
        for i in range(para.Element_definition[0]['NDOFN']):
            for j in str([int(1)]*para.Element_definition[0]['NDOFN']):
                if not(j=='[' or j==']'):
                    fobj_out.write(j)
            fobj_out.write('\n')
        
        fobj_out.write('*User_Elements, N=')
        fobj_out.write(str(len(para.Element_definition)))
        fobj_out.write('\n')
        for i in range(len(para.Element_definition)):
            fobj_out.write('*Element'+str(i+1)+', TYPE='+str(para.Element_definition[i]['JTYPE'])+'\n')
            fobj_out.write('*Element_dof, N='+str(len(para.Element_definition[i]['element_dof']))+'\n')
            for k in para.Element_definition[i]['element_dof']:
                for j in str(k):
                    if not(j=='[' or j==']'):
                        fobj_out.write(j)
                fobj_out.write('\n')
            fobj_out.write('*NNODE, N='+str(len(para.Element_definition[i]['element_to_node'][0]))+'\n')
            fobj_out.write('*NSVARS, N='+str(para.Element_definition[i]['NSTATEV'])+'\n')
            fobj_out.write('*PROPS, N='+str(len(para.Element_definition[i]['PROPS']))+'\n')
            for j in str(para.Element_definition[i]['PROPS']):
                if not(j=='[' or j==']'):
                    fobj_out.write(j)
            fobj_out.write('\n')
            fobj_out.write('*JPROPS, N=0\n') #assume zero integer element properties, could be set manually in the inputfile
            fobj_out.write('*n_additional_hyper_outputs, N='+str(para.Element_definition[i]['NTENS']+1)+'\n') #output stress and stress power
            fobj_out.write('*End_of_Element\n')
        fobj_out.write('*End_of_User_Elements\n')
        
        fobj_out.write('*Node, dimens='+str(para.dimens)+'\n')
        for i in para.COORDSGLOBAL:
            for j in str(i):
                if not(j=='[' or j==']'):
                    fobj_out.write(j)
            fobj_out.write('\n')
        
        fobj_out.write('*Element_to_Node, N='+str(len(para.Element_definition))+'\n')
        for i in range(len(para.Element_definition)):
            fobj_out.write('*Element_to_Node_Assignments'+str(i+1)+', N='+str(len(para.Element_definition[i]['element_to_node']))+'\n')
            for j in para.Element_definition[i]['element_to_node']:
                for k in str(j):
                    if not(k=='[' or k==']'):
                        fobj_out.write(k)
                fobj_out.write('\n')
            fobj_out.write('*End_of_Element_to_Node_Assignments\n')
        fobj_out.write('*End_of_Element_to_Node\n')
        
        fobj_out.write('*Reaction_Force_dof, N='+str(len(para.reaction_force_dof))+'\n')
        for i in para.reaction_force_dof:
            for j in str(i):
                if not(j=='[' or j==']'):
                    fobj_out.write(j)
            fobj_out.write('\n')
        
        fobj_out.write('*Additional_dof, N='+str(len(para.additional_dof))+'\n')
        for i in para.additional_dof:
            for j in str(i):
                if not(j=='[' or j==']'):
                    fobj_out.write(j)
            fobj_out.write('\n')
        
        fobj_out.write('*Equations, N='+str(len(para.equations))+'\n')
        for i in range(len(para.equations)):
            fobj_out.write('*EQUATION'+str(i+1)+', N='+str(len(para.equations[i]))+'\n')
            for j in str(para.equations[i]):
                if not(j=='[' or j==']'):
                    fobj_out.write(j)
            fobj_out.write('\n')
        fobj_out.write('*End_of_Equations\n')
        
        fobj_out.write('*End_of_File')
        
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
    
    #get number of shear components NSHR
    if len(strains)>0:
        NSHR=len(strains[1])-3
    else:
        NSHR=len(stresses[1])-3
    
    with open(jobname+'.inp',"w") as f:
        ######################## basic informations ##################################
        f.write('*Material\nMonolithFE2\n1\n1\n')
        f.write('*Depvar\n1\n')
        f.write('*NDI\n3\n*NSHR\n'+str(NSHR)+'\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 set_FE2_Analysisparameters(scheme,convergence_ratio,max_iters,save_soe,
                               indefinite_matrix,symmetric_matrix,solving_process,
                               training_data='NO',simulation_time=1.0,dtime_data_dump=0.1):
    
    direct_path=os.path.abspath(".")
    fobj_out = open(direct_path+"/"+"FE2_Analysisparameters.cfg","w")
    
    #output Version number
    fobj_out.write('**MonolithFE2 Version 3.0\n')
    
    #general warning
    fobj_out.write('**Do not modify this file, otherwise it may not work properly!\n')
    #staggered/monolithic algorithm
    fobj_out.write('*Algorithm\n')
    if scheme=='staggered':
        fobj_out.write('staggered\n')
    else:
        fobj_out.write('monolithic\n')
    #save factorization in monolithic scheme?
    if not(scheme=='staggered'):
        fobj_out.write('*Factorization\n')
        if save_soe=='FALSE':
            fobj_out.write('notsaved')
        else:
            fobj_out.write('saved')
        fobj_out.write('\n')
    #max. number of equilibrium iterations in staggered scheme
    else:
        fobj_out.write('*Equilibrium_iterations\n')
        if max_iters>0 and max_iters<21:
            fobj_out.write(str(max_iters))
            fobj_out.write('\n')
        else:
            raise Exception('max. number of iterations in staggered scheme is supposed to be between 1 and 20!')
        fobj_out.write('*Convergence_ratio\n')
        s = '%1.10f' % convergence_ratio
        if convergence_ratio>0.0 and convergence_ratio<1.0:
            fobj_out.write(s)
            fobj_out.write('\n')
        else:
            raise Exception('convergence_ratio is supposed to be between 0.0 and 1.0!')
    #indefinite matrix?
    fobj_out.write('*Systemmatrix_definiteness\n')
    if indefinite_matrix=='indefinite':
        fobj_out.write('indefinite')
    else:
        fobj_out.write('positiv')
    fobj_out.write('\n')
    
    #symmetric matrix?
    fobj_out.write('*Systemmatrix_symmetry\n')
    if symmetric_matrix=='symm':
        fobj_out.write('symm')
    else:
        fobj_out.write('unsymm')
    fobj_out.write('\n')
    
    #solving process restricted by ROM modes or full?
    fobj_out.write('*Solving_Process\n')
    if solving_process=='full':
        fobj_out.write('full')
    elif solving_process=='reduced':
        fobj_out.write('reduced')
    else:
        fobj_out.write('hyperreduced')
    fobj_out.write('\n')
        
    #generate training data or not
    fobj_out.write('*Training_data\n')
    if training_data=='ROM':
        fobj_out.write('ROM\n')
        #dump the data at equally distributed time slots
        n_data_dumps=int(simulation_time/dtime_data_dump)
        fobj_out.write('*Data_dump, N='+str(n_data_dumps+1)+'\n')
        for n in range(n_data_dumps+1):
            fobj_out.write(str(dtime_data_dump*float(n)))
            if n!=n_data_dumps:
                fobj_out.write(',')
            else:
                fobj_out.write('\n')
    else:
        fobj_out.write('NO')

    fobj_out.close()

def MonolithFE2Module(mode='',Modelname='',Partname='',RVE_Volume='',jobname='',rve_number='',
            information='',PBC='',set_boundary='',period_vectors='',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='',NELEM='',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':
        
        #do some consistency checks
        if not float(RVE_Volume)>0.0:
            WARNUNG=getWarningReply('RVE Volume has to be greater than Zero.',(CANCEL,))
            raise Exception ('Execution stopped because of some problem')
        else:
            para.RVE_Volume=RVE_Volume
        
        if not information.strip():
            WARNUNG=getWarningReply('RVE description shall not be empty.',(CANCEL,))
            raise Exception ('Execution stopped because of some problem')
        
        #-> get mesh coordinates and element to node connection
        para.get_data_from_abaqus()
    
        #-> get periodic boundary conditions
        if PBC=='Create Periodic BCs':
            #get the periodicity vectors
            try:
                para.period_vectors=[np.fromstring(vec,dtype=float,sep=',') for vec in period_vectors.split(';')]
            except:
                WARNUNG=getWarningReply('Could not understand the periodicity vectors.',(CANCEL,))
                raise Exception ('Execution stopped because of some problem')
            if any([len(vec)!=para.dimens for vec in para.period_vectors]):
                WARNUNG=getWarningReply('The dimension of the periodicity vectors shall be identical to the dimension of the problem.',(CANCEL,))
                raise Exception ('Execution stopped because of some problem')
            para.get_equations_PBC(set_boundary)
        
        #-> 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,NELEM,ncpus,Method)
    
    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)
