#
# Calibrate constituent properties of a UD composite using Mori-Tanaka
# 
import numpy as np
from scipy.optimize import minimize
from scipy import linalg
#
# Utility functions
# 
#Get stiffness matrix from engineering constants
def getC_engconsts(elasconsts):
    if elasconsts.size == 5:
        E1 = elasconsts[0]
        E2 = elasconsts[1]
        v12 = elasconsts[2]
        v23 = elasconsts[3]
        G12 = elasconsts[4]   
    else:
        E1 = elasconsts[0]
        E2 = E1
        v12 = elasconsts[1]
        v23 = v12
        G12 = E1*0.5/(1+v12)
    G23 = E2*0.5/(1+v23)    
    A = np.matrix([[1/E1,-v12/E1,-v12/E1],[-v12/E1,1/E2,-v23/E2],[-v12/E1,-v23/E2,1/E2]])
    Ainv = np.linalg.inv(A)
    u = np.concatenate((Ainv, np.zeros((3,3))), axis=1)
    shear = np.array([[G12,0,0],[0,G12,0],[0,0,G23]])
    b = np.concatenate((np.zeros((3,3)),np.array([[G12,0,0],[0,G12,0],[0,0,G23]])), axis=1)
    C = np.concatenate((u,b),axis=0)
    return C
#
#Get engineering constants from stiffness matrix
#
def getengconsts(C):
    Cinv = np.linalg.inv(C)
    E1 = 1/Cinv[0,0]
    E2 = 1/Cinv[1,1]   
    G12 = 1/Cinv[3,3]   
    G23 = 1/Cinv[5,5]
    gamma12 = -Cinv[1,0]*E1   
    gamma23 = -Cinv[2,1]*E2
    y = np.array([E1,E2,gamma12,gamma23,G12])
    return y
#
#Compute effective stiffness matrix of UD composite using Mori-Tanaka homogenization
#
def MT_cylinder(elasconsts_matrix,elasconsts_fiber,vff):    
    #volume fraction vff, vfm
    vfm = 1-vff
    v0 = elasconsts_matrix[1]
    Cm = getC_engconsts(elasconsts_matrix)    
    Cf = getC_engconsts(elasconsts_fiber)   
    Cinv = np.linalg.inv(Cm)
    #Eshelby tensor for cylindrical inclusion
    s1 = (5-4*v0)/(8*(1-v0))
    s2 = (4*v0-1)/(8*(1-v0))
    s3 = v0/(2*(1-v0))
    s4 = (3-4*v0)/(8*(1-v0))
    s5 = 0.25
    s = np.matrix([[0,0,0,0,0,0],
          [s3,s1,s2,0,0,0],
          [s3,s2,s1,0,0,0],
          [0,0,0,2*s5,0,0],
          [0,0,0,0,2*s5,0],
          [0,0,0,0,0,2*s4]])
    I = np.matrix([[1,0,0,0,0,0],
     [0,1,0,0,0,0],
     [0,0,1,0,0,0],
     [0,0,0,1,0,0],
     [0,0,0,0,1,0],
     [0,0,0,0,0,1]])
    AI = np.linalg.inv(s.dot(Cinv.dot(Cf)-I)+I)
    t1 = vfm*Cm+vff*Cf.dot(AI)
    t2 = np.linalg.inv(vfm*I+vff*AI)
    Ceff = t1.dot(t2)
    return getengconsts(Ceff)
#
# Function
#
def f_with_args(design_vars,elasconsts_matrix,elasconsts_fiber,vff,elasconsts_lamina,isDesignVar,weights):
    arr = np.concatenate([np.concatenate([elasconsts_matrix,elasconsts_fiber]),np.array([vff])])   
    arr2 = arr
    index = 0
    for i in range(arr.size):
        if isDesignVar[i] == 1:
            arr2[i] = design_vars[index]
            index = index + 1 
    y = MT_cylinder(arr2[0:2],arr2[2:7],arr2[7])
    err = np.linalg.norm(np.divide(np.multiply(weights,y-elasconsts_lamina),elasconsts_lamina))
    return err
#
# User Input
#
#Lamina constants (target values)
elasconsts_lamina = np.array([55000,9500,0.33,0.45,5500])
#weights of targets
weights = np.array([1,1,1,1,1])
#Matrix elastic constants
#Em, nu_m
elasconsts_matrix = np.array([3500,0.35])
#Fiber elastic constants (transversely isotropic)
#E1 E2 v12 v23 G12
elasconsts_fiber = np.array([74000,74000,0.2,0.2,30800])
#Fiber volume fraction
vff = 0.68
#
#flags to mark design variables
#Em, nu_m, E1, E2, v12, v23, G12, vff -- set to zero to keep fixed
isDesignVar = np.array([1,1,1,1,1,1,1,0])
#
# Initialize and optimize
#
#Starting value
arr = np.concatenate([np.concatenate([elasconsts_matrix,elasconsts_fiber]),np.array([vff])])
var0 = np.array([])
for i in range(arr.size):
        if isDesignVar[i] == 1:
            var0 = np.concatenate([var0,np.array([arr[i]])])

#Optimize
res = minimize(f_with_args, var0, method='nelder-mead',args=(elasconsts_matrix,elasconsts_fiber,vff,elasconsts_lamina,isDesignVar,weights), options={'xatol': 1e-8, 'disp': True})
var = res.x

#Output
arr = np.concatenate([np.concatenate([elasconsts_matrix,elasconsts_fiber]),np.array([vff])]) 
arr2 = arr
index = 0
for i in range(arr.size):
        if isDesignVar[i] == 1:
            arr2[i] = var[index]
            index = index + 1 
y = MT_cylinder(arr2[0:2],arr2[2:7],arr2[7])
print ("Lamina Elastic Constants (calibrated)")
print ("E1=",y[0])
print ("E2=",y[1])
print ("nu12=",y[2])
print ("nu23=",y[3])
print ("G12=",y[4])
##Results
print ("Volume fraction")
print ("vf=",arr2[7])
print ("Matrix Elastic Constants")
print ("E=",arr2[0])
print ("nu=",arr2[1])
print ("Fiber Elastic Constants")
print ("E1=",arr2[2])
print ("E2=",arr2[3])
print ("nu12=",arr2[4])
print ("nu23=",arr2[5])
print ("G12=",arr2[6])