added plot.py to create angle-curvature plot from Q and B

microstructure_testsuite/plot.py

+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Wed Jul  6 13:17:28 2022
+
+@author: stefan
+"""
+import numpy as np
+import matplotlib.pyplot as plt
+from helper_functions import *
+
+
+def energy(kappa,alpha,Q,B)  :
+    G=kappa*np.array([[np.cos(alpha)**2],[np.sin(alpha)**2],[np.sqrt(2)*np.cos(alpha)*np.sin(alpha)]])-B
+    return np.matmul(np.transpose(G),np.matmul(Q,G))[0,0]
+
+def xytokappaalpha(x,y):
+   
+    if y>0:
+        return [np.sqrt(x**2+y**2), np.abs(np.arctan2(y,x))]
+    else:
+        return [-np.sqrt(x**2+y**2), np.abs(np.arctan2(y,x))]
+
+
+#
+case = -1
+if case == 0:        # Beispiel Figure 10 (c)    
+    q1=1; q2=2; q12=1/2; q3=0.5*(np.sqrt(4*q1*q2)-q12);
+    Q=np.array([[q1,q12/2,0],[q12/2,q2,0],[0,0,q3]])
+    B=np.array([[0.491],[0.347],[0]])
+elif case == 1: # Beispiel Figure 10 (b)    
+    q1=1; q2=2; q12=0; q3=1;
+    Q=np.array([[q1,q12/2,0],[q12/2,q2,0],[0,0,q3]])
+    B=np.array([[2],[1.5],[0]])       
+elif case == 2:  # Beispiel Figure 10 (a)    
+    q1=1; q2=2; q12=0; q3=1;
+    Q=np.array([[q1,q12/2,0],[q12/2,q2,0],[0,0,q3]])
+    B=np.array([[0],[0.5],[0]])
+elif case == 3: # Homogener Fall
+    q1=1; q2=1; q12=0; q3=1;
+    Q=np.array([[q1,q12/2,0],[q12/2,q2,0],[0,0,q3]])
+    B=np.array([[.5],[.5],[0]])
+elif case == 4: # Section 6.1
+    theta=0.2
+    rho=1
+    q1=1/(3*(2-theta)); q2=(1+theta)/6; q12=0; q3=q1;
+    Q=np.array([[q1,q12/2,0],[q12/2,q2,0],[0,0,q3]])
+    B=np.array([[rho*3/2*(1-theta)],[rho*3/2*(1-theta)/(1+theta)],[0]])
+elif case==-1: # Read from outputs
+    QFilePath = os.path.dirname(os.getcwd()) + '/outputs/QMatrix.txt'
+    BFilePath = os.path.dirname(os.getcwd())+ '/outputs/BMatrix.txt'
+    Q, B = ReadEffectiveQuantities(QFilePath,BFilePath)
+    B=np.transpose([B])
+# 
+
+#    
+N=200
+h=0.0005
+X=np.zeros([N,N])
+for i in range(0,N): 
+    for j in range(0,N):     
+        x=(i-N/2)*h
+        y=(j-N/2)*h
+        K=xytokappaalpha(x,y)
+        X[N-1-j,i]=energy(K[0],K[1],Q,B) # N-1-j damit die x,y Axen stimmen
+
+fig = plt.figure(figsize=(6,6))
+plt.imshow(np.log(X-np.min(X)+0.0001)) # normalize to min = 0 and log scale to emphasize energy landscape
+# TODO: Beschriftung der Axen sollte von [-h*N/2, h*N/2] sein!
+# Kreis mit Radius 1 einzeichnen
+plt.show()
+