# coding: utf-8
#
# Copyright © 2020 MSC.Software Corporation
#  
# Permission to use Scripts
#  
# Permission is granted to MSC.Software customers with valid licenses of the MSC Apex™ software 
# to use script files posted on this page (collectively “Scripts”) solely for customers’ own authorized 
# internal use in conjunction with licensed use of the MSC Apex software.  
# No other use is permitted.  No license to MSC Apex or other software is granted under this Permission.
#  
# SCRIPTS ARE PROVIDED ON AN "AS-IS" BASIS, WITH NO WARRANTIES OF ANY KIND.  USE OF THE SCRIPTS IS AT CUSTOMER'S SOLE RISK.
#  
# IN NO EVENT SHALL MSC.SOFTWARE CORPORATION, ITS AFFILIATES, SUBSIDIARIES, SUPPLIERS, AGENTS, REPRESENTATIVES OR EMPLOYEES BE LIABLE FOR ANY 
# DIRECT, INDIRECT, INCIDENTAL, SPECIAL OR CONSEQUENTIAL DAMAGES HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, TORT, 
# STRICT LIABILITY, OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION USE OF THE SCRIPTS, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.  
#  
# Customer shall indemnify, defend and hold MSC.Software Corporation, its affiliates, subsidiaries, suppliers, agents, 
# employees or representatives (collectively, “Indemnified Parties”) harmless from and against any third-party claims 
# brought against Indemnified Parties rising out of or in connection with Customer’s use of the Scripts.
#


import apex
from apex.construct import Point3D, Point2D

import math
from operator import itemgetter
import tempfile
import logging
import time


def createCbushProps(fastenerDia, fastenerMatl, elm1, elm2, thick1, thick2, E1, E2, singleShear, length):
    Ef = fastenerMatl.getElasticModulus()

    # a and b for bolted metallic
    # a = 2/3, b = 3

    #a and b for rivet metallic
    # a = 2/5, b = 2.2

    a = 2./3.
    b = 3.0

    # n = 1 for single shear, n = 2 for double shear
    if singleShear:
        n = 1.0
    else:
        n = 2.0

    if length <= 0.0: length = 1.0  # just to make sure i have something.  i will add better error checking later
    s1 = math.pi*fastenerDia**2*Ef/length

    f = ((thick1+thick2)/2/fastenerDia)**a
    f = f + (1./thick1/E1 + 1./n/thick2/E2 + 1./2./thick1/Ef + 1./2./n/thick2/Ef)*b/n

    s2 = 1./f
    s3 = s2

    s4 = 100.0
    s5 = s6 = 1.0E8

    connProps = apex.attribute.createConnectorProperties(
          stiffness               = 0.0,
          damping                 = 0.0,
          translationalStiffnessX = s1,
          translationalStiffnessY = s2,
          translationalStiffnessZ = s3,
          rotationalStiffnessX    = s4,
          rotationalStiffnessY    = s5,
          rotationalStiffnessZ    = s6,
          translationalDampingX   = 0.0,
          translationalDampingY   = 0.0,
          translationalDampingZ   = 0.0,
          rotationalDampingX      = 0.0,
          rotationalDampingY      = 0.0,
          rotationalDampingZ      = 0.0,
          diameter                = 0.0,
          linkMaterial            = None)
    return connProps


def calcDist(n1, n2):
    dx = n2.getX() - n1.getX()
    dy = n2.getY() - n1.getY()
    dz = n2.getZ() - n1.getZ()

    dist = dx**2 + dy**2 + dz**2
    if dist > 0.0: dist = math.sqrt(dist)

    return dist


def setCurrentPart(part):
    myModel = apex.currentModel()

    assy = part.getParent()
    if assy == None:
        partPath = "/"+part.getName()
    else:
        partPath = part.getPath()[part.getPath().find("/")+1:]+"/"+part.getName()

    myModel.setCurrentPart(partPath)


def createRBE3(fastenerName, parentPart, midNode, attachmentNodes):

    #dt = apex.attribute.getDiscreteTies(target=[{"path":parentPart.getPath()+"/"+parentPart.getName()}])
    #if dt.len() == 0:
    #	createDiscreteTie = True
    #else:
    #	createDiscreteTie = False

    # for the moment, I am assuming that the only discrete ties are the ones that I have created
    #if createDiscreteTie:
    #	discreteTie = apex.attribute.createDiscreteTie(
    #             name              = "", 
    #            dependentDoF      = "123456", 
    #              distributionType  = apex.attribute.DistributionType.Compliant,
    #              distributionMode  = apex.attribute.DistributionMode.Custom,
    #              parentInstance    = parentPart)

    path = parentPart.getPath()+"/"+parentPart.getName()
    dofAndWtFactors = {}
    for node in attachmentNodes:
        dist = calcDist(midNode, node)
        dofAndWtFactors[path+":"+node.getId()] = {"123":dist,}

#	_region = apex.EntityCollection()
#	_region.extend(midPoint.getVertices())
    _referenceRegion = midNode
    _attachmentRegion = apex.EntityCollection()
    _attachmentRegion.extend(attachmentNodes)

    nodeTie = apex.attribute.createNodeTie(
              name                 = fastenerName,
              distributionType     = apex.attribute.DistributionType.Compliant,
              distributionMode     = apex.attribute.DistributionMode.Custom,
              dofDefinitionMode    = apex.attribute.DOFDefinitionMode.All,
              referencePoint       = _referenceRegion,
              attachmentRegions    = _attachmentRegion,
              dependentDoF         = "123456",
              DOFAndWeightFactors  = dofAndWtFactors,
              description          = "")
    #print(nodeTie.getName())
    return nodeTie


def calcTriaAreaFromNodes(n1, n2, n3):
    s1 = calcDist(n1, n2)
    s2 = calcDist(n2, n3)
    s3 = calcDist(n3, n1)

    s = 0.5*(s1 + s2 + s3)
    area = math.sqrt(s*(s-s1)*(s-s2)*(s-s3))

    return area


def calcAreaTria(elm):
    elmNodes = elm.getNodes()
    n1 = elmNodes[0]
    n2 = elmNodes[1]
    n3 = elmNodes[2]

    area = calcTriaAreaFromNodes(n1, n2, n3)

    return area


def calcAreaQuad2(elm):
    nodes = elm.getNodes()

    xAvg = yAvg = zAvg = 0.0
    for node in nodes:
        xAvg += node.getX()
        yAvg += node.getY()
        zAvg += node.getZ()
    xAvg = xAvg/4.0
    yAvg = yAvg/4.0
    zAvg = zAvg/4.0
    c = apex.Coordinate(xAvg, yAvg, zAvg)

    area  = calcTriaAreaFromNodes(c, nodes[0], nodes[1])
    area += calcTriaAreaFromNodes(c, nodes[1], nodes[2])
    area += calcTriaAreaFromNodes(c, nodes[2], nodes[3])
    area += calcTriaAreaFromNodes(c, nodes[3], nodes[0])

    return area


def calcAreaQuad(elm):
    elmNodes = elm.getNodes()
    n1 = elmNodes[0]
    n2 = elmNodes[1]
    n3 = elmNodes[2]
    n4 = elmNodes[3]

    # I am calculating the area of triangles in the quad
    # I am averaging area from 2 possible dissections of the quad
    # I will try to make this better soon
    a1 = calcTriaAreaFromNodes(n1, n2, n3)
    a2 = calcTriaAreaFromNodes(n1, n3, n4)
    area = a1+a2

    a1 = calcTriaAreaFromNodes(n1, n2, n4)
    a2 = calcTriaAreaFromNodes(n2, n3, n4)
    area = 0.5*(area + (a1+a2))

    return area


def calcElmArea(elm):
    elmTopology = elm.getTopology()
    if elmTopology == apex.mesh.ElementTopology.Tria3:
        area = calcAreaTria(elm)
    if elmTopology == apex.mesh.ElementTopology.Quad4:
        area = calcAreaQuad(elm)

    return area


def calcNodeArea(node, elm):
    nodes = elm.getNodes()
    elmTopology = elm.getTopology()
    area = 0.0

    if elmTopology == apex.mesh.ElementTopology.Tria3:
        area  = calcTriaAreaFromNodes(node, nodes[0], nodes[1])
        area += calcTriaAreaFromNodes(node, nodes[1], nodes[2])
        area += calcTriaAreaFromNodes(node, nodes[2], nodes[0])
    else:
        a1 = calcTriaAreaFromNodes(node, nodes[0], nodes[1])
        a2 = calcTriaAreaFromNodes(node, nodes[1], nodes[2])
        a3 = calcTriaAreaFromNodes(node, nodes[2], nodes[3])
        a4 = calcTriaAreaFromNodes(node, nodes[3], nodes[0])

        area = a1+a2+a3+a4

    return area


def getPiercedElm(nearNode, connectionPoint):
    elms = nearNode.getElements()

    elmData = []
    for elm in elms:
        baseArea = calcElmArea(elm)

        nodeArea = calcNodeArea(connectionPoint, elm)

        elmData.append([elm, abs(baseArea-nodeArea)])

    sortedData = sorted(elmData, key=itemgetter(1), reverse=False) # sort on area difference (ascending order)

    return sortedData[0][0]


def main(dict={}):
    startTime = time.process_time()

    apex.disableUndoRedo()
    apex.disableRecoveryFileRecording()

    myModel = apex.currentModel()

    scriptName = "auSpiderFast.py"
    logFileName = myModel.getName()+"_spiderFast.log"

    try:
        tmp = tempfile.gettempdir()
        logfile=tmp+"\\"+logFileName
    except:
        logfile="C:\\"+logFileName

    try:
        logging.basicConfig(filename=logfile, filemode='w', level=logging.DEBUG)
        print("Output data being written to: "+logfile)
        logEnabled = True
    except:
        logEnabled = False

    if logEnabled:
        logging.info("Script: "+scriptName)
        logging.info("Send questions/comments to larry.pearce@mscsoftware.com")
        logging.info("#####################################################")

        logging.info("Apex model: "+myModel.getName())
        logging.info(time.strftime("%a, %d %b %Y %H:%M:%S", time.localtime()))

    main2(dict, logEnabled)

    endTime = time.process_time()

    seconds = endTime-startTime
    min, sec = divmod(seconds, 60)

    msg = "Approximate elapsed time: "+"%02d : %02d" % (min, sec)

    logging.info(msg)
    print("### Normal exit")
    print(msg)

    if logEnabled: logging.shutdown()

    apex.enableUndoRedo()
    apex.enableRecoveryFileRecording()


def myEcho(msg, logEnabled):
    print(msg)
    if logEnabled: logging.info(msg)


def main2(dict, logEnabled):
    apex.disableShowOutput()
    #apex.setScriptUnitSystem(unitSystemName = r'''in-slinch-s-lbf''')
    myModel = apex.currentModel()

    matlName     = dict["material"]
    fastenerDia  = float(dict["diameter"])
    fastenerName = dict["name"]

    if fastenerName == "": fastenerName = "Connector"

    myEcho("Fastener Diameter:"  + str(fastenerDia), logEnabled)
    myEcho("Fastener Material:"  + str(matlName),    logEnabled)
    myEcho("Name: "              + fastenerName,     logEnabled)

    if fastenerDia <= 0.0:
        myEcho("Fastener diameter is less than zero!", logEnabled)
        return
    fastenerRadius = 0.5*fastenerDia

    myMatl = None
    matls = apex.catalog.getMaterials([{"path":myModel.getName()+"/Materials"}])
    for matl in matls:
        if matl.getName() == matlName:
            myMatl = matl
            break
    if myMatl == None:
        print("Unknown material?")
        return

    parts = myModel.getParts(True)

    _selected = apex.selection.getCurrentSelection()
    curves = apex.EntityCollection()
    for sel in _selected:
        if sel.getEntityType() == apex.EntityType.Curve:
            curves.append(sel)

    # if no curves are selected, I can get all curves in the db.  perhaps this is not a good idea
    if curves.len() == 0:
        for part in parts:
            curves += part.getCurves()

    if curves.len() == 0:
        myEcho("No curves in the database!", logEnabled)
        return

    surfs = apex.EntityCollection()
    for part in parts: surfs += part.getSurfaces()
    for part in parts: surfs += part.getFacetedSurfaces()

    myEcho("Parts to be processed: "+str(parts.len()),    logEnabled)
    myEcho("Surfaces to be processed: "+str(surfs.len()), logEnabled)
    myEcho("Curves to be processed: "+str(curves.len()),  logEnabled)

    if surfs.len() < 2:
        myEcho("Less than 2 Surfaces to be connected!", logEnabled)
        return

    numPointsCreated  = 0
    numSpidersCreated = 0

    curveData = {} # Dictionary of lists of data for each fastener curve
    intCurves = [] # list of curves that intersect surfaces
    for curve in curves:
        edges = curve.getEdges()
        for surf in surfs:
            elms = surf.getElements()
            if elms.len() == 0: continue # if the surface has no elements then skip the surface

            setCurrentPart(surf.getParent())
            faces = surf.getFaces()

            # Create points at the intersection of the bolt curves and the surface faces
            points = apex.geometry.createPointCurveSurface(targetCurves = edges, intersectingSurfaces = faces)
            if points.len() != 1:
                if points.len() > 0: apex.deleteEntities(points)
                continue

            # At this point, the curve does intersect the surface
            connectionPoint = points[0]
            numPointsCreated += 1

            cData = {"surf":surf, "connectionPoint":connectionPoint, "rbeCreated":False, "midNode":None, "nodes":None, "elm":None, "singleShear":True, "skipCurve":False, "thick":-1.0, "E":0.0, "length":0.0}
            if curve.getParent().getName()+"-"+curve.getName() in curveData:
                curveData[curve.getParent().getName()+"-"+curve.getName()].append(cData)
            else:
                curveData[curve.getParent().getName()+"-"+curve.getName()] = [cData]
                intCurves.append(curve)

    numFasteners = numApexConnectors = 0

    # at this point i have a point at all curve-surface intersections

    # find the pierced element
    for curve in intCurves:
        for data in curveData[curve.getParent().getName()+"-"+curve.getName()]:
            connectionPoint = data["connectionPoint"]
            nodes = data["surf"].getNodes()
            if nodes.len() == 0: # The surface is apparently not meshed
                data["skipCurve"] = True
                break
            mySearch = apex.utility.ProximitySearch()
            mySearch.insertCollection(nodes)
            zNearbyObject = mySearch.findNearestObject(connectionPoint.getLocation())
            nearbyObject = zNearbyObject.getNearestObject()

            nearNode   = nearbyObject.asNode()
            piercedElm = getPiercedElm(nearNode, connectionPoint)

            data["nodes"] = piercedElm.getNodes()
            data["elm"]   = piercedElm
            data["thick"] = piercedElm.getSection().getThickness()
            # this is to work around an api issue
            try:
                data["E"] = piercedElm.getMaterial().getElasticModulus()
            except:
                data["E"] = 68947.5 # MPa for aluminum


    for curve in intCurves:
        if curve.getParent().getName()+"-"+curve.getName() not in curveData: continue # Then this curve has no surface intersections
        # i dont like this
        skipFlag = False
        for data in curveData[curve.getParent().getName()+"-"+curve.getName()]:
            if data["skipCurve"] == True:
                break
                skipFlag = True
        if skipFlag:  continue
        # plz improve

        numConnections = len(curveData[curve.getParent().getName()+"-"+curve.getName()]) - 1 # len(curveData) must be at least 2, i.e., the curve must intersect 2 geometry bodies for 1 connection
        if numConnections <= 0: # The curve intersects only 1 or possibly zero geometry bodies
            #myEcho("No connection will be created for curve: "+curve.getName(), logEnabled)
            continue

        if numConnections > 1: # The curve intersects more than 2 geometry bodies
            # If the curve intersects 3 or more geometry bodies, then I need to get the correct connection order
            sortData = []
            k = 0
            for data in curveData[curve.getParent().getName()+"-"+curve.getName()]:
                u = curve.evaluatePointOnCurve(data["connectionPoint"]).getParametricCoordinates()
                sortData.append([u, data["connectionPoint"], k])
                k += 1
                data["singleShear"] = False
            locData = sorted(sortData, key=itemgetter(0), reverse=False) # ascending sort on u

        else: # The curve intersects exactly 2 geometry bodies
            locData = [[0.0, curveData[curve.getParent().getName()+"-"+curve.getName()][0]["connectionPoint"], 0], [1.0, curveData[curve.getParent().getName()+"-"+curve.getName()][1]["connectionPoint"], 1]]
            for data in curveData[curve.getParent().getName()+"-"+curve.getName()]:
                data["singleShear"] = True

        # At this point I am getting ready to create the fasteners for the curve
        #nodeTies = []
        for data in curveData[curve.getParent().getName()+"-"+curve.getName()]:
            setCurrentPart(curve.getParent())
            midNode = apex.mesh.createNodeByPickLocation( target = data["connectionPoint"].getVertices()[0],
                          location = data["connectionPoint"].getLocation())
            data["midNode"] = midNode

        # calculate fastener length
        n1 = curveData[curve.getParent().getName()+"-"+curve.getName()][0]["midNode"]
        n2 = curveData[curve.getParent().getName()+"-"+curve.getName()][-1]["midNode"]
        dist = calcDist(n1, n2)
        dist = 0.5*curveData[curve.getParent().getName()+"-"+curve.getName()][0]["thick"] + 0.5*curveData[curve.getParent().getName()+"-"+curve.getName()][-1]["thick"]
        for data in curveData[curve.getParent().getName()+"-"+curve.getName()]:
            data["length"] = dist

        connections = []
        for k in range(numConnections):
            connections.append([locData[k][2], locData[k+1][2]])
        connStr = ""
        for connect in connections:
            d1 = curveData[curve.getParent().getName()+"-"+curve.getName()][connect[0]]
            d2 = curveData[curve.getParent().getName()+"-"+curve.getName()][connect[1]]

            _endPoint1 = apex.construct.createLocationByEntity(entity = d1["midNode"])
            _endPoint2 = apex.construct.createLocationByEntity(entity = d2["midNode"])

            cbushProps = createCbushProps(fastenerDia, myMatl, d1["elm"], d2["elm"], d1["thick"], d2["thick"], d1["E"], d2["E"], d1["singleShear"], d1["length"])

            myConnector = apex.attribute.createConnector(
                          name = fastenerName,
                          connectorType = apex.attribute.ConnectorType.Bushing,
                          connectorProperties = cbushProps,
                          applicationMethod1 = apex.attribute.ApplicationMethod.Direct,
                          applicationMethod2 = apex.attribute.ApplicationMethod.Direct,
                          end1InterfacePoint = _endPoint1,
                          end2InterfacePoint = _endPoint2,
                          description = "")
            numApexConnectors += 1
            if connStr == "":
                connStr = myConnector.getName()
            else:
                connStr = connStr + ":" + myConnector.getName()

            name = myConnector.getName()

            if d1["rbeCreated"] == False:
                nodeTieName = name+"-d1"
                if d1["surf"]  != None: parentPart = d1["surf"].getParent()
                nodeTie = createRBE3(nodeTieName, parentPart, d1["midNode"], d1["nodes"])
                d1["rbeCreated"]=True
                numSpidersCreated += 1
            if d2["rbeCreated"] == False:
                nodeTieName = name+"-d2"
                if d2["surf"]  != None: parentPart = d2["surf"].getParent()
                nodeTie = createRBE3(nodeTieName, parentPart, d2["midNode"], d2["nodes"])
                d2["rbeCreated"]=True
                numSpidersCreated += 1
        numFasteners += 1

    myEcho("Number of Points created: "+str(numPointsCreated), logEnabled)
    myEcho("Number of RBE3 spiders created:"+str(numSpidersCreated), logEnabled)
    myEcho(" ", logEnabled)
    myEcho("Number of Fasteners created: "+str(numFasteners), logEnabled)
    myEcho("Number of Apex Connections created: "+str(numApexConnectors), logEnabled)



