# -*- coding: mbcs -*-
"""
Incremental post-processing of 3D HSK contact analyses

This Abaqus script evaluates a completed 3D ODB and extracts the quantities
used in the numerical investigations of Chapter 4. The script reads contact
pressure along the plane and taper interfaces, computes integrated contact
forces for the available load states, and exports the aggregated as well as
incremental results to text files.

Parameter interface
-------------------
g_Jobname : str
    Abaqus job name provided by a calling batch script. A trailing line break is
    tolerated.

Notes
-----
The computational workflow and exported quantities were retained. Revisions are
limited to readability, sectioning and removal of redundant commented code.
"""

from part import *
from material import *
from section import *
from assembly import *
from step import *
from interaction import *
from load import *
from mesh import *
from optimization import *
from job import *
from sketch import *
from visualization import *
from connectorBehavior import *
from operator import itemgetter
import math
import operator
import sys

import displayGroupOdbToolset as dgo
import numpy as np
import odbAccess
from abaqusConstants import *
from odbAccess import *
import __main__

# ---------------------------------------------------------------------------
# Evaluation settings and ODB access
# ---------------------------------------------------------------------------

reib = 0.15
odbName = 'HexMC_3d.odb'
myOdb = openOdb(odbName, readOnly=True)

# ---------------------------------------------------------------------------
# Create evaluation paths
# ---------------------------------------------------------------------------
#Auswertepfade erstellen
Liste_Flaechenpressung=[]
Liste_Kegel=[]
Liste_Kegel_coordinates = []
Liste_Kegel_Halter=[]
Zaehler=0
Laufkoordinate=0.0
Laufknotennummer=0

#Pfad Flaechenpressung
for Zaehler in range(len(myOdb.rootAssembly.nodeSets['FLAECHENPRESSUNG-PLAN'].nodes[0])):
    Laufkoordinate=myOdb.rootAssembly.nodeSets['FLAECHENPRESSUNG-PLAN'].nodes[0][Zaehler].coordinates[0]
    Laufknotennummer=myOdb.rootAssembly.nodeSets['FLAECHENPRESSUNG-PLAN'].nodes[0][Zaehler].label
    Liste_Flaechenpressung.append((Laufkoordinate, Laufknotennummer))
 
Zaehler=0
Laufkoordinate=0.0
Laufknotennummer=0  
 
#Pfad Kegel Aufnahme
for Zaehler in range(len(myOdb.rootAssembly.nodeSets['KEGEL'].nodes[0])):
    Laufkoordinate=myOdb.rootAssembly.nodeSets['KEGEL'].nodes[0][Zaehler].coordinates[1]
    Laufknotennummer=myOdb.rootAssembly.nodeSets['KEGEL'].nodes[0][Zaehler].label
    Liste_Kegel.append((Laufkoordinate, Laufknotennummer))

Zaehler=0
Laufkoordinate=0.0
Laufknotennummer=0 

#Pfad Kegel Halter
for Zaehler in range(len(myOdb.rootAssembly.nodeSets['KEGEL_HALTER'].nodes[0])):
    Laufkoordinate=myOdb.rootAssembly.nodeSets['KEGEL_HALTER'].nodes[0][Zaehler].coordinates[1]
    Laufkoordinate=Laufkoordinate*-1
    Laufknotennummer=myOdb.rootAssembly.nodeSets['KEGEL_HALTER'].nodes[0][Zaehler].label
    Liste_Kegel_Halter.append((Laufkoordinate, Laufknotennummer))

Liste_Flaechenpressung_sortiert=sorted(Liste_Flaechenpressung, key=itemgetter(0))
Liste_Kegel_sortiert=sorted(Liste_Kegel, key=itemgetter(0))
Liste_Kegel_Halter_sortiert=sorted(Liste_Kegel_Halter, key=itemgetter(0))


#Auslesen der Sortierten Knoten
Knotenliste_Flaechenpressung_sortiert=zip(*Liste_Flaechenpressung_sortiert)[1]
Knotenliste_Kegel_sortiert=zip(*Liste_Kegel_sortiert)[1]
Knotenliste_Kegel_Halter_sortiert=zip(*Liste_Kegel_Halter_sortiert)[1]

# print('Anzahl der Knoten auf der Planflaeche: %F '% (len(Knotenliste_Flaechenpressung_sortiert)))
# print('Anzahl der Knoten auf dem Kegel Halterseite: %F '% (len(Knotenliste_Kegel_Halter_sortiert)))
# print('Anzahl der Knoten auf dem Kegel: %F '% (len(Knotenliste_Kegel_sortiert)))

#Auswertepfad erzeugen
session.Path(name = 'Planflaeche', type = NODE_LIST, expression =('AUFNAHME_PARAMETER', Knotenliste_Flaechenpressung_sortiert))
path_Plan = session.paths['Planflaeche']
session.Path(name = 'Kegelflaeche', type = NODE_LIST, expression =('AUFNAHME_PARAMETER', Knotenliste_Kegel_sortiert))
path_Kegel = session.paths['Kegelflaeche']
session.Path(name = 'Kegelflaeche Halter', type = NODE_LIST, expression =('HALTER_PARAMETER', Knotenliste_Kegel_Halter_sortiert))
path_Kegel_Halter = session.paths['Kegelflaeche Halter']
#print('Auswertepfade erzeugt')

#alte Plots loeschen
#del session.xyDataObjects['RB-P']
#del session.xyDataObjects['RB-K']
#del session.xyDataObjects['AK-P']
#del session.xyDataObjects['AK-K']

#Create xy-Data
#Flaechenpressung PA Vorspannkraft
odb = session.odbs[odbName]
session.viewports['Viewport: 1'].setValues(displayedObject=odb)
session.viewports['Viewport: 1'].odbDisplay.setPrimaryVariable(variableLabel='CPRESS', outputPosition=ELEMENT_NODAL)
x = len(myOdb.steps['Vorspannkraft'].frames)
if x < 1:
	Ergebnis = 0
else: 
	Ergebnis = 1
	session.viewports['Viewport: 1'].odbDisplay.setFrame(step='Vorspannkraft', frame=x-1)
	pth = session.paths['Planflaeche']
	session.XYDataFromPath(name='FP_PA_0Umin', path=pth, includeIntersections=False, 
		projectOntoMesh=False, pathStyle=PATH_POINTS, numIntervals=10, 
		projectionTolerance=0, shape=UNDEFORMED, labelType=X_COORDINATE, 
		removeDuplicateXYPairs=True, includeAllElements=False)
	

print('Flachenpressung_PA_0Umin')

F_PA_ink_0Umin = []
x_val_1_pa = []
x_val_2_pa = []
p_val_pa_0 = []
FP_PA_0 = session.xyDataObjects['FP_PA_0Umin']
x = len(FP_PA_0)
for i in range(0,x-1,1):
	val_1 = FP_PA_0[i]
	val_2 = FP_PA_0[i+1]
	x_val_1 = val_1[0]
	x_val_1_pa.append(x_val_1)
	x_val_2 = val_2[0]
	x_val_2_pa.append(x_val_2)
	p_val = val_2[1]
	p_val_pa_0.append(p_val)
	F_val = (x_val_2*x_val_2 - x_val_1*x_val_1)*3.14*p_val
	F_PA_ink_0Umin.append(F_val)

F_PA_ges_0Umin=sum(F_PA_ink_0Umin)	

#Flaechenpressung PA Radialbeschleunigung
# x = len(myOdb.steps['Radialbeschleunigung'].frames)
# if x < 1:
	# Ergebnis = 0
# else: 
	# Ergebnis = 1
	# session.viewports['Viewport: 1'].odbDisplay.setFrame(step='Radialbeschleunigung', frame=x-1)
	# pth = session.paths['Planflaeche']
	# session.XYDataFromPath(name='FP_PA_25kUmin', path=pth, includeIntersections=False, 
		# projectOntoMesh=False, pathStyle=PATH_POINTS, numIntervals=10, 
		# projectionTolerance=0, shape=UNDEFORMED, labelType=X_COORDINATE, 
		# removeDuplicateXYPairs=True, includeAllElements=False)

F_PA_ink_25kUmin = []
# FP_PA_25 = session.xyDataObjects['FP_PA_25kUmin']
# p_val_pa_25 = []
# x = len(FP_PA_25)
# for i in range(0,x-1,1):
	# val_1 = FP_PA_25[i]
	# val_2 = FP_PA_25[i+1]
	# x_val_1 = val_1[0]
	# x_val_2 = val_2[0]
	# p_val = val_2[1]
	# p_val_pa_25.append(p_val)
	# F_val = (x_val_2*x_val_2 - x_val_1*x_val_1)*3.14*p_val
	# F_PA_ink_25kUmin.append(F_val)

F_PA_ges_25kUmin=sum(F_PA_ink_25kUmin)	

##Flaechenpressung Konus Vorspannkraft
odb = session.odbs[odbName]
session.viewports['Viewport: 1'].setValues(displayedObject=odb)
session.viewports['Viewport: 1'].odbDisplay.setPrimaryVariable(variableLabel='CPRESS', outputPosition=ELEMENT_NODAL)
x = len(myOdb.steps['Vorspannkraft'].frames)
if x < 1:
	Ergebnis = 0
else: 
	Ergebnis = 1
	session.viewports['Viewport: 1'].odbDisplay.setFrame(step='Vorspannkraft', frame=x-1)
	pth = session.paths['Kegelflaeche']
	session.XYDataFromPath(name='FP_K_0Umin', path=pth, includeIntersections=False, 
		projectOntoMesh=False, pathStyle=PATH_POINTS, numIntervals=10, 
		projectionTolerance=0, shape=UNDEFORMED, labelType=Y_COORDINATE, 
		removeDuplicateXYPairs=True, includeAllElements=False)

Knoten_KONUS_SPINDEL = myOdb.rootAssembly.nodeSets['KEGEL'].nodes[0]
Liste_Knoten_KONUS_SPINDEL = []
for i in Knoten_KONUS_SPINDEL:
	Liste_Knoten_KONUS_SPINDEL.append((i.label, i.coordinates[0], i.coordinates[1], i. coordinates[2]))
Liste_Knoten_KONUS_SPINDEL_sort = sorted(Liste_Knoten_KONUS_SPINDEL, key = itemgetter(1))

F_K_ink_0Umin = []
x_val_1_K = []
x_val_2_K = []
y_val_1_K = []
y_val_2_K = []
p_val_K_0 = []
FP_K_0 = session.xyDataObjects['FP_K_0Umin']
x = len(FP_K_0)
for i in range(0,x-1,1):
	val_1 = FP_K_0[i]
	val_2 = FP_K_0[i+1]
	coordinates_1 = Liste_Knoten_KONUS_SPINDEL_sort[i]
	coordinates_2 = Liste_Knoten_KONUS_SPINDEL_sort[i+1]
	x_val_1 = coordinates_1[1]
	x_val_1_K.append(x_val_1)
	x_val_2 = coordinates_2[1]
	x_val_2_K.append(x_val_2)
	y_val_1 = coordinates_1[2]
	y_val_1_K.append(y_val_1)
	y_val_2 = coordinates_2[2]
	y_val_2_K.append(y_val_2)
	h = abs(y_val_2-y_val_1)
	d_r = x_val_2-x_val_1
	m = sqrt(d_r ** 2 + h ** 2)
	A_m = (x_val_1 + x_val_2) * m * 3.14
	p_val = val_2[1]
	p_val_K_0.append(p_val)
	F_val = A_m*p_val
	F_K_ink_0Umin.append(F_val)

F_K_N_ges_0Umin=sum(F_K_ink_0Umin)
	
#Anteile Normalkraft
F_K_N_y_0Umin = cos(radians(2.8625)) * F_K_N_ges_0Umin
F_K_N_x_0Umin = sin(radians(2.8625)) * F_K_N_ges_0Umin

F_K_R_x_0Umin = cos(radians(2.8625)) * F_K_N_ges_0Umin * reib
F_K_x_ges_0Umin = F_K_N_x_0Umin  + F_K_R_x_0Umin
F_x_ges_0Umin = F_PA_ges_0Umin + F_K_N_x_0Umin + F_K_R_x_0Umin
F_lost_percentage_0Umin = (18000-F_x_ges_0Umin)/18000
F_ratio_0Umin_1 = F_PA_ges_0Umin/F_x_ges_0Umin
F_ratio_0Umin_2 = F_PA_ges_0Umin/18000

#Flaechenpressung Konus Radialbeschleunigung
# x = len(odb.steps['Radialbeschleunigung'].frames)
# if x < 1:
	# Ergebnis = 0
# else: 
	# Ergebnis = 1
	# session.viewports['Viewport: 1'].odbDisplay.setFrame(step='Radialbeschleunigung', frame=x-1)
	# pth = session.paths['Kegelflaeche']
	# session.XYDataFromPath(name='FP_K_25Umin', path=pth, includeIntersections=False, 
		# projectOntoMesh=False, pathStyle=PATH_POINTS, numIntervals=10, 
		# projectionTolerance=0, shape=UNDEFORMED, labelType=Y_COORDINATE, 
		# removeDuplicateXYPairs=True, includeAllElements=False)

F_K_ink_25kUmin = []
p_val_K_25 = []
# FP_K_25 = session.xyDataObjects['FP_K_25Umin']
# x = len(FP_K_25)
# for i in range(0,x-1,1):
	# val_1 = FP_K_25[i]
	# val_2 = FP_K_25[i+1]
	# coordinates_1 = Liste_Knoten_KONUS_SPINDEL_sort[i]
	# coordinates_2 = Liste_Knoten_KONUS_SPINDEL_sort[i+1]
	# x_val_1 = coordinates_1[1]
	# x_val_2 = coordinates_2[1]
	# y_val_1 = coordinates_1[2]
	# y_val_2 = coordinates_2[2]
	# h = abs(y_val_2-y_val_1)
	# d_r = x_val_2-x_val_1
	# m = sqrt(d_r ** 2 + h ** 2)
	# A_m = (x_val_1 + x_val_2) * m * 3.14
	# p_val = val_2[1]
	# p_val_K_25.append(p_val)
	# F_val = A_m*p_val
	# F_K_ink_25kUmin.append(F_val)

F_K_N_ges_25kUmin=sum(F_K_ink_25kUmin)	
F_K_N_y_25kUmin = cos(radians(2.8625)) * F_K_N_ges_25kUmin
F_K_N_x_25kUmin = sin(radians(2.8625)) * F_K_N_ges_25kUmin
F_K_R_x_25kUmin = cos(radians(2.8625)) * F_K_N_ges_25kUmin * reib
F_K_x_ges_25kUmin = F_K_N_x_25kUmin  + F_K_R_x_25kUmin
F_K_x_ges_25kUmin = F_K_N_x_25kUmin  + F_K_R_x_25kUmin
F_x_ges_25kUmin = F_PA_ges_25kUmin + F_K_N_x_25kUmin + F_K_R_x_25kUmin
F_lost_percentage_25kUmin = (18000-F_x_ges_25kUmin)/18000
F_ratio_25kUmin_1 = 1 #F_PA_ges_25kUmin/F_x_ges_25kUmin
F_ratio_25kUmin_2 = F_PA_ges_25kUmin/18000


### Torsionsmoment bzw. Verdrehung auswerten
# odb = session.odbs[odbName]
# session.viewports['Viewport: 1'].setValues(displayedObject=odb)
# session.viewports['Viewport: 1'].odbDisplay.setPrimaryVariable(variableLabel='UR', outputPosition=NODAL, refinement=(COMPONENT, 'UR2'), )
# session.viewports['Viewport: 1'].odbDisplay.setPrimaryVariable(variableLabel='CPRESS', outputPosition=ELEMENT_NODAL)

# xydata_temp = session.xyDataListFromField(odb=odb, outputPosition=NODAL, variable=((
    # 'UR', NODAL, ((COMPONENT, 'UR2'), )), ), nodeSets=("REFPOINT", ))
# xy_data = session.xyDataObjects['UR:UR2 PI: ASSEMBLY N: 1']
# x = len(xy_data)
# t_rot = []
# UR2 = []
# for i in range(0,x-1,1):
	# data = xy_data[i]
	# x = data[0]
	# y = data[1]
	# t_rot.append(x)
	# UR2.append(y)
	
# UR2_max = max(UR2)	
##close current odb
# del session.xyDataObjects['UR:UR2 PI: ASSEMBLY N: 1']
# myOdb.close()

# odb = session.odbs[odbName]
# session.viewports['Viewport: 1'].setValues(displayedObject=odb)
# session.viewports['Viewport: 1'].odbDisplay.setPrimaryVariable(variableLabel='UR', outputPosition=NODAL, refinement=(COMPONENT, 'UR2'), )
# #session.viewports['Viewport: 1'].odbDisplay.setPrimaryVariable(variableLabel='CPRESS', outputPosition=ELEMENT_NODAL)

# xydata_temp = session.xyDataListFromField(odb=odb, outputPosition=NODAL, variable=((
    # 'U', NODAL, ((COMPONENT, 'U1'), )), ), nodeSets=("REFPOINT", ))
# xy_data = session.xyDataObjects['U:U1 PI: ASSEMBLY N: 1']
# x = len(xy_data)
# t_rot = []
# UR2 = []
# for i in range(0,x-1,1):
	# data = xy_data[i]
	# x = data[0]
	# y = data[1]
	# t_rot.append(x)
	# UR2.append(y)
	
# UR2_max = max(UR2)	
# ##close current odb
# del session.xyDataObjects['U:U1 PI: ASSEMBLY N: 1']
# myOdb.close()



#save data
odbName = odbName[0:-4]

parameter = ['F_PA_ges_0Umin', 'F_PA_ges_25kUmin', 'F_K_N_ges_0Umin', 'F_K_N_y_0Umin', 'F_K_N_x_0Umin', 'F_K_R_x_0Umin',  'F_K_x_ges_0Umin', 'F_x_ges_0Umin','F_lost_percentage_0Umin',
'F_ratio_0Umin_1', 'F_ratio_0Umin_2', 'F_K_N_ges_25kUmin', 'F_K_N_y_25kUmin', 'F_K_N_x_25kUmin', 'F_K_R_x_25kUmin', 'F_K_x_ges_25kUmin', 'F_x_ges_25kUmin', 'F_lost_percentage_25kUmin', 'F_ratio_25kUmin_1','F_ratio_25kUmin_2', 'UR2_max']

x_val_1_pa = []
x_val_2_pa = []
p_val_pa_0 = []
p_val_pa_25 = []
x_val_1_K = []
x_val_2_K = []
y_val_1_K = []
y_val_2_K = []
p_val_K_0 = []
p_val_K_25 = []

data_ink_str = ['F_PA_ink_0Umin','F_PA_ink_25kUmin','F_K_ink_0Umin','F_K_ink_25kUmin','x_val_1_pa','x_val_2_pa','p_val_pa_0','p_val_pa_25',
'x_val_1_K','x_val_2_K','y_val_1_K','y_val_2_K','p_val_K_0','p_val_K_25', 't_rot', 'UR2']

data_ink = [F_PA_ink_0Umin, F_PA_ink_25kUmin, F_K_ink_0Umin, F_K_ink_25kUmin, x_val_1_pa, x_val_2_pa, p_val_pa_0, p_val_pa_25,
x_val_1_K, x_val_2_K, y_val_1_K, y_val_2_K, p_val_K_0, p_val_K_25, t_rot, UR2]

for i in range(0, len(data_ink),1):
	filename = data_ink_str[i] + odbName +'.txt'
	paramsFile = open (filename,'w')
	data_temp = data_ink[i]
	print(len(data_temp))
	for x in range(0, len(data_temp),1):
		paramsFile.write(str(data_temp[x]) + '\n')
	paramsFile.close()
	

if Ergebnis == 0:
	results = np.zeros(len(parameter))
else:
	results = [F_PA_ges_0Umin, F_PA_ges_25kUmin, F_K_N_ges_0Umin, F_K_N_y_0Umin, F_K_N_x_0Umin, F_K_R_x_0Umin,  F_K_x_ges_0Umin, F_x_ges_0Umin,F_lost_percentage_0Umin,
	F_ratio_0Umin_1, F_ratio_0Umin_2, F_K_N_ges_25kUmin, F_K_N_y_25kUmin, F_K_N_x_25kUmin, F_K_R_x_25kUmin, F_K_x_ges_25kUmin, F_x_ges_25kUmin, F_lost_percentage_25kUmin, F_ratio_25kUmin_1,F_ratio_25kUmin_2, UR2_max]

filename = odbName +'.txt'
paramsFile = open (filename,'w')
for i in range(0, len(parameter),1):
	paramsFile.write(parameter[i]+'\t'+ str(results[i]) + '\n' )
paramsFile.close()

for i in range(0, len(parameter),1):
	filename = parameter[i]+'.txt'
	paramsFile = open (filename,'a')
	paramsFile.write(odbName+'\t'+ str(results[i]) + '\n' )
	paramsFile.close()
