# 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.
#

# Author: Larry Pearce
# Revised by: Larry Pearce
# Revision History
# Date         Change                                 By
# 2020-10-13   Original code (1.1)                    Larry Pearce

import apex

import math


def getFaceData(face):
    # gets the area and center or centroid/center location
    nodes = face.getNodes()

    if nodes.len() == 3 or nodes.len() == 6:
        data = triaData(nodes[0:3])
    else:
        data = quadData(nodes[0:4])

    return data


def quadData(nodes):
    # gets the area and centroid/center for quad elements
    nn = [nodes[0], nodes[1], nodes[3]]
    data = triaData(nn)
    area = data[0]
    nn = [nodes[1], nodes[2], nodes[3]]
    data = triaData(nn)
    area += data[0]

    #
    v12 = createVector3D(nodes[0], nodes[1])
    v13 = createVector3D(nodes[0], nodes[2])

    v = crossProduct(v12, v13)

    vNorm   = v.getLength()
    v12Norm = v12.getLength()
    v13Norm = v13.getLength()

    sin_beta = vNorm/v12Norm/v13Norm
    # this is carry over from PCL, not sure it is needed here
    if sin_beta > 1.0:  sin_beta = 1.0
    if sin_beta < -1.0: sin_beta = -1.0

    beta = math.asin(sin_beta) # beta is the angle between v12 and v13
    # use cross product to determine angle (gamma) between v21 and v24
    v21 = createVector3D(nodes[1], nodes[0])
    v24 = createVector3D(nodes[1], nodes[3])

    v = crossProduct(v24, v21)

    vNorm = v.getLength()
    v21Norm = v21.getLength()
    v24Norm = v24.getLength()

    sin_gamma = vNorm/v21Norm/v24Norm
    # this is carry over from PCL, not sure it is needed here
    if sin_gamma > 1.0:  sin_gamma = 1.0
    if sin_gamma < -1.0: sin_gamma = -1.0

    gamma = math.asin(sin_gamma) # angle between v21 and v24

    # calculate phi from beta and gamma
    phi = math.pi - (beta + gamma) # angle between v13 and v24 or v31 and v42
    sin_phi = math.sin(phi)

    # g_dist is opposite the angle gamma
    # law of sines
    g_dist = v12Norm*sin_gamma/sin_phi

    v13 = makeUnitVector(v13)
    cx = nodes[0].getX() + g_dist*v13.getX()
    cy = nodes[0].getY() + g_dist*v13.getY()
    cz = nodes[0].getZ() + g_dist*v13.getZ()

    centerPoint = apex.Coordinate(cx, cy, cz)

    return [area, centerPoint]


def triaData(nodes):
    # gets the area and centroid/center for tria elements
    v1 = createVector3D(nodes[0], nodes[1])
    v2 = createVector3D(nodes[0], nodes[2])

    vv = crossProduct(v1, v2)

    area = 0.5*vv.getLength()

    cx = cy = cz = 0.0
    for node in nodes:
        cx += node.getX()
        cy += node.getY()
        cz += node.getZ()

    cx = cx/3.0
    cy = cy/3.0
    cz = cz/3.0

    centerPoint = apex.Coordinate(cx, cy, cz)


    return [area, centerPoint]


def createVector3D(p1 , p2):
    dx = p2.getX() - p1.getX()
    dy = p2.getY() - p1.getY()
    dz = p2.getZ() - p1.getZ()

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


def getRadForEdges(edges):

    edgeRadius = []
    for edge in edges:
        dia = apex.measureDiameter(edge)
        radius = 0.5*dia
        edgeRadius.append(radius)

    # uses a formula to calculate radius based on triangle edge lengths and area
    # Radius = lenA * lenB * lenC /4/area
    #edgeRadius = []
    #for edge in edges:
    #	mp    = edge.getMidPoint()
    #	verts = edge.getVertices()
    #
    #	vc = createVector3D(verts[0], verts[1])
    #
    #	va = createVector3D(verts[0], mp)
    #
    #	vb = createVector3D(verts[1], mp)
    #
    #	vv = crossProduct(vc, va)
    #
    #	area = 0.5*vv.getLength()
    #
    #	edgeRadius.append(va.getLength()*vb.getLength()*vc.getLength()/4.0/area)

    return edgeRadius


def getEdgeVector(edge):
    # creates a vector based on the positions of the end vertices of an edge
    verts = edge.getVertices()
    if verts.len() != 2: return None

    edgeVec = createVector3D(verts[0], verts[1])

    edgeVec = makeUnitVector(edgeVec)

    return edgeVec


def areEdgesParallel(edges):
    # used dot product to check to see if 2 edges are parallel
    edgeVectors = []
    for edge in edges:
        edgeVec = getEdgeVector(edge)
        if edgeVec is None: return False
        edgeVectors.append(edgeVec)

    dot = 0.0
    dot += edgeVectors[0].getX()*edgeVectors[1].getX()
    dot += edgeVectors[0].getY()*edgeVectors[1].getY()
    dot += edgeVectors[0].getZ()*edgeVectors[1].getZ()

    if (abs(dot) - 1.0) <= .0001: return True

    return False


def isEdgeStraight(edge):
    # checks to see if edge is straight by calculating the dot product of 2 vectors from midpoint to end vertices
    verts = edge.getVertices()

    mp = edge.getMidPoint()

    v1 = createVector3D(mp, verts[0])
    v2 = createVector3D(mp, verts[1])
    cosTheta = dotProduct(v1, v2)
    cosTheta = abs(cosTheta)

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


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 dotProduct(Va, Vb):
    dot  = Va.getX()*Vb.getX()
    dot += Va.getY()*Vb.getY()
    dot += Va.getZ()*Vb.getZ()

    cosTheta = dot/Va.getLength()/Vb.getLength()

    return cosTheta


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 getElmFacesForFace(gFace):

    elmTopo = {}
    elmTopo[apex.mesh.ElementTopology.Bar2]    = False
    elmTopo[apex.mesh.ElementTopology.Tria3]   = False
    elmTopo[apex.mesh.ElementTopology.Tria6]   = False
    elmTopo[apex.mesh.ElementTopology.Quad4]   = False
    elmTopo[apex.mesh.ElementTopology.Quad8]   = False
    elmTopo[apex.mesh.ElementTopology.Tetra4]  = True
    elmTopo[apex.mesh.ElementTopology.Tetra10] = True
    elmTopo[apex.mesh.ElementTopology.Penta6]  = True
    elmTopo[apex.mesh.ElementTopology.Penta15] = True
    elmTopo[apex.mesh.ElementTopology.Hexa8]   = True
    elmTopo[apex.mesh.ElementTopology.Hexa20]  = True

    faceNodes = apex.entityCollection()
    faceNodes += gFace.getExteriorNodes()
    faceNodes += gFace.getInteriorNodes()

    if faceNodes.len() == 0: return None

    faceNodeIds = []
    for node in faceNodes:
        if node.getId() in faceNodeIds: continue
        faceNodeIds.append(node.getId())

    nodeElms   = apex.entityCollection()
    nodeElmIds = []
    for node in faceNodes:
        elms = node.getElements()
        for elm in elms:
            if elm.getId() in nodeElmIds: continue
            if not elmTopo[elm.getTopology()]: continue
            nodeElms.append(elm)
            nodeElmIds.append(elm.getId())

    elmFaces = []
    for elm in nodeElms:
        faces = elm.getFaces()
        for face in faces:
            nodes = face.getNodes()
            if nodes.len() == 3 or nodes.len() == 6: limit = 3
            if nodes.len() == 4 or nodes.len() == 8: limit = 4
            foundFace = True
            for node in nodes[0:limit]:  # I am only checking corner nodes
                if node.getId() in faceNodeIds: continue
                foundFace = False
                break
            if foundFace: elmFaces.append(face)

    return elmFaces


def main(dict={}):
#def main():
    unitSystem = r'''in-slinch-s-lbf'''
    print("Using unit system: "+unitSystem)
    myModel = apex.currentModel()

    load = dict["load"]
    if load == "":
        print("Load is BLANK!")
        return
    load = float(load)

    lbcName = dict["name"]
    if lbcName == "": lbcName = "Pressure"

    distrType = dict["distrType"]

    areaAdjustment = False
    if dict["areaAdjust"].lower() == "true":
        areaAdjustment = True

    print("Load: "+ str(load))
    print("LBC Name:" +lbcName)
    print("Distribution Type: "+distrType)
    print("Geom/FE area adjustment: "+str(areaAdjustment))


    _selected = apex.selection.getCurrentSelection()

    selFaces = apex.entityCollection()
    for sel in _selected:
        if sel.getEntityType() == apex.EntityType.Face: selFaces.append(sel)
    if selFaces.len() == 0:
        print("No faces selected!")
        return
    if selFaces.len() > 1:
        print("Only 1 face may be selected!  Processing the first face...")

    selFace = selFaces[0]

    edges = selFace.getEdges()
    if edges.len() != 4:
        print("Selected face must have only 4 edges!  The selected face has "+str(edges.len())+" edges ...")
        return

    straightEdges = []
    curvedEdges   = []
    for edge in edges:
        if isEdgeStraight(edge):
            straightEdges.append(edge)
        else:
            curvedEdges.append(edge)
    if len(straightEdges) != 2:
        print("Expecting 2 straight edges! "+str(straightEdges.len())+" straight edges found ...")
        return

    if not areEdgesParallel(straightEdges):
        print("Straight edges are not parallel!")
        return

    if abs(straightEdges[1].getLength() - straightEdges[0].getLength()) > .001:
        print("Different lengths for 2 straight edges!")
        print("Length 1: "+straightEdges[0].getLength())
        print("Length 2: "+straightEdges[1].getLength())
        return

    thick = straightEdges[0].getLength()

    edgeRad = getRadForEdges(curvedEdges) # this will return the radius for each of the 2 curved edges

    if abs(edgeRad[1] - edgeRad[0]) > .001:
        print("Different radii for the 2 curved edges!")
        print("Radius 1: "+str(edgeRad[0]))
        print("Radius 2: "+str(edgeRad[1]))
        return

    arcLength = curvedEdges[0].getLength()
    arcTheta = arcLength/edgeRad[0]
    if arcTheta-.1 > math.pi:
        print("Application area should span no more than 180 degrees!")
        print("Application span is: "+str(arcTheta*180./math.pi))
        return
    if distrType == "Gencoz":
        if abs(arcTheta) - math.pi > .1:
            print("Gencoz distribution is only supported for 180 degree span!")
            return
    theta2 =  arcTheta/2.0
    theta1 = -theta2

    #tmpPart = myModel.createPart("tmpPart")
    _target = apex.entityCollection()
    _target.append(curvedEdges[0])
    nodes = apex.mesh.createNodeByCurveArcCenter(_target)
    p1 = apex.Coordinate(nodes[0].getX(), nodes[0].getY(), nodes[0].getZ())
    #nn = []
    #for node in nodes:
    #	nn.append(node.getId())

    _target = apex.entityCollection()
    _target.append(curvedEdges[1])
    nodes = apex.mesh.createNodeByCurveArcCenter(_target)
    p2 = apex.Coordinate(nodes[0].getX(), nodes[0].getY(), nodes[0].getZ())
    #for node in nodes:
    #	nn.append(node.getId())

    #solid = selFace.getBody()
    #part  = solid.getParent()
    #meshes = part.getMeshes()
    #meshesToDelete = apex.entityCollection()
    #for mesh in meshes:
    #	elms = mesh.getElements()
    #	if elms.len() > 0: continue
    #	nodes = mesh.getNodes()
    #	for node in nodes:
    #		if node.getId() in nn:
    #			meshesToDelete.append(mesh)
    #			break
    #if meshesToDelete.len() > 0:
    #	apex.deleteEntities(meshesToDelete)

    # set up a coordinate system using the first edge
    verts = curvedEdges[0].getVertices()
    Vs = createVector3D(p1, verts[0])  # vector is from p1 to verts[0]
    sv = makeUnitVector(Vs)

    Vn = createVector3D(p2, p1) # vector is from p2 to p1
    nv = makeUnitVector(Vn)

    tv = crossProduct(nv, sv)

    elmFaces = getElmFacesForFace(selFace)  # this is a work-around, no method to get 3d elm faces from solid face
    if elmFaces is None: return

    geomFaceArea = selFace.getArea()
    elmFaceArea  = 0.0
    elmFaceData  = []
    for face in elmFaces:
        data = getFaceData(face)
        elmFaceArea += data[0]
        elmFaceData.append([face, data[0], data[1]])

    #pBar = 2*load/math.pi/thick/edgeRad[0]
    lugDescription = ""
    lugDescription = "Lug distribution by: "+distrType+" - Lug load: "+str(load)
    if distrType == "Cosine^2":
        pBar = load/thick/edgeRad[0]
        d    = math.sin(theta2)*(math.cos(theta2)**2 + 2.0)/3.0
        d   -= math.sin(theta1)*(math.cos(theta1)**2 + 2.0)/3.0
        pBar = pBar/d
    elif distrType == "Gencoz":
        pBar = 2*load/thick/edgeRad[0]/math.pi
    else:
        pBar = 2*load/thick/edgeRad[0]
        pBar = pBar/((.5*math.sin(2*theta2) + theta2) - (0.5*math.sin(2*theta1) + theta1))

    if areaAdjustment:
        pBar = pBar*geomFaceArea/elmFaceArea

    print("Lug thickness: "+str(thick))
    print("Lug radius: "+str(edgeRad[0]))
    print("Application region: "+str(arcTheta*180./math.pi))
    print("Geometry face area: "+str(geomFaceArea))
    print("FE face area: "+str(elmFaceArea))
    print("Elm faces: "+str(len(elmFaces)))
    if areaAdjustment:
        print("Area adjustment factor: "+str(geomFaceArea/elmFaceArea))

    for data in elmFaceData:
        # data = [face, area, centerPoint]

        Vc = createVector3D(p1, data[2])

        ds = dotProduct(sv, Vc) # dd is the component of Vc on the sv axis and can be used to calculate theta
        ds = ds*Vc.getLength()

        dt = dotProduct(tv, Vc)
        dt = dt*Vc.getLength()

        # now do dot between sv and ds?
        length = ds*ds + dt*dt
        length = math.sqrt(length)

        ds = ds/length
        dt = dt/length

        cosTheta = ds
        theta = math.acos(cosTheta)
        #theta -= math.pi/2.0
        theta -= theta2

        if distrType == "Cosine^2":
            cosTheta = math.cos(theta)
            p = pBar*cosTheta**2
        elif distrType == "Gencoz":
            d  = math.cos(theta)
            nn = [5, 9, 13, 17, 21]
            for n in nn:
                d -= 5.*math.cos(n*theta)/14./(n-1)/(n-8)
            nn = [3, 7, 11, 15, 19]
            for n in nn:
                d -= 2.*math.cos(n*theta)/5./(4-n)**2
            p = pBar*d
        else:
            cosTheta = math.cos(theta)
            p = pBar*cosTheta

        _target = apex.entityCollection()
        _target.append(data[0])
        pressure = apex.environment.createPressure(
          name          = lbcName,
          target        = _target, 
          pressureValue = p,
          description   = lugDescription)




