update

geometries/material_neukamm.py

@@ -5,25 +5,40 @@ import math
 # x[0] : y1-component
 # x[1] : y2-component
 # x[2] : x3-component
+#    --- replace with your definition of indicatorFunction:
+###############
+# Wood
+###############
+# def f(x):
+#     theta=0.25
+    # if ((abs(x[0]) < theta/2) and x[2]<0.25):
+    #     return 1    #latewood
+    # elif ((abs(x[0]) > theta/2) and x[2]<0.25):
+    #     return 2    #earlywood
+    # else :
+    #     return 0    #Phase3
+
+# def b1(x):
+#     return [[.5, 0, 0], [0,1,0], [0,0,0]]
+
+# def b2(x):
+#     return [[.4, 0, 0], [0,.4,0], [0,0,0]]
+
+# def b3(x):
+#     return [[0, 0, 0], [0,0,0], [0,0,0]]
+###############
+# Cross
+###############
 def f(x):
     theta=0.25
     factor=1
-    # --- replace with your definition of indicatorFunction:
-    #
-    # if ((abs(x[0]) < theta/2) and x[2]<0):
-    #     return 1    #Phase1
-    # else :
-    #     return 0    #Phase3
-    
-    # 
-    if ((abs(x[0]) < theta/2) and x[2]<-1/2+theta):
+    if (x[0] <-1/2+theta and x[2]<-1/2+theta):
         return 1    #Phase1
-    elif ((abs(x[1]) < factor*theta/2) and x[2]>1/2-factor*theta):
+    elif (x[1]< -1/2+theta and x[2]>1/2-theta):
         return 2    #Phase2
     else :
         return 0    #Phase3
 
-
 def b1(x):
     return [[1, 0, 0], [0,1,0], [0,0,1]]
 

inputs/cellsolver.parset

@@ -27,7 +27,7 @@ cellDomain=1
 ## {start,finish} computes on all grid from 2^(start) to 2^finish refinement
 #----------------------------------------------------
 
-numLevels=2 4
+numLevels=2 5
 #numLevels =  1 1   # computes all levels from first to second entry
 #numLevels =  2 2   # computes all levels from first to second entry
 #numLevels =  1 3   # computes all levels from first to second entry
@@ -64,13 +64,14 @@ lambda1=1.0
 
 
 # better: pass material parameters as a vector
-mu=80 80 60
-lambda=80 80 25
+# Poisson ratio nu = 1/4 => lambda = 0.4*mu
+mu=1.2 1.2 1
+lambda=0.48 0.48 0.4
 #rho=1.0 0
 
 
 # ---volume fraction  (default value = 1.0/4.0)
-theta=0.25
+theta=0.5
 #theta = 0.3
 #theta = 0.75
 #theta = 0.125
@@ -96,8 +97,8 @@ material_prestrain_imp= "material_neukamm"   #(Python-indicator-function with sa
 
 
 # --- (Optional output) write Material / prestrain / Corrector functions to .vtk-Files (default=false):
-#write_materialFunctions = true
-#write_prestrainFunctions = true  # VTK norm of B ,
+write_materialFunctions = true
+write_prestrainFunctions = true  # VTK norm of B ,
 #write_VTK = true
 
 

microstructure_testsuite/plot.py

@@ -48,15 +48,17 @@ elif case == 4: # Section 6.1
     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'
+    DataPath = os.path.dirname(os.getcwd())+'/outputs'
+    #DataPath='/home/stefan/DUNE/dune-microstructure/outputs/Results/4x4/1'
+    QFilePath = DataPath + '/QMatrix.txt'
+    BFilePath = DataPath + '/BMatrix.txt'
     Q, B = ReadEffectiveQuantities(QFilePath,BFilePath)
     Q=0.5*(np.transpose(Q)+Q) # symmetrize
     B=np.transpose([B])
 # 
 
 #    
-length=80
+length=0.4
 N=200
 h=length/N
 E=np.zeros([N,N])
@@ -75,7 +77,7 @@ ax.set_aspect(1)
 ax.set_xticks([-length/4,0,length/4])
 ax.set_yticks([])
 #pcm = plt.pcolor(X,Y,E, norm=colors.LogNorm(vmin=E.min(), vmax=E.max()), cmap='winter', shading='auto')
-pcm = plt.pcolor(X,Y,E, norm=colors.PowerNorm(gamma=0.25), cmap='brg')
+pcm = plt.pcolor(X,Y,E, norm=colors.PowerNorm(gamma=0.2), cmap='brg')
 plt.colorbar(pcm, extend='max')
 #plt.imshow(np.log(E-np.min(E)+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!

outputs/Results/2x2/1/BMatrix.txt

-1 1 -0.0144287202371896055
-1 2 0.014428690223795131
-1 3 -2.0320334857961725e-11

outputs/Results/2x2/1/QMatrix.txt

-1 1 12.4813265250416485
-1 2 2.0831795458888962
-1 3 3.83811740666402797e-10
-2 1 2.0831795513179836
-2 2 12.481326525238547
-2 3 1.78942192122814843e-10
-3 1 -6.40561063029264521e-10
-3 2 -1.00945173229121783e-09
-3 3 10.3922965257922275

outputs/Results/2x2/1/output.txt

-material_prestrain used: material_neukamm
------ Input Parameters -----: 
-alpha: 8
-gamma: 1
-theta: 0.25
-beta: 3
-material parameters: 
-mu1: 1
-mu2: 3
-lambda1: 1
-lambda2: 3
-----------------------------: 
-Number of Elements in each direction: [4,4,4]
-size of FiniteElementBasis: 240
-Solver-type used:  CG-Solver
----------- OUTPUT ----------
- --------------------
-Corrector-Matrix M_1: 
--0.00922212 -2.49142e-10 0
--2.49142e-10 -0.00494017 0
-0 0 0
-
- --------------------
-Corrector-Matrix M_2: 
-0.00494016 -2.03726e-10 0
--2.03726e-10 0.00922211 0
-0 0 0
-
- --------------------
-Corrector-Matrix M_3: 
-1.90697e-10 -4.74653e-11 0
--4.74653e-11 1.93256e-10 0
-0 0 0
-
- --------------------
-Effective Matrix Q: 
-12.4813 2.08318 3.83812e-10
-2.08318 12.4813 1.78942e-10
--6.40561e-10 -1.00945e-09 10.3923
-
---- Prestrain Output --- 
-B_hat: -0.150032 0.150032 -2.16498e-10
-B_eff: -0.0144287 0.0144287 -2.03203e-11 (Effective Prestrain)
------------------------- 
-q1=12.4813
-q2=12.4813
-q3=10.3923
-q12=2.08318
-b1=-0.0144287
-b2=0.0144287
-b3=-2.03203e-11
-b1_hat: -0.150032
-b2_hat: 0.150032
-b3_hat: -2.16498e-10
-mu_gamma=10.3923
-q_onetwo=2.083180
----------------------------------------------------------------------------------------------------------
-  Levels     |      q1      |      q2      |      q3      |      b1      |      b2      |      b3      | 
----------------------------------------------------------------------------------------------------------
-     2       & 1.24813e+01  & 1.24813e+01  & 1.03923e+01  & -1.44287e-02 & 1.44287e-02  & -2.03203e-11 & 
----------------------------------------------------------------------------------------------------------

outputs/Results/2x2/1/parameter

-mu=80 60
-lambda=80 25
-rho=1.0 0
-
-def f(x):
-    theta=0.25
-    # --- replace with your definition of indicatorFunction:
-    if ((abs(x[0]) < theta/2) and x[2]<-1/2+theta):
-        return 1    #Phase1
-    elif ((abs(x[1]) < theta/2) and x[2]>1/2-theta):
-        return 1    #Phase1
-    else :
-        return 0    #Phase2
-

outputs/Results/2x2/2/BMatrix.txt

-1 1 -0.0144287202371896055
-1 2 0.014428690223795131
-1 3 -2.0320334857961725e-11

outputs/Results/2x2/2/QMatrix.txt

-1 1 12.4813265250416485
-1 2 2.0831795458888962
-1 3 3.83811740666402797e-10
-2 1 2.0831795513179836
-2 2 12.481326525238547
-2 3 1.78942192122814843e-10
-3 1 -6.40561063029264521e-10
-3 2 -1.00945173229121783e-09
-3 3 10.3922965257922275

outputs/Results/2x2/2/output.txt

-material_prestrain used: material_neukamm
------ Input Parameters -----: 
-alpha: 8
-gamma: 1
-theta: 0.25
-beta: 3
-material parameters: 
-mu1: 1
-mu2: 3
-lambda1: 1
-lambda2: 3
-----------------------------: 
-Number of Elements in each direction: [4,4,4]
-size of FiniteElementBasis: 240
-Solver-type used:  CG-Solver
----------- OUTPUT ----------
- --------------------
-Corrector-Matrix M_1: 
--0.00922212 -2.49142e-10 0
--2.49142e-10 -0.00494017 0
-0 0 0
-
- --------------------
-Corrector-Matrix M_2: 
-0.00494016 -2.03726e-10 0
--2.03726e-10 0.00922211 0
-0 0 0
-
- --------------------
-Corrector-Matrix M_3: 
-1.90697e-10 -4.74653e-11 0
--4.74653e-11 1.93256e-10 0
-0 0 0
-
- --------------------
-Effective Matrix Q: 
-12.4813 2.08318 3.83812e-10
-2.08318 12.4813 1.78942e-10
--6.40561e-10 -1.00945e-09 10.3923
-
---- Prestrain Output --- 
-B_hat: -0.150032 0.150032 -2.16498e-10
-B_eff: -0.0144287 0.0144287 -2.03203e-11 (Effective Prestrain)
------------------------- 
-q1=12.4813
-q2=12.4813
-q3=10.3923
-q12=2.08318
-b1=-0.0144287
-b2=0.0144287
-b3=-2.03203e-11
-b1_hat: -0.150032
-b2_hat: 0.150032
-b3_hat: -2.16498e-10
-mu_gamma=10.3923
-q_onetwo=2.083180
----------------------------------------------------------------------------------------------------------
-  Levels     |      q1      |      q2      |      q3      |      b1      |      b2      |      b3      | 
----------------------------------------------------------------------------------------------------------
-     2       & 1.24813e+01  & 1.24813e+01  & 1.03923e+01  & -1.44287e-02 & 1.44287e-02  & -2.03203e-11 & 
----------------------------------------------------------------------------------------------------------

outputs/Results/2x2/2/parameter

-mu=80 60
-lambda=80 25
-rho=1.0 0
-
-def f(x):
-    theta=0.25
-    factor=0.8
-    # --- replace with your definition of indicatorFunction:
-    if ((abs(x[0]) < theta/2) and x[2]<-1/2+theta):
-        return 1    #Phase1
-    elif ((abs(x[1]) < factor*theta/2) and x[2]>1/2-theta):
-        return 1    #Phase1
-    else :
-        return 0    #Phase2
-

outputs/Results/2x2/3/BMatrix.txt

-1 1 -0.0144287202371896055
-1 2 0.014428690223795131
-1 3 -2.0320334857961725e-11

outputs/Results/2x2/3/QMatrix.txt

-1 1 12.4813265250416485
-1 2 2.0831795458888962
-1 3 3.83811740666402797e-10
-2 1 2.0831795513179836
-2 2 12.481326525238547
-2 3 1.78942192122814843e-10
-3 1 -6.40561063029264521e-10
-3 2 -1.00945173229121783e-09
-3 3 10.3922965257922275

outputs/Results/2x2/3/output.txt

-material_prestrain used: material_neukamm
------ Input Parameters -----: 
-alpha: 8
-gamma: 1
-theta: 0.25
-beta: 3
-material parameters: 
-mu1: 1
-mu2: 3
-lambda1: 1
-lambda2: 3
-----------------------------: 
-Number of Elements in each direction: [4,4,4]
-size of FiniteElementBasis: 240
-Solver-type used:  CG-Solver
----------- OUTPUT ----------
- --------------------
-Corrector-Matrix M_1: 
--0.00922212 -2.49142e-10 0
--2.49142e-10 -0.00494017 0
-0 0 0
-
- --------------------
-Corrector-Matrix M_2: 
-0.00494016 -2.03726e-10 0
--2.03726e-10 0.00922211 0
-0 0 0
-
- --------------------
-Corrector-Matrix M_3: 
-1.90697e-10 -4.74653e-11 0
--4.74653e-11 1.93256e-10 0
-0 0 0
-
- --------------------
-Effective Matrix Q: 
-12.4813 2.08318 3.83812e-10
-2.08318 12.4813 1.78942e-10
--6.40561e-10 -1.00945e-09 10.3923
-
---- Prestrain Output --- 
-B_hat: -0.150032 0.150032 -2.16498e-10
-B_eff: -0.0144287 0.0144287 -2.03203e-11 (Effective Prestrain)
------------------------- 
-q1=12.4813
-q2=12.4813
-q3=10.3923
-q12=2.08318
-b1=-0.0144287
-b2=0.0144287
-b3=-2.03203e-11
-b1_hat: -0.150032
-b2_hat: 0.150032
-b3_hat: -2.16498e-10
-mu_gamma=10.3923
-q_onetwo=2.083180
----------------------------------------------------------------------------------------------------------
-  Levels     |      q1      |      q2      |      q3      |      b1      |      b2      |      b3      | 
----------------------------------------------------------------------------------------------------------
-     2       & 1.24813e+01  & 1.24813e+01  & 1.03923e+01  & -1.44287e-02 & 1.44287e-02  & -2.03203e-11 & 
----------------------------------------------------------------------------------------------------------

outputs/Results/2x2/3/parameter

-mu=80 60
-lambda=80 25
-rho=1.0 0
-
-#Indicator function that determines both phases
-# x[0] : y1-component
-# x[1] : y2-component
-# x[2] : x3-component
-def f(x):
-    theta=0.5
-    factor=1
-    # --- replace with your definition of indicatorFunction:
-    if ((abs(x[0]) < theta/2) and x[2]<-1/2+theta):
-        return 1    #Phase1
-    elif ((abs(x[1]) < factor*theta/2) and x[2]>1/2-factor*theta):
-        return 1    #Phase1
-    else :
-        return 0    #Phase2
-

outputs/Results/2x2/4/BMatrix.txt

-1 1 -2.74206114740746898
-1 2 -3.10270304357049476
-1 3 3.0398257121603124e-09

outputs/Results/2x2/4/QMatrix.txt

-1 1 14.4975598095779166
-1 2 3.98327256772200711
-1 3 -6.56940620062560124e-09
-2 1 3.98327260564529295
-2 2 14.8275228823281235
-2 3 -6.51275988527633122e-09
-3 1 -2.09762192341099332e-09
-3 2 -1.83361834006557389e-09
-3 3 10.6556275724235974

outputs/Results/2x2/4/output.txt

-material_prestrain used: material_neukamm
------ Input Parameters -----: 
-alpha: 8
-gamma: 1
-theta: 0.25
-beta: 3
-material parameters: 
-mu1: 1
-mu2: 3
-lambda1: 1
-lambda2: 3
-----------------------------: 
-Number of Elements in each direction: [4,4,4]
-size of FiniteElementBasis: 240
-Solver-type used:  CG-Solver
----------- OUTPUT ----------
- --------------------
-Corrector-Matrix M_1: 
--0.277147 1.10701e-09 0
-1.10701e-09 -0.334813 0
-0 0 0
-
- --------------------
-Corrector-Matrix M_2: 
--0.478197 8.17764e-10 0
-8.17764e-10 -0.188937 0
-0 0 0
-
- --------------------
-Corrector-Matrix M_3: 
-2.06316e-11 0.115943 0
-0.115943 3.27431e-10 0
-0 0 0
-
- --------------------
-Effective Matrix Q: 
-14.4976 3.98327 -6.56941e-09
-3.98327 14.8275 -6.51276e-09
--2.09762e-09 -1.83362e-09 10.6556
-

outputs/Results/2x2/4/parameter

-mu=80 60
-lambda=80 25
-rho=1.0 0
-
-#Indicator function that determines both phases
-# x[0] : y1-component
-# x[1] : y2-component
-# x[2] : x3-component
-def f(x):
-    theta=0.5
-    factor=0.5
-    # --- replace with your definition of indicatorFunction:
-    if ((abs(x[0]) < theta/2) and x[2]<-1/2+theta):
-        return 1    #Phase1
-    elif ((abs(x[1]) < factor*theta/2) and x[2]>1/2-factor*theta):
-        return 1    #Phase1
-    else :
-        return 0    #Phase2
-

outputs/Results/4x4/1/BMatrix.txt

-1 1 0.660848323245719738
-1 2 -0.66084827993852846
-1 3 -6.068513378040036e-10