# Note:
#  The deformed coordinates of the nodes are required
#  Request COORD as ODB Field Output in A/CAE
#  Use *Node Output when adding the output request into the .inp
#  Provide a instance name and optional a set name for a subset of the instance
###################################################################################################################


#### User Input

# Name of instance
instanceName='PART-1-1'

#Assign elementSet='' if the element normal for all the elements in the instance has to be computed
elementSet=''

# Output file
outputFile='ElementNormal.txt'

####End of User Inputs

###################################################################################################################

from abaqus import *
from abaqusConstants import *


vps = session.viewports[session.currentViewportName]
vpname = vps.name

odbName = vps.displayedObject.name
odb = session.odbs[odbName]
currentstepnumber = vps.odbDisplay.fieldFrame[0]
currentstepname = odb.steps.keys()[currentstepnumber]
currentframe = vps.odbDisplay.fieldFrame[1]




#######################################################
#Method to find the cross product of two vectors

def crossProduct(a, b):
    c = [a[1]*b[2] - a[2]*b[1],a[2]*b[0] - a[0]*b[2],a[0]*b[1] - a[1]*b[0]]
    return c

###################################################################################################################

try:
    o=session.odbs[odbName]
except KeyError:
    raise KeyError, ('The odb should be displayed in the viewport: '+odbName)

#Check if the instance exists
try:
    i=o.rootAssembly.instances[instanceName]
except KeyError:
    raise KeyError, ('The instance %s does not exist.'%instanceName)

#Check if the Element set exists
try:
    if elementSet != '':
        e=i.elementSets[elementSet].elements
    else:
        e=i.elements
except KeyError:
    raise KeyError, ('The element set %s does not exist.'%elementSet)

#Check if there are elements in the element set
if not e:
    raise ValueError, ('There are no elements in the element set %s. Please check the element set contents'%elementSet)

#n=i.nodes


###################################################################################################################


f=open(outputFile,'w')
print '\n************Unit Normal of Elements in Set: '+instanceName+' '+elementSet+'**************************\n'
f.write('************Unit Normal of Elements in Set: '+instanceName+' '+elementSet+'**************************\n')

coords = odb.steps[currentstepname].frames[currentframe].fieldOutputs['COORD']
for elem in e:
    c=elem.connectivity
    if len(c)>2:
    	
        #Get deformed coordinates
        node0=i.getNodeFromLabel(c[0])
        x = coords.getSubset(region=node0)
        n0=x.values[0].data

        node1=i.getNodeFromLabel(c[1])
        x = coords.getSubset(region=node1)
        n1=x.values[0].data
 
        node2=i.getNodeFromLabel(c[2])
        x = coords.getSubset(region=node2)
        n2=x.values[0].data
 

        #Contruct a vector between the first two nodes
        v1=[n1[0]-n0[0],n1[1]-n0[1],n1[2]-n0[2]]
        #Contruct a vector between the second & third node
        v2=[n2[0]-n1[0],n2[1]-n1[1],n2[2]-n1[2]]
        
        #Calculate the normal and unit normal vector
        normal=crossProduct(v1,v2)
        unitNormal=normal/(normal[0]**2+normal[1]**2+normal[2]**2)**0.5
        
        print "Element: " + str(elem.label) + "\t ",  unitNormal
        f.write("Element: " + str(elem.label) + "\t "+str(unitNormal)+'\n')

f.close()
print '\nOutput also written to file:', outputFile