# 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
import time


def getFaceForEdge(selEdge):
    # this is a work-around because edge.getReferncedBy() is not working in IL-FP1

    tol = .001

    targetVerts    = selEdge.getVertices()
    targetMidPoint = selEdge.getMidPoint()
    targetLength   = selEdge.getLength()

    myModel = apex.currentModel()
    parts = myModel.getParts(True)

    allFaces = apex.entityCollection()
    for part in parts:
        surfs = part.getSurfaces()
        for surf in surfs:
            allFaces += surf.getFaces()

    mySearch = apex.utility.ProximitySearch()
    mySearch.insertCollection(targetVerts)
    for face in allFaces:
        edges = face.getEdges()
        for edge in edges:
            if abs(edge.getLength() - targetLength) > tol: continue
            mp = edge.getMidPoint()
            dist = calcDist(mp, targetMidPoint)
            if dist > tol: continue
            verts = edge.getVertices()
            numVertsMatch = 0
            for vert in verts:
                nearObj = mySearch.findNearestObject(vert)
                dist = nearObj.getDistance()
                if dist > tol: break
                numVertsMatch += 1
            if numVertsMatch == 2: return face

    return None


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 processInputForce(f_in):
    if f_in == "":
        f = 0.0
    else:
        try:
            f = float(f_in)
        except:
            f = 0.0
    return f


def isEdgeStraight(edge):
    verts = edge.getVertices()

    mp = edge.getMidPoint()

    v1x = verts[0].getX() - mp.getX()
    v1y = verts[0].getY() - mp.getY()
    v1z = verts[0].getZ() - mp.getZ()
    dot = math.sqrt(v1x**2 + v1y**2 + v1z**2)
    v1x = v1x/dot
    v1y = v1y/dot
    v1z = v1z/dot

    v2x = verts[1].getX() - mp.getX()
    v2y = verts[1].getY() - mp.getY()
    v2z = verts[1].getZ() - mp.getZ()
    dot = math.sqrt(v2x**2 + v2y**2 + v2z**2)
    v2x = v2x/dot
    v2y = v2y/dot
    v2z = v2z/dot

    cosTheta = v1x*v2x + v1y*v2y + v1z*v2z
    cosTheta = abs(cosTheta)

    if (1.0 - cosTheta) <= .01:
        return True
    else:
        return False


def createForce(lbcName, forceMomentRep_, _target, _loadEventOrientation):
    forcemoment = apex.environment.createForceMoment(
       name = lbcName,
       forceMomentRep = forceMomentRep_,
       applicationMethod = apex.attribute.ApplicationMethod.Direct,
       target = _target,
       orientation = _loadEventOrientation,
       description = "")


def crossProduct(v1, v2):
    a1 = v1.getX()
    a2 = v1.getY()
    a3 = v1.getZ()

    b1 = v2.getX()
    b2 = v2.getY()
    b3 = v2.getZ()

    s1 = a2*b3 - a3*b2
    s2 = a3*b1 - a1*b3
    s3 = a1*b2 - a2*b1

    v = apex.construct.Vector3D(x = s1, y = s2, z = s3)

    return v


def makeUnitVector(aVec):
    d = aVec.getX()**2 + aVec.getY()**2 + aVec.getZ()**2
    if d == 0:
        vec = apex.construct.Vector3D(x = 0., y = 0., z = 0.)
    else:
        d = math.sqrt(d)
        dx = aVec.getX()/d
        dy = aVec.getY()/d
        dz = aVec.getZ()/d

        vec = apex.construct.Vector3D(x = dx, y = dy, z = dz)

    return vec


def main(dict={}):

    unitSystems = {}
    unitSystems["inch"] = r'''in-slinch-s-lbf-deg'''
    unitSystems["feet"] = r'''ft-slug-s-lbf-deg'''
    unitSystems["mm"]   = r'''mm-t-s-N-deg'''
    unitSystems["m"]    = r'''m-kg-s-N-deg'''

    apex.setScriptUnitSystem(unitSystemName = unitSystems[dict["unit"]])

    print("Selected base unit: "+dict["unit"])

    myModel = apex.currentModel()

    lbcName = dict["name"]
    if lbcName == "": lbcName = "Distr Edge Load"
    f1Val = dict["f1"]
    f2Val = dict["f2"]
    f3Val = dict["f3"]

    f1 = processInputForce(f1Val)
    f2 = processInputForce(f2Val)
    f3 = processInputForce(f3Val)

    if f1+f2+f3 == 0.0:
        apex.session.displayStatusMessage("No Load values entered?")
        time.sleep(5)
        print("No Force values entered?")
        return

    _selected = apex.selection.getCurrentSelection()

    if _selected.len() == 0:
        apex.session.displayStatusMessage("Nothing selected?")
        time.sleep(5)
        print("Nothing selected?")
        return

    totalDist = 0.0
    for sel in _selected:
        if sel.getEntityType() == apex.EntityType.Edge:
            totalDist += sel.getLength()

    for sel in _selected:
        if sel.getEntityType() == apex.EntityType.Edge:
            selNodes = sel.getNodes()
        else:
            continue

        if isEdgeStraight(sel) == False:
            apex.session.displayStatusMessage("Selected edge is not straight?")
            time.sleep(5)
            print("Selected edge is not straight?")
            return

        print("Selected nodes: "+str(selNodes.len()))
        if selNodes.len() == 0:
            apex.session.displayStatusMessage("No nodes were selected?")
            time.sleep(5)
            print("No nodes were selected?")
            return

        selElms = apex.entityCollection()
        for node in selNodes:
            selElms += node.getElements()

        nodeIds = []
        nodes = apex.entityCollection()
        for node in selNodes:
            if node.getId() in nodeIds: continue
            nodeIds.append(node.getId())
            nodes.append(node)

        elmIds = []
        elms = apex.entityCollection()
        for elm in selElms:
            if elm.getId() in elmIds: continue
            elmIds.append(elm.getId())
            elms.append(elm)

        nodeLoads = {}
        for node in nodes:
            nodeLoads[node.getId()] = [node, 0, 0.0, 0.0, 0.0]

        totalLoad = [0.0, 0.0, 0.0]
        for elm in elms:  # was selElms
            elmNodes = elm.getNodes()
            elmTop   = elm.getTopology()

            if elmTop == apex.mesh.ElementTopology.Quad4:
                edges = [[0, 1], [1, 2], [2, 3], [3, 0]]
            elif elmTop == apex.mesh.ElementTopology.Tria3:
                edges = [[0, 1], [1, 2], [2, 0]]
            else:
                continue

            for edge in edges:
                eNodes = []
                ids = []
                ids.append(elmNodes[edge[0]].getId())
                ids.append(elmNodes[edge[1]].getId())
                if ids[0] in nodeIds and ids[1] in nodeIds:
                    eNodes.append(elmNodes[edge[0]])
                    eNodes.append(elmNodes[edge[1]])
                    break
            if len(eNodes) == 0: continue

            dist = calcDist(eNodes[0], eNodes[1])

            ff1 = f1*dist/totalDist
            ff2 = f2*dist/totalDist
            ff3 = f3*dist/totalDist

            totalLoad[0] += ff1
            totalLoad[1] += ff2
            totalLoad[2] += ff3

            print("f1 load for elm: "+str(elm.getId())+" = "+str(ff1))
            print("f2 load for elm: "+str(elm.getId())+" = "+str(ff2))
            print("f3 load for elm: "+str(elm.getId())+" = "+str(ff3))

            nData = nodeLoads[eNodes[0].getId()]
            nData[1] += 1
            nData[2] += 0.5*ff1
            nData[3] += 0.5*ff2
            nData[4] += 0.5*ff3

            nData = nodeLoads[eNodes[1].getId()]
            nData[1] += 1
            nData[2] += 0.5*ff1
            nData[3] += 0.5*ff2
            nData[4] += 0.5*ff3

        verts = sel.getVertices()

        dx = verts[1].getX() - verts[0].getX()
        dy = verts[1].getY() - verts[0].getY()
        dz = verts[1].getZ() - verts[0].getZ()
        dist = math.sqrt(dx**2 + dy**2 + dz**2)
        dx = dx/dist
        dy = dy/dist
        dz = dz/dist

        yVec = apex.construct.Vector3D(x = dx, y=dy, z=dz)  # was yVec

        face = getFaceForEdge(sel)  # edge.getReferencedBy() is not returning anything in IL-FP1
        if face == None:
            print("No face?")
            return

        loc = face.evaluatePointOnFace(verts[0])
        zVec = face.evaluateNormal(loc.getU(), loc.getV())  # was zVec

        xVec = crossProduct(yVec, zVec)

        d = zVec.getX()**2 + zVec.getY()**2
        if d == 0.0:
            _aa = 0.0
            if zVec.getZ() > 0.0:
                _bb = 0.0
            else:
                _bb = 0.0
            _gg = -1.0*math.atan2(xVec.getY(), xVec.getX())
        else:
            d = math.sqrt(d)

            zGlobal = apex.construct.Vector3D(x=0.0, y=0.0, z=1.0)
            xp = crossProduct(zVec, zGlobal)

            x1px = xp.getX()*xVec.getX() + xp.getY()*xVec.getY() + xp.getZ()*xVec.getZ()
            x1py = xp.getX()*yVec.getX() + xp.getY()*yVec.getY() + xp.getZ()*yVec.getZ()

            _aa = math.atan2(x1py, x1px)
            _bb = math.atan2(d, zVec.getZ())
            #_gg = -1.0*math.atan2(xVec.getY(), xVec.getX())
            _gg = -1.0*math.atan2(xp.getY(), xp.getX())

        _aa = _aa*180./math.pi
        _bb = _bb*180./math.pi
        _gg = _gg*180./math.pi

        _aa = -1.0*_aa
        _bb = -1.0*_bb
        _gg = -1.0*_gg

        _alpha = _gg
        _beta  = _bb
        _gamma = _aa

        print("alpha: "+str(_alpha))
        print("beta:  "+str(_beta))
        print("gamma: "+str(_gamma))


        _orientation = apex.construct.createOrientation(alpha=_alpha, beta=_beta, gamma=_gamma)
        #_orientation.setXAxis(xVec)
        #_orientation.setYAxis(yVec)
        #_orientation.setZAxis(zVec)

        coordSys = apex.construct.createCoordinateSystemByEulerMethod(
            name = "Coordinate System",
            description = "",
            coordType = apex.construct.CoordinateSystemType.Cartesian,
            origin = verts[0],
            orientation=_orientation)

        _loadEventOrientation = apex.construct.createOrientationByCoordinateSystem(coordSys)

        # this is not precise but is the best i can come up with for the moment
        # this is to determine the sense of the edge
        mp = face.getCentroid()
        dx = mp.getX() - verts[0].getX()
        dy = mp.getY() - verts[0].getY()
        dz = mp.getZ() - verts[0].getZ()
        mpVec = apex.construct.Vector3D(x = dx, y = dy, z = dz)
        mpLength = mpVec.getLength()

        testLength = 0.5*mpLength

        xv = makeUnitVector(xVec)

        testLocX = verts[0].getX() + xv.getX()*testLength
        testLocY = verts[0].getY() + xv.getY()*testLength
        testLocZ = verts[0].getZ() + xv.getZ()*testLength

        dx = mp.getX() - testLocX
        dy = mp.getY() - testLocY
        dz = mp.getZ() - testLocZ
        newVec = apex.construct.Vector3D(x = dx, y = dy, z = dz)
        dist = newVec.getLength()

        if dist > mpLength:
            f = 1.0
        else:
            f = -1.0
        # end

        if nodes.len() > 0:
            for node in nodes:
                if nodeLoads[node.getId()][1] == 0: continue
                fx = f*nodeLoads[node.getId()][2]
                fy = nodeLoads[node.getId()][3]
                fz = nodeLoads[node.getId()][4]

                forceMomentRep_ = apex.environment.createForceMomentStaticRepByComponent(
                   name   = lbcName,
                   forceX = fx,
                   forceY = fy,
                   forceZ = fz,
                   description = "")
                _target = apex.entityCollection()
                _target.append(node)
                createForce(lbcName, forceMomentRep_, _target, _loadEventOrientation)

    print("summation of f1: "+str(totalLoad[0]))
    print("summation of f2: "+str(totalLoad[1]))
    print("summation of f3: "+str(totalLoad[2]))



