Forked from
Klaus Böhnlein / dune-microstructure
594 commits behind, 124 commits ahead of the upstream repository.
-
Klaus Böhnlein authoredKlaus Böhnlein authored
ClassifyMin.py 13.92 KiB
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 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...")
CaseCount = 0
epsilon = sys.float_info.epsilon #Machine epsilon
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]])
# print('Matrix A:', A)
# print('Matrix B:', B)
# print('Matrix rank of A:', np.linalg.matrix_rank(A))
# print('shape of A:', A.shape)
# print('shape of B:', B.shape)
# print('Matrix [A,B]:', np.c_[A, B])
# print('Matrix rank of [A,B]:', np.linalg.matrix_rank(np.c_[A, B]))
# print('shape of [A,B]:', C.shape)
#
# x = np.linalg.solve(A, B) # works only if A is not singular!!
# print("Solution LGS:", x)
# # print("sym Matrix", sym.Matrix(([A],[B])))
#
# # Test
# C = np.array([[1, 0], [0, 0]])
# d = np.array([5, 0])
# y = np.linalg.lstsq(C, d)[0]
# print("Solution LGS:", y)
# T = np.c_[C, d]
# print('T:', T)
# Trref = sym.Matrix(T).rref()[0]
# Out = np.array(Trref, dtype=float)
# print('rref:', Out)
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 print_Cases: print('P : parabolic case (determinant equal zero)')
print('P : parabolic case (determinant equal zero)')
# if print_Cases: print('P : parabolic case (determinant equal zero)')
# check if B is in range of A
# check if rank(A) == rank([A,B])
# OK this way? (or use Sympy?)
if np.linalg.matrix_rank(A) == np.linalg.matrix_rank(np.c_[A, B]):
if print_Cases: print('P1 (B is in the range of A)')
if (q12+q3)/2.0 <= -1.0*epsilon:
print('should not happen(not compatible with det = 0)')
if (abs(B[0]) < epsilon and abs(B[1]) < epsilon) and (q12+q3)/2.0 >= epsilon:
if print_Cases: print('P1.1')
a1 = 0.0
a2 = 0.0
type = 3
CaseCount += 1
if (abs(B[0]) >= epsilon or abs(B[1]) >= epsilon) and (q12+q3)/2.0 >= epsilon:
# Continuum of minimizers located along a line of negative slope in Lambda
if print_Cases: print('P1.2 (Continuum of minimizers located along a line of negative slope in Lambda) ')
# Just solve Aa* = b (alternatively using SymPy ?)
# we know that A is singular (det A = 0 ) here..
# & we know that there are either infinitely many solutions or a unique solution ...
# ---- determine one via Least Squares
# "If A is square and of full rank, then x (but for round-off error)
# is the “exact” solution of the equation. Else, x minimizes the
# Euclidean 2-norm || b-Ax ||. If there are multiple minimizing solutions,
# the one with the smallest 2-norm is returned. ""
# a = np.linalg.lstsq(A, B)[0] # TODO check is this Ok ?
a = np.linalg.lstsq(A, -B/2)[0] # TODO check is this Ok ? (UPDATE 13-10-21)
print("Solution LGS: a =", a)
a1 = a[0]
a2 = a[1]
type = 2
CaseCount += 1
else:
if print_Cases: print('P2 (B is not in the range of A)')
# local Minimizers occur on the boundary of Lambda...
# There are at most two local minima and they are either
# (b1_star, 0) or (0, b2_star)
# could also outsource this to another function..
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
# ------------------------------------ Elliptic Case -----------------------------------
if determinant >= epsilon:
if print_Cases: print('E : elliptic case (determinant greater zero)')
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
if prod >= epsilon:
if print_Cases: print('(E1) - inside Lambda ')
a1 = a1_star
a2 = a2_star
type = 1 # non-axial Minimizer
CaseCount += 1
if abs(prod) < epsilon: # same as prod = 0 ? or better use <=epsilon ?
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)))]
angle = math.atan2(e[1], e[0])
# 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, 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)
# -----------------------------------------------------------------------------------------------------------------------------------------------------------------