import numpy as np
import matplotlib.pyplot as plt
import sympy as sym
import math
import os
import subprocess
import fileinput
import re
import matlab.engine
import sys
# print(sys.executable)
# from subprocess import Popen, PIPE
# --------------------------------------------------
# 'classifyMin' classifies Minimizers by utilizing the result of
# Lemma1.6
#
#
#
#
# 'classifyMin_ana': (Special-Case : Lemma1.4)
# ..additionally assumes Poisson-ratio=0 => q12==0
#
#
#
# Output : MinimizingMatrix, Angle, Type, Curvature
def get_gstar(q1,q2,q12,q3,b1,b2):
# a = np.array([a1,a2])
b = np.array([b1,b2])
H = np.array([[2*q1, q12+2*q3], [q12+2*q3,2*q2] ])
A = np.array([[q1,q12/2], [q12/2,q2] ])
# print('det(H)=', np.linalg.det(H))
# check if g* is in G^*_R^2
tmp = A.dot(b)
## compute inverse of H :
inv_H = np.linalg.inv(H)
g_star = 2*inv_H.dot(tmp)
# print('g_star=', g_star)
return g_star
def determine_b(q1,q2,q12,q3,g_star):
##
# Input: g_star
# Output : b such that g_star is minimizer, i.e A*b = (1/2)*H*g_star
# q1=1;
# q2=2;
# q12=1/2;
# q3=((4*q1*q2)**0.5-q12)/2;
# H=[2*q1,q12+2*q3;q12+2*q3,2*q2];
H = np.array([[2*q1, q12+2*q3], [q12+2*q3,2*q2] ])
A = np.array([[q1,1/2*q12], [1/2*q12,q2] ])
abar = np.array([q12+2*q3, 2*q2])
rhs = (1/2)*H.dot(g_star)
print('rhs:', rhs)
b = np.linalg.lstsq(A, rhs)[0]
print('b',b)
b1=b[0]
b2=b[1]
return b
## ---------------
def get_minimizer(q1,q2,q3,b1,b2):
# In the case if q12 == 0:
quotient = (q1*q2-q3**2)
g_star = np.array([(q1*q2*b1-q3*q2*b2)/quotient, (q1*q2*b2-q3*q1*b1)/quotient])
print('g_star=', g_star)
return g_star
def energy(a1,a2,q1,q2,q12,q3,b1,b2):
a = np.array([a1,a2])
b = np.array([b1,b2])
H = np.array([[2*q1, q12+2*q3], [q12+2*q3,2*q2] ])
A = np.array([[q1,1/2*q12], [1/2*q12,q2] ])
tmp = H.dot(a)
# print('H',H)
# print('A',A)
# print('b',b)
# print('a',a)
# print('tmp',tmp)
tmp = (1/2)*a.dot(tmp)
# print('tmp',tmp)
tmp2 = A.dot(b)
# print('tmp2',tmp2)
tmp2 = 2*a.dot(tmp2)
# print('tmp2',tmp2)
energy = tmp - tmp2
# print('energy',energy)
# energy_axial1.append(energy_1)
return energy
def determinant(q1,q2,q3,q12): # TODO General:Matrix
return q1*q2 - (q3**2 + 2*q3*q12 + q12**2)
def harmonicMean(mu_1, beta, theta):
return mu_1*(beta/(theta+((1-theta)*beta)))
def arithmeticMean(mu_1, beta, theta):
return mu_1*((1-theta)+theta*beta)
def prestrain_b1(rho_1, beta, alpha, theta):
return (3.0*rho_1/2.0)*(1-(theta*(1+alpha)))
# return (3.0*rho_1/2.0)*beta*(1-(theta*(1+alpha)))
def prestrain_b2(rho_1, beta, alpha, theta):
return (3.0*rho_1/(2.0*((1.0-theta) + theta*beta)))*(1-theta*(1+beta*alpha))
# return (3.0*rho_1/(4.0*((1.0-theta) + theta*beta)))*(1-theta*(1+beta*alpha))
# Define function to be minimized
def f(a1, a2, q1, q2, q3, q12, b1, b2):
A = np.array([[q1, q3 + q12/2.0], [q3 + q12/2.0, q2]])
B = np.array([-2.0*q1*b1-q12*b2, -2.0*q2*b2-q12*b1])
a = np.array([a1, a2])
tmp = np.dot(A, a)
tmp2 = np.dot(a, tmp)
tmpB = np.dot(B, a)
return tmp2 + tmpB + q1*(b1**2) + q2*(b2**2) + q12*b1*b2
# ---- Alternative Version using alpha,beta,theta ,mu_1,rho_1
def classifyMin_ana(alpha,beta,theta,q3,mu_1,rho_1,print_Cases=False, print_Output=False):
# Assumption of Classification-Lemma1.6:
# 1. [b3 == 0]
# 2. Q is orthotropic i.e. q13 = q31 = q23 = q32 == 0
# 3. This additionally assumes that Poisson-Ratio = 0 => q12 == 0
q12 = 0.0
q1 = (1.0/6.0)*harmonicMean(mu_1, beta, theta)
q2 = (1.0/6.0)*arithmeticMean(mu_1, beta, theta)
# print('q1: ', q1)
# print('q2: ', q2)
b1 = prestrain_b1(rho_1, beta, alpha,theta)
b2 = prestrain_b2(rho_1, beta, alpha,theta)
# print('alpha:',alpha)
# print('beta:',beta)
# print('theta:',theta)
return classifyMin(q1, q2, q3, q12, b1, b2, print_Cases, print_Output)
# Matrix Version that just gets matrices Q & B
def classifyMin_mat(Q,B,print_Cases=False, print_Output=False):
q1 = Q[0][0]
q2 = Q[1][1]
q3 = Q[2][2]
q12 = Q[0][1]
b1 = B[0]
b2 = B[1]
b3 = B[2]
return classifyMin(q1, q2, q3, q12, b1, b2, print_Cases, print_Output)
# --------------------------------------------------------------------
# Classify Type of minimizer 1 = R1 , 2 = R2 , 3 = R3 # before : destinction between which axis.. (4Types )
# where
# R1 : unique local (global) minimizer which is not axial
# R2 : continuum of local (global) minimizers which are not axial
# R3 : one or two local (global) minimizers which are axial
# Partition given by
# R1 = E1
# R2 = P1.2
# R3 = E2 U E3 U P1.1 U P2 U H
# -------------------------------------------------------------------
def classifyMin(q1, q2, q3, q12, b1, b2, print_Cases=False, print_Output=False): #ClassifyMin_hom?
# Assumption of Classification-Lemma1.6:
# 1. [b3 == 0]
# 2. Q is orthotropic i.e. q13 = q31 = q23 = q32 == 0
# TODO: check if Q is orthotropic here - assert()
if print_Output: print("Run ClassifyMin_NEW...")
CaseCount = 0
epsilon = sys.float_info.epsilon #Machine epsilon
# print('epsilon:',epsilon)
b = np.array([b1,b2])
H = np.array([[2*q1, q12+2*q3], [q12+2*q3,2*q2] ])
A = np.array([[q1,q12/2], [q12/2,q2] ])
# B = np.array([-2.0*q1*b1-q12*b2, -2.0*q2*b2-q12*b1])
# A = np.array([[q1, q3 + q12/2.0], [q3 + q12/2.0, q2]])
determinant = np.linalg.det(H)
# print('determinant:',determinant)
# determinant = q1*q2 - (q3**2 + 2*q3*q12 + q12**2)
if print_Cases: print("determinant:", determinant)
# Define values for axial-Solutions (b1*,0) & (0,b2*)
# b1_star = (2.0*q1*b1 + b2*q12)/(2*q1)
# b2_star = (2.0*q2*b2 + b1*q12)/(2*q2)
# ------------------------------------ Parabolic Case -----------------------------------
# if abs(determinant) < epsilon:
if abs(determinant) < epsilon:
if print_Cases: print('P : parabolic case (determinant equal zero)')
print('P : PARABOLIC CASE (determinant equal zero)')
# ------------------------------------ Elliptic Case -----------------------------------
if determinant >= epsilon:
if print_Cases: print('E : elliptic case (determinant greater zero)')
g_star = get_gstar(q1,q2,q12,q3,b1,b2)
# a1_star = (b1*(q12**2) + 2*b1*q3*q12 - 4*b1*q1*q2 + 4*b2*q2*q3) / \
# (4*(q3**2) + 4*q3*q12 + (q12**2) - 4*q1*q2)
# a2_star = (b2*(q12**2) + 2*b2*q3*q12 + 4*b1*q1*q3 - 4*b2*q1*q2) / \
# (4*(q3**2) + 4*q3*q12 + (q12**2) - 4*q1*q2)
# prod = a1_star*a2_star
prod = g_star[0]*g_star[1]
if prod >= epsilon:
if print_Cases: print('(E1) - inside Lambda ')
# a1 = a1_star
# a2 = a2_star
# type = 1 # non-axial Minimizer
# CaseCount += 1
a1 = g_star[0]
a2 = g_star[1]
CaseCount += 1
if prod < epsilon: # same as prod = 0 ? or better use <=epsilon ?
if abs(((b2**2)/q1 - (b1**2)/q2)) < epsilon :
print('((b2**2)/q1 - (b1**2)/q2)', ((b2**2)/q1 - (b1**2)/q2))
print('two minimizers')
a1 = b1
a2 = 0.0
# type = 3 # 1
CaseCount += 1
else:
#compare energy values
if energy(b1, 0, q1, q2, q3, q12, b1, b2) > energy(0, b2, q1, q2, q3, q12, b1, b2):
a1 = 0
a2 = b2
CaseCount += 1
else :
a1 = b1
a2 = 0
CaseCount += 1
else :
# g_star = get_gstar(q1,q2,q12,q3,b1,b2)
# prod = g_star[0]*g_star[1]
if abs((b2**2)/q1 - (b1**2)/q2) < epsilon :
print('two minimizers')
a1 = b1
a2 = 0.0
# type = 3 # 1
CaseCount += 1
else:
#compare energy values
if energy(b1, 0, q1, q2, q3, q12, b1, b2) > energy(0, b2, q1, q2, q3, q12, b1, b2):
a1 = 0
a2 = b2
CaseCount += 1
else :
a1 = b1
a2 = 0
CaseCount += 1
#
# if print_Cases: print('(E2) - on the boundary of Lambda ')
# a1 = a1_star
# a2 = a2_star
# type = 3 # could check which axis: if a1_star or a2_star close to zero.. ?
# CaseCount += 1
#
# if prod <= -1.0*epsilon:
# if print_Cases: print('(E3) - Outside Lambda ')
# if f(b1_star, 0, q1, q2, q3, q12, b1, b2) < f(0, b2_star, q1, q2, q3, q12, b1, b2):
# a1 = b1_star
# a2 = 0.0
# type = 3 # 1
# CaseCount += 1
# if f(b1_star, 0, q1, q2, q3, q12, b1, b2) > f(0, b2_star, q1, q2, q3, q12, b1, b2):
# a1 = 0
# a2 = b2_star
# type = 3 # 2
# CaseCount += 1
#
# # TODO Problem: angle depends on how you choose... THE angle is not defined for this case
# if f(b1_star, 0, q1, q2, q3, q12, b1, b2) == f(0, b2_star, q1, q2, q3, q12, b1, b2):
# # Two Minimizers pick one
# a1 = b1_star
# a2 = 0.0
# type = 3 # 4
# CaseCount += 1
# ------------------------------------ Hyperbolic Case -----------------------------------
# if determinant <= -1.0*epsilon:
# if print_Cases: print('H : hyperbolic case (determinant smaller zero)')
# # One or two minimizers wich are axial
# type = 3 # (always type 3)
# if f(b1_star, 0, q1, q2, q3, q12, b1, b2) < f(0, b2_star, q1, q2, q3, q12, b1, b2):
# a1 = b1_star
# a2 = 0.0
# # type = 3 # 1
# CaseCount += 1
# if f(b1_star, 0, q1, q2, q3, q12, b1, b2) > f(0, b2_star, q1, q2, q3, q12, b1, b2):
# a1 = 0
# a2 = b2_star
# # type = 3 # 2
# CaseCount += 1
# # TODO can add this case to first or second ..
# if f(b1_star, 0, q1, q2, q3, q12, b1, b2) == f(0, b2_star, q1, q2, q3, q12, b1, b2):
# # Two Minimizers pick one
# a1 = b1_star
# a2 = 0.0
# # type = 3 # 4
# CaseCount += 1
# ---------------------------------------------------------------------------------------
if (CaseCount > 1):
print('Error: More than one Case happened!')
# compute a3
# a3 = math.sqrt(2.0*a1*a2) # never needed?
# print('a1:', a1)
# print('a2:', a2)
# compute the angle <(e,e_1) where Minimizer = kappa* (e (x) e)
# e = [math.sqrt((a1/(a1+a2))), math.sqrt((a2/(a1+a2)))]
e = [((a1/(a1+a2)))**0.5, ((a2/(a1+a2)))**0.5]
angle = math.atan2(e[1], e[0])
type = 1 # ToDO..
# compute kappa
kappa = (a1 + a2)
# Minimizer G
# Minimizer = np.array([[a1, math.sqrt(a1*a2)], [math.sqrt(a1*a2), a2]],dtype=object)
Minimizer = np.array([[a1, (a1*a2)**0.5], [(a1*a2)**0.5, a2]],dtype=object)
# Minimizer = np.array([[a1, math.sqrt(a1*a2)], [math.sqrt(a1*a2), a2]])
# MinimizerVec = np.array([a1, a2],dtype=object)
MinimizerVec = np.array([a1, a2])
if print_Output:
print('--- Output ClassifyMin ---')
print("Minimizing Matrix G:")
print(Minimizer)
print("angle = ", angle)
print("type: ", type)
print("kappa = ", kappa)
return Minimizer, angle, type, kappa
# return MinimizerVec, angle, type, kappa #return Minimizer Vector instead
# ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# ---------------------------------------------- Main ---------------------
# --- Input Parameters ----
# mu_1 = 1.0
# rho_1 = 1.0
# alpha = 9.0
# beta = 2.0
# theta = 1.0/8.0
# # define q1, q2 , mu_gamma, q12
# # 1. read from Cell-Output
# # 2. define values from analytic formulas (expect for mu_gamma)
# q1 = (1.0/6.0)*harmonicMean(mu_1, beta, theta)
# q2 = (1.0/6.0)*arithmeticMean(mu_1, beta, theta)
# # TEST
# q12 = 0.0 # (analytical example)
# # q12 = 12.0 # (analytical example)
# set mu_gamma to value or read from Cell-Output
# mu_gamma = q1 # TODO read from Cell-Output
# b1 = prestrain_b1(rho_1, beta, alpha, theta)
# b2 = prestrain_b2(rho_1, beta, alpha, theta)
# print('---- Input parameters: -----')
# print('mu_1: ', mu_1)
# print('rho_1: ', rho_1)
# print('alpha: ', alpha)
# print('beta: ', beta)
# print('theta: ', theta)
# print("q1: ", q1)
# print("q2: ", q2)
# print("mu_gamma: ", mu_gamma)
# print("q12: ", q12)
# print("b1: ", b1)
# print("b2: ", b2)
# print('----------------------------')
# # print("machine epsilon", sys.float_info.epsilon)
#
#
# # ------- Options --------
# print_Cases = True
# print_Output = True
# G, angle, type, kappa = classifyMin(q1, q2, mu_gamma, q12, b1, b2, print_Cases, print_Output)
#
# G, angle, type, kappa = classifyMin_ana(alpha, beta, theta, mu_gamma, q12, print_Cases, print_Output)
#
# Out = classifyMin_ana(alpha, beta, theta, mu_gamma, q12, print_Cases, print_Output)
# print('TEST:')
# Out = classifyMin_ana(alpha, beta, theta)
# print('Out[0]', Out[0])
# print('Out[1]', Out[1])
# print('Out[2]', Out[2])
# print('Out[3]', Out[3])
# #supress certain Outout..
# _,_,T,_ = classifyMin_ana(alpha, beta, theta, mu_gamma, q12, print_Cases, print_Output)
# print('Output only type..:', T)
# Test = f(1,2 ,q1,q2,mu_gamma,q12,b1,b2)
# print("Test", Test)
# -----------------------------------------------------------------------------------------------------------------------------------------------------------------