# -*- 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.
#
"""
Created on Thu Apr 21 12:55:09 2016
Failure processing routines for composite design explorer
@author: e1851
"""
try:
    import numpy
except:
    print("numpy not installed")
import apex
import math
try:
    import matplotlib.pyplot as plt
except:
    print("matplotlib not installed")
import classicalLaminateTheory as CLT
import failureTheoriesCDE
# reload modules during development
from importlib import reload
reload(failureTheoriesCDE)


def failureEnvelope(plies,
                    abdInv,
                    mat1_dict,
                    mat8_dict,
                    failureTheory,
                    loading):
    ''' create failure results as 2d numpy array.
    theta,ply dimensions'''
    num_plies = len(plies)
    num_theta = 361
    theta = 0.0
    fxy = 0.0
    SR_table = numpy.full([num_theta, num_plies], 0.0)
    theta_count = 0
    # get total thickness
    total_thickness = 0.0
    for ply in plies:
        total_thickness += ply.getThickness()
    if loading == "membrane":
        print(f"Membrane failure using {failureTheory}")
    else:
        print(f"Bending failure using {failureTheory}")
    # initialise lists
    # loop around 360 degrees
    while theta <= 360.0:
        theta_rad = math.radians(theta)
        # apply unit force at angle theta
        fx = math.cos(theta_rad)
        fy = math.sin(theta_rad)
        if loading == "membrane":
            f_tensor = numpy.array([fx, fy, fxy, 0.0, 0.0, 0.0]).transpose()
        else:
            f_tensor = numpy.array([0.0, 0.0, 0.0, fx, fy, fxy]).transpose()
        # multiply force by abdInv
        strainTensor = numpy.dot(abdInv, f_tensor)
        # strainTensor2 = strainTensor[0]
        # print(f"element strain: {strainTensor}")
        # print(f"element strain2: {strainTensor2}")
        # calculate total thickness
        ply_count = 0
        # get mid ply strains for each ply starting z at -t/2
        z_bot = -total_thickness/2.0
        for ply in plies:
            z_top = z_bot + ply.getThickness()
            z_mid = (z_bot + z_top)/2.0
            ply_strain_x = float(strainTensor[0]) - z_mid * float(strainTensor[3])
            ply_strain_y = float(strainTensor[1]) - z_mid * float(strainTensor[4])
            ply_strain_xy = float(strainTensor[2]) - z_mid * float(strainTensor[5])
            # Rotate ply strains from xy to 12 system
            # Converting shear strain from engineering to tensor by
            # multiplying by 0.5
            angle = ply.getAngle()
            angle_rad = math.radians(angle)
            mohrs_return = CLT.mohrs_rotate(ply_strain_x,
                                            ply_strain_y,
                                            0.5 * ply_strain_xy,
                                            angle_rad)
            ply_strain_12 = numpy.array([mohrs_return[0],
                                         mohrs_return[1],
                                         mohrs_return[2]*-2.0])
            # print(f"Ply strain: {ply_strain_12}")
            # get material Q
            mat_id = ply.getMaterial().getName()
            material = mat1_dict.get(mat_id)
            if material is None:
                material = mat8_dict.get(mat_id)
                q = material.q
            # print(material)
            # stress in ply axes
            ply_stress_12 = numpy.dot(ply_strain_12, q).transpose()
            # print(f"Ply stress: {ply_stress_12}")
            z_bot = z_top
            # stresses and strains calculated, calculate failure
            # failure theory called with relevent stress/strain values
            # returning strength ratio
            if failureTheory == 'Max Strain':
                strength_ratio = failureTheoriesCDE.max_strain(ply_strain_12,
                                                               material)
            elif failureTheory == "Hill":
                strength_ratio = failureTheoriesCDE.hill(ply_stress_12,
                                                         material)
            elif failureTheory == "Tsai-Wu":
                strength_ratio = failureTheoriesCDE.tsai_wu(ply_stress_12,
                                                            material)

            SR_table[theta_count, ply_count] = strength_ratio
            ply_count += 1
        theta += 1
        theta_count += 1
    # print(f"SR_table {SR_table}")
    plt_data = processEnvelope(SR_table, failureTheory, loading)
    return plt_data


def processEnvelope(SR_table, failureTheory, loading):
    table_dims = (SR_table.shape)
    rows = table_dims[0]
    columns = table_dims[1]

    # Initialise settings
    angle_increment = 360 / (rows-1)
    angle = 0.0
    angle_list = []
    SR_list = []
    crit_ply_list = []
    row = 0
    while angle < 360.0:
        # for each angle scan for minimum SR and add to list
        column = 0
        angle_list.append(math.radians(angle))
        min_SR = 1.0e10
        while column < columns:
            SR = SR_table[row, column]
            if SR < min_SR:
                min_SR = SR
                crit_ply = column
            column += 1
        SR_list.append(min_SR)
        crit_ply_list.append(crit_ply)
        row += 1
        angle += angle_increment
    plt_data = (angle_list, SR_list)
    return plt_data
    # plt.show()


def list_process_min_rf(input_list):
    """ scans list of results to find minimum reserve factor and its associated
    element and ply number.\n
    Input list is made of tuple of EID,PLYID,RF.\n
    Output is tuple of EID,PLYID,RF for the minimum result."""
    rf_min = 1e33
    eid_min = 0
    plyid_min = 0
    for i in range(0, len(input_list)):
        eid = input_list[i][0]
        plyid = input_list[i][1]
        rf = input_list[i][2]
        if rf < rf_min:
            rf_min = rf
            eid_min = eid
            plyid_min = plyid
#    print(f"Element: {eid_min} Ply: {plyid_min} RF: {rf_min}")
    return[eid_min, plyid_min, rf_min]


def list_process_worst_ply(input_list):
    """ scans list of results to create list of worst plyid for each element.\n
    Input list is made of tuple of EID,PLYID,RF.\n
    Output is list of EID,WORST_PLYID for the minimum RF result."""
    output_list = []
    eid = input_list[0][0]  # initialise with first element id in list
    worst_ply = 0
    rf_min = 1e33
    for i in range(0, len(input_list)):
        if input_list[i][0] != eid:  # new element in list
            output_list.append([eid, worst_ply, rf_min])
            eid = input_list[i][0]
            worst_ply = 0  # initialise to 0 for each element
            rf_min = 1e33
        if input_list[i][2] < rf_min:
            rf_min = input_list[i][2]
            worst_ply = input_list[i][1]
    output_list.append([input_list[i][0], worst_ply, rf_min])
#    print(output_list)
    return output_list
