Merge remote-tracking branch 'origin/development' into dislotwin-fix-TWIP-TRIP

PRIVATE

@@ -1 +1 @@
-Subproject commit 5774122bf48d637704bb4afb10b87c34a4dbcaba
+Subproject commit 0d639a9ba41db279b0d2825c8e8eddf0ccd91326

examples/config/phase/AISI304.yaml

@@ -1,6 +1,7 @@
 references:
-  - H.M. Ledbetter
+  - H.M. Ledbetter,
-    physica status solidi (a) 85(1):89-96, 1984
+    physica status solidi (a) 85(1):89-96, 1984,
     https://doi.org/10.1002/pssa.2210850111
+
 lattice: cF
 rho: 7937.0

examples/config/phase/Ag.yaml

@@ -1,4 +1,5 @@
 references:
   - https://en.wikipedia.org/wiki/Silver
+
 lattice: cF
 rho: 10490.0

examples/config/phase/Al.yaml

@@ -1,4 +1,5 @@
 references:
   - https://en.wikipedia.org/wiki/Aluminium
+
 lattice: cF
 rho: 2700.0

examples/config/phase/Au.yaml

@@ -1,4 +1,5 @@
 references:
   - https://en.wikipedia.org/wiki/Gold
+
 lattice: cF
 rho: 19300.0

examples/config/phase/Cu.yaml

@@ -1,4 +1,5 @@
 references:
   - https://en.wikipedia.org/wiki/Copper
+
 lattice: cF
 rho: 8960.0

examples/config/phase/Fe.yaml

@@ -1,4 +1,5 @@
 references:
   - https://en.wikipedia.org/wiki/Iron
+
 lattice: cI
 rho: 7874.0

examples/config/phase/Mg.yaml

@@ -2,6 +2,7 @@ references:
   - D. Tromans,
     International Journal of Recent Research and Applied Studies 6(4):462-483, 2011,
     https://www.arpapress.com/Volumes/Vol6Issue4/IJRRAS_6_4_14.pdf
+
 lattice: hP
 c/a: 1.62350
 rho: 1740.0

examples/config/phase/Ni.yaml

@@ -1,4 +1,5 @@
 references:
   - https://en.wikipedia.org/wiki/Nickel
+
 lattice: cF
 rho: 8908.0

examples/config/phase/Sn-beta.yaml

@@ -0,0 +1,9 @@
+references:
+  - J.A. Rayne and B.S. Chandrasekhar,
+    Physical Review 120(5):1658-1663, 1960,
+    https://doi.org/10.1103/PhysRev.120.1658
+  - https://en.wikipedia.org/wiki/Tin
+
+lattice: tI
+c/a: 0.5458 # T=300K (c=31.83nm, a=5.832nm)
+rho: 7265.0

examples/config/phase/Ti.yaml

@@ -1,6 +1,7 @@
 references:
   - https://www.totalmateria.com/page.aspx?ID=CheckArticle&site=ktn&NM=221
   - https://en.wikipedia.org/wiki/Titanium
+
 lattice: hP
 c/a: 1.587
 rho: 4506.0

examples/config/phase/W.yaml

@@ -1,4 +1,5 @@
 references:
   - https://en.wikipedia.org/wiki/Tungsten
+
 lattice: cI
 rho: 19300.0

examples/config/phase/bct.yaml

@@ -0,0 +1,2 @@
+lattice: tI
+c/a: 0.55

examples/config/phase/damage/anisobrittle_cubic.yaml

@@ -1,11 +1,13 @@
 type: anisobrittle
+
+output: [f_phi]
+
 N_cl: [3]
+
 g_crit: [0.5e+7]
 s_crit: [0.006666]
 dot_o: 1.e-3
 q: 20
 
-output: [f_phi]
-
 D_11: 1.0
 mu: 0.001

examples/config/phase/damage/isobrittle_generic.yaml

@@ -1,7 +1,7 @@
 type: isobrittle
-W_crit: 1400000.0
 
 output: [f_phi]
 
+W_crit: 1400000.0
 D_11: 1.0
 mu: 0.001

examples/config/phase/hcp.yaml

@@ -0,0 +1,2 @@
+lattice: hP
+c/a: 1.6333

examples/config/phase/mechanical/eigen/thermalexpansion_AISI304.yaml

@@ -1,7 +1,13 @@
 type: thermalexpansion
+
 references:
-  - R.H. Bogaard et al.
+  - R.H. Bogaard et al.,
-    Thermochimica Acta 218:373-393, 1993
+    Thermochimica Acta 218:373-393, 1993,
-    https://doi.org/10.1016/0040-6031(93)80437-F
+    https://doi.org/10.1016/0040-6031(93)80437-F,
-A_11: 15.0e-6
+    fit to Fig. 6 (T_min=100K, T_max=1400K)
-T_ref: 300.0
+
+A_11: 2.068e-08
+A_11,T: 1.579e-09
+A_11,T^2: 3.449e-13
+
+T_ref: 293.15

examples/config/phase/mechanical/eigen/thermalexpansion_Al.yaml

@@ -1,5 +1,7 @@
 type: thermalexpansion
+
 references:
-  - https://en.wikipedia.org/wiki/Thermal_expansion
+  - https://en.wikipedia.org/wiki/Thermal_expansion,
+    293.15K
+
 A_11: 23.1e-6
-T_ref: 293.15

examples/config/phase/mechanical/eigen/thermalexpansion_Au.yaml

@@ -1,5 +1,7 @@
 type: thermalexpansion
+
 references:
-  - https://en.wikipedia.org/wiki/Thermal_expansion
+  - https://en.wikipedia.org/wiki/Thermal_expansion,
+    293.15K
+
 A_11: 14.e-6
-T_ref: 293.15

examples/config/phase/mechanical/eigen/thermalexpansion_C35E.yaml

@@ -1,8 +1,11 @@
 type: thermalexpansion
+
 references:
   - https://commons.wikimedia.org/wiki/File:Coefficient_dilatation_lineique_aciers.svg,
-    fitted from image description (Scilab code)
+    fit to image description (Scilab code)
+
 A_11: 12.70371e-6
 A_11,T: 7.54e-9
 A_11,T^2: -1.0e-11
+
 T_ref: 273.0

examples/config/phase/mechanical/eigen/thermalexpansion_Cu.yaml

@@ -1,5 +1,7 @@
 type: thermalexpansion
+
 references:
-  - https://en.wikipedia.org/wiki/Thermal_expansion
+  - https://en.wikipedia.org/wiki/Thermal_expansion,
+    293.15K
+
 A_11: 17.e-6
-T_ref: 293.15

examples/config/phase/mechanical/eigen/thermalexpansion_Fe.yaml

@@ -1,5 +1,7 @@
 type: thermalexpansion
+
 references:
-  - https://en.wikipedia.org/wiki/Thermal_expansion
+  - https://en.wikipedia.org/wiki/Thermal_expansion,
+    293.15K
+
 A_11: 11.8e-6
-T_ref: 293.15

examples/config/phase/mechanical/eigen/thermalexpansion_Sn-beta.yaml

@@ -0,0 +1,17 @@
+type: thermalexpansion
+
+references:
+  - V.T. Deshpande and D.B. Sirdeshmukh,
+    Acta Crystallographica 15:294-295, 1962,
+    https://doi.org/10.1107/S0365110X62000742,
+    fit to Tab. 2 (T_min=30ºC, T_max=210ºC)
+
+A_11: 1.639e-05
+A_11,T: 1.799e-08
+A_11,T^2: 1.734e-10
+
+A_33: 3.263e-05
+A_33,T: 1.387e-08
+A_33,T^2: 5.794e-10
+
+T_ref: 293.15

examples/config/phase/mechanical/eigen/thermalexpansion_W.yaml

@@ -1,5 +1,7 @@
 type: thermalexpansion
+
 references:
-  - https://en.wikipedia.org/wiki/Thermal_expansion
+  - https://en.wikipedia.org/wiki/Thermal_expansion,
+    293.15K
+
 A_11: 4.5e-6
-T_ref: 293.15

examples/config/phase/mechanical/eigen/thermalexpansion_X20Cr13.yaml

@@ -1,7 +1,10 @@
 type: thermalexpansion
+
 references:
-  - https://commons.wikimedia.org/wiki/File:Coefficient_dilatation_lineique_aciers.svg
+  - https://commons.wikimedia.org/wiki/File:Coefficient_dilatation_lineique_aciers.svg,
-    fitted from image description (Scilab code)
+    fit to image description (Scilab code)
+
 A_11: 11.365e-6
 A_11,T: 5.0e-9
+
 T_ref: 273.0

examples/config/phase/mechanical/elastic/Hooke_AISI304.yaml

@@ -1,8 +1,10 @@
 type: Hooke
+
 references:
-  - H.M. Ledbetter
+  - H.M. Ledbetter,
-    physica status solidi (a) 85(1):89-96, 1984
+    physica status solidi (a) 85(1):89-96, 1984,
     https://doi.org/10.1002/pssa.2210850111
+
 C_11: 204.6e+9
 C_12: 137.7e+9
 C_44: 126.2e+9

examples/config/phase/mechanical/elastic/Hooke_Ag.yaml

@@ -1,4 +1,5 @@
 type: Hooke
+
 references:
   - J.R. Neighbours and G.A. Alers,
     Physical Review 111:707-712, 1958,
@@ -7,8 +8,6 @@ references:
     Journal of Applied Physics 37:3567-3572, 1966,
     https://doi.org/10.1063/1.1708903
 
-T_ref: 300
-
 C_11: 122.9e+9
 C_11,T: -313.5e+5
 C_11,T^2: -107.3e+2
@@ -20,3 +19,5 @@ C_12,T^2: -681.6e+1
 C_44: 42.63e+9
 C_44,T: -180.5e+5
 C_44,T^2: -353.8e+1
+
+T_ref: 300

examples/config/phase/mechanical/elastic/Hooke_Al.yaml

@@ -1,22 +1,25 @@
 type: Hooke
+
 references:
   - G.N. Kamm and G.A. Alers,
     Journal of Applied Physics 35:327-330, 1964,
-    https://doi.org/10.1063/1.1713309
+    https://doi.org/10.1063/1.1713309,
+    fit to Tab. I (T_min=100K, T_max=300K)
   - D. Gerlich and E.S. Fisher,
     Journal of Physics and Chemistry of Solids 30:1197-1205, 1969
-    https://doi.org/10.1016/0022-3697(69)90377-1
+    https://doi.org/10.1016/0022-3697(69)90377-1,
+    fit to Tab. 2 (T_min=293K, T_max=900K)
 
-T_ref: 300
+C_11: 106.9e+9
+C_11,T: -3.685e+7
+C_11,T^2: -1.016e+4
 
-C_11: 106.1e+9
+C_12: 60.55e+9
-C_11,T: -359.3e+5
+C_12,T: -8.307e+6
-C_11,T^2: -152.7e+2
+C_12,T^2: -4.353e+3
 
-C_12: 57.83e+9
+C_44: 28.37e+9
-C_12,T: -781.6e+4
+C_44,T: -1.418e+7
-C_12,T^2: -551.3e+1
+C_44,T^2: -3.253e+3
 
-C_44: 24.31e+9
+T_ref: 293.15
-C_44,T: -142.9e+5
-C_44,T^2: -404.6e+1

examples/config/phase/mechanical/elastic/Hooke_Au.yaml

@@ -1,9 +1,11 @@
 type: Hooke
+
 references:
   - J.P. Hirth and J. Lothe,
     Theory of Dislocations, 1982,
     John Wiley & Sons,
     page 837
+
 C_11: 186.e+9
 C_12: 157.e+9
 C_44: 42.e+9

examples/config/phase/mechanical/elastic/Hooke_Cu.yaml

@@ -1,7 +1,21 @@
 type: Hooke
+
 references:
-  - https://www.mit.edu/~6.777/matprops/copper.htm,
+  - W.C. Overton, Jr. and J. Gaffney,
-    fixed typo
+    Physical Review 98(4):969-977, 1955,
-C_11: 168.3e+9
+    https://doi.org/10.1103/PhysRev.98.969,
-C_12: 122.1e+9
+    fit to Tab. I (T_min=100K, T_max=300K)
-C_44: 75.7e+9
+
+C_11: 168.6e+9
+C_11,T: -3.779e+7
+C_11,T^2: -2.536e+4
+
+C_12: 121.5e+9
+C_12,T: -1.632e+7
+C_12,T^2: -1.116e+4
+
+C_44: 75.59e+9
+C_44,T: -2.912e+7
+C_44,T^2: -1.669e+4
+
+T_ref: 293.15

examples/config/phase/mechanical/elastic/Hooke_Fe.yaml

@@ -1,19 +1,21 @@
 type: Hooke
+
 references:
   - D.J. Dever,
     Journal of Applied Physics 43(8):3293-3301, 1972,
     https://doi.org/10.1063/1.1661710
+    fit to Tab. II (T_min=25ºC, T_max=880ºC)
 
-T_ref: 300
+C_11: 232.1e+9
+C_11,T: -4.678e+7
+C_11,T^2: -5.762e+4
 
-C_11: 231.7e+9
+C_12: 135.9e+9
-C_11,T: -47.6e+6
+C_12,T: -1.695e+7
-C_11,T^2: -54.4e+3
+C_12,T^2: 1.555e+3
 
-C_12: 135.8e+9
+C_44: 117.0e+9
-C_12,T: -12.9e+6
+C_44,T: -2.047e+7
-C_12,T^2: -7.3e+3
+C_44,T^2: -2.814e+2
 
-C_44: 116.8e+9
+T_ref: 293.15
-C_44,T: -19.4e+6
-C_44,T^2: -2.5e+3

examples/config/phase/mechanical/elastic/Hooke_Mg.yaml

@@ -1,10 +1,29 @@
 type: Hooke
+
 references:
-  - D. Tromans,
+  - L.J. Slutsky and C.W. Garland,
-    International Journal of Recent Research and Applied Studies 6(4):462-483, 2011,
+    Physical Review 107(4):972-976, 1957,
-    https://www.arpapress.com/Volumes/Vol6Issue4/IJRRAS_6_4_14.pdf
+    https://doi.org/10.1103/PhysRev.107.972,
-C_11: 59.3e+9
+    fit to Tab. I (T_min=100K, T_max=300K)
-C_33: 61.5e+9
+
-C_44: 16.4e+9
+C_11: 59.50e+9
-C_12: 25.7e+9
+C_11,T: -1.930e+7
-C_13: 21.4e+9
+C_11,T^2: -1.215e+4
+
+C_33: 61.72e+9
+C_33,T: -2.175e+7
+C_33,T^2: -5.624e+3
+
+C_44: 16.46e+9
+C_44,T: -1.006e+7
+C_44,T^2: -7.692e+3
+
+C_12: 25.62e+9
+C_12,T: -2.216e+6
+C_12,T^2: -4.138e+3
+
+C_13: 21.46e+9
+C_13,T: -1.921e+6
+C_13,T^2: -4.283e+3
+
+T_ref: 293.15

examples/config/phase/mechanical/elastic/Hooke_Ni.yaml

@@ -1,9 +1,21 @@
 type: Hooke
+
 references:
-  - J.P. Hirth and J. Lothe,
+  - G.A. Alers,
-    Theory of Dislocations, 1982,
+    Journal of Physics and Chemistry of Solids 13(1-2):40-55, 1960,
-    John Wiley & Sons,
+    https://doi.org/10.1016/0022-3697(60)90125-6,
-    page 837
+    fit to Tab. 2 (T_min=100K, T_max=700K)
-C_11: 246.5e+9
+
-C_12: 147.3e+9
+C_11: 251.0e+9
-C_44: 124.7e+9
+C_11,T: -4.998e+7
+C_11,T^2: -2.952e+4
+
+C_12: 150.0e+9
+C_12,T: -4.269e+6
+C_12,T^2: -6.429e+3
+
+C_44: 123.7e+9
+C_44,T: -3.618e+7
+C_44,T^2: -7.024e+3
+
+T_ref: 293.15

examples/config/phase/mechanical/elastic/Hooke_SAE1050-martensite.yaml

@@ -1,8 +1,10 @@
 type: Hooke
+
 references:
   - S.A. Kim and W.L. Johnson,
     Materials Science & Engineering A 452-453:633-639, 2007,
     https://doi.org/10.1016/j.msea.2006.11.147
+
 C_11: 268.1e+9
 C_12: 111.2e+9
 C_44: 79.06e+9

examples/config/phase/mechanical/elastic/Hooke_Sn-beta.yaml

@@ -0,0 +1,32 @@
+type: Hooke
+
+references:
+  - J.A. Rayne and B.S. Chandrasekhar,
+    Physical Review 120(5):1658-1663, 1960,
+    https://doi.org/10.1103/PhysRev.120.1658,
+    fit to Fig. 2 (T_min=100K, T_max=300K) and Tab. IV (C_13, T_min=77K, T_max=300K)
+
+C_11: 72.90e+9
+C_11,T: -4.399e+7
+C_11,T^2: -2.645e+4
+
+C_12: 59.27e+9
+C_12,T: 1.058e+7
+C_12,T^2: 1.002e+4
+
+C_13: 35.80e+9
+C_13,T: -2.870e+6
+
+C_33: 88.78e+9
+C_33,T: -5.250e+7
+C_33,T^2: 3.546e+3
+
+C_44: 22.26e+9
+C_44,T: -1.982e+7
+C_44,T^2: -8.711e+3
+
+C_66: 24.18e+9
+C_66,T: -1.806e+7
+C_66,T^2: -4.112e+3
+
+T_ref: 293.15

examples/config/phase/mechanical/elastic/Hooke_TWIP-steel.yaml

@@ -1,4 +1,5 @@
 type: Hooke
+
 references:
   - D. Music et al.,
     Applied Physics Letters 99(19):191904, 2007,
@@ -6,6 +7,7 @@ references:
   - S.L. Wong et al.,
     Acta Materialia 118:140-151, 2016,
     https://doi.org/10.1016/j.actamat.2016.07.032
+
 C_11: 175.0e+9
 C_12: 115.0e+9
 C_44: 135.0e+9

examples/config/phase/mechanical/elastic/Hooke_Ti.yaml

@@ -1,10 +1,33 @@
 type: Hooke
+
 references:
-  - L. Wang et al.,
+  - E.S. Fisher and C.J. Renken,
-    Acta Materialia 132:598-610, 2017,
+    Physical Review 135(2A):A482-A494, 1964,
-    https://doi.org/10.1016/j.actamat.2017.05.015
+    https://doi.org/10.1103/PhysRev.135.A482,
-C_11: 162.4e+9
+    fit to Tab. IV (T_min=150K, T_max=250K)
-C_33: 181.6e+9
+  - H. Ogi et al.,
-C_44: 47.2e+9
+    Acta Materialia 52(7):2075-2080, 2004,
-C_12: 92.e+9
+    https://doi.org/10.1016/j.actamat.2004.01.002,
-C_13: 69.e+9
+    fit to Fig. 3 (T_min=300K, T_max=900K)
+
+C_11: 162.6e+9
+C_11,T: -6.150e+7
+C_11,T^2: -5.557e+2
+
+C_33: 183.3e+9
+C_33,T: -1.655e+07
+C_33,T^2: -1.022e+04
+
+C_44: 45.80e+9
+C_44,T: -2.936e+07
+C_44,T^2: 7.120e+02
+
+C_12: 89.97e+9
+C_12,T: 2.776e+6
+C_12,T^2: -2.389e+4
+
+C_13: 69.53e+9
+C_13,T: 1.057e+7
+C_13,T^2: -2.966e+3
+
+T_ref: 293.15

examples/config/phase/mechanical/elastic/Hooke_W.yaml

@@ -1,8 +1,20 @@
 type: Hooke
+
 references:
-  - D. Cereceda et al.,
+  - F.H. Featherston and J.R. Nieghbours,
-    International Journal of Plasticity 78:242-265, 2016,
+    Physical Review 130(4):1324-1333,
-    https://doi.org/10.1016/j.ijplas.2015.09.002
+    https://doi.org/10.1103/PhysRev.130.1324,
-C_11: 523.e+9
+    fit to Tab. III (T_min=100K, T_max=300K)
-C_12: 202.e+9
+
-C_44: 161.e+9
+C_11: 523.6e+9
+C_11,T: -7.607e+7
+C_11,T^2: -1.551e+5
+
+C_12: 205.1e+9
+C_12,T: -2.843e+6
+
+C_44: 160.8e+9
+C_44,T: -1.057e+7
+C_44,T^2: 9.933e+3
+
+T_ref: 293.15

examples/config/phase/mechanical/elastic/Hooke_vanishing-Poisson-ratio.yaml

@@ -1,8 +1,10 @@
 type: Hooke
+
 references:
   - T. Maiti and P. Eisenlohr,
     Scripta Materialia 145:37-40, 2018,
     https://doi.org/10.1016/j.scriptamat.2017.09.047
+
 C_11: 1.e+8
 C_12: 1.e+6
 C_44: 4.95e+7

examples/config/phase/mechanical/plastic/dislotungsten_W.yaml

@@ -1,4 +1,5 @@
 type: dislotungsten
+
 references:
   - D. Cereceda et al.,
     International Journal of Plasticity 78:242-265, 2016,
@@ -6,8 +7,11 @@ references:
   - R. Gröger et al.,
     Acta Materialia 56(19):5412-5425, 2008,
     https://doi.org/10.1016/j.actamat.2008.07.037
+
 output: [Lambda_sl]
+
 N_sl: [12]
+
 b_sl: [2.72e-10]
 rho_mob_0: [1.0e+9]  # estimated from section 3.2
 rho_dip_0: [1.0]     # not given

examples/config/phase/mechanical/plastic/dislotwin_IF-steel.yaml

@@ -1,4 +1,5 @@
 type: dislotwin
+
 references:
   - K. Sedighiani et al.,
     International Journal of Plasticity 134:102779, 2020,
@@ -6,8 +7,11 @@ references:
   - K. Sedighiani et al.,
     Mechanics of Materials, 164:104117, 2022,
     https://doi.org/10.1016/j.mechmat.2021.104117
+
 output: [rho_dip, rho_mob]
+
 N_sl: [12, 12]
+
 b_sl: [2.49e-10, 2.49e-10]
 rho_mob_0: [2.81e12, 2.8e+12]
 rho_dip_0: [1.0, 1.0] # not given

examples/config/phase/mechanical/plastic/isotropic_free-surface.yaml

@@ -1,9 +1,12 @@
 type: isotropic
+
 references:
   - T. Maiti and P. Eisenlohr,
     Scripta Materialia 145:37-40, 2018,
     https://doi.org/10.1016/j.scriptamat.2017.09.047
+
 output: [xi]
+
 dot_gamma_0: 0.001
 n: 20.
 xi_0: 0.3e+6

examples/config/phase/mechanical/plastic/nonlocal_Al.yaml

@@ -1,10 +1,13 @@
 type: nonlocal
+
 references:
-  C. Kords,
+  - C. Kords,
-  On the role of dislocation transport in the constitutive description of crystal plasticity,
+    On the role of dislocation transport in the constitutive description of crystal plasticity,
-  RWTH Aachen 2013,
+    RWTH Aachen 2013,
-  http://publications.rwth-aachen.de/record/229993/files/4862.pdf
+    http://publications.rwth-aachen.de/record/229993/files/4862.pdf
+
 output: [rho_u_ed_pos, rho_b_ed_pos, rho_u_ed_neg, rho_b_ed_neg, rho_u_sc_pos, rho_b_sc_pos, rho_u_sc_neg, rho_b_sc_neg, rho_d_ed, rho_d_sc]
+
 N_sl: [12]
 
 b_sl: [2.86e-10]

examples/config/phase/mechanical/plastic/nonlocal_Ni.yaml

@@ -1,10 +1,13 @@
 type: nonlocal
+
 references:
-  C. Kords,
+  - C. Kords,
-  On the role of dislocation transport in the constitutive description of crystal plasticity,
+    On the role of dislocation transport in the constitutive description of crystal plasticity,
-  RWTH Aachen 2013,
+    RWTH Aachen 2013,
-  http://publications.rwth-aachen.de/record/229993/files/4862.pdf
+    http://publications.rwth-aachen.de/record/229993/files/4862.pdf
+
 output: [rho_u_ed_pos, rho_b_ed_pos, rho_u_ed_neg, rho_b_ed_neg, rho_u_sc_pos, rho_b_sc_pos, rho_u_sc_neg, rho_b_sc_neg, rho_d_ed, rho_d_sc]
+
 N_sl: [12]
 
 b_sl: [2.48e-10]

examples/config/phase/mechanical/plastic/phenopowerlaw_Al.yaml

@@ -1,4 +1,5 @@
 type: phenopowerlaw
+
 references:
   - W.F. Hosford et al.,
     Acta Metallurgica 8(3):187-199, 1960,
@@ -7,8 +8,11 @@ references:
   - U.F. Kocks,
     Metallurgical and Materials Transactions B 1:1121–1143, 1970,
     https://doi.org/10.1007/BF02900224
+
 output: [xi_sl, gamma_sl]
+
 N_sl: [12]
+
 n_sl: 20
 a_sl: 3.1
 h_0_sl-sl: 1.7e+8

examples/config/phase/mechanical/plastic/phenopowerlaw_Au.yaml

@@ -1,4 +1,5 @@
 type: phenopowerlaw
+
 references:
   - D. Ma et al.,
     Acta Materialia 103:796-808, 2016,
@@ -9,8 +10,11 @@ references:
   - U.F. Kocks,
     Metallurgical and Materials Transactions B 1:1121–1143, 1970,
     https://doi.org/10.1007/BF02900224
+
 output: [xi_sl, gamma_sl]
+
 N_sl: [12]
+
 n_sl: 83.3
 a_sl: 1.0
 h_0_sl-sl: 75.0e+6

examples/config/phase/mechanical/plastic/phenopowerlaw_Cu.yaml

@@ -1,4 +1,5 @@
 type: phenopowerlaw
+
 references:
   - T Takeuchi,
     Transactions of the Japan Institute of Metals 16(10):629-640, 1975,
@@ -7,8 +8,11 @@ references:
   - U.F. Kocks,
     Metallurgical and Materials Transactions B 1:1121–1143, 1970,
     https://doi.org/10.1007/BF02900224
+
 output: [xi_sl, gamma_sl]
+
 N_sl: [12]
+
 n_sl: 20
 a_sl: 1.0
 h_0_sl-sl: 2.4e+8

examples/config/phase/mechanical/plastic/phenopowerlaw_DP-steel-ferrite.yaml

@@ -1,4 +1,5 @@
 type: phenopowerlaw
+
 references:
   - C.C. Tasan et al.,
     Acta Materialia 81:386-400, 2014,
@@ -6,8 +7,11 @@ references:
   - U.F. Kocks,
     Metallurgical and Materials Transactions B 1:1121–1143, 1970,
     https://doi.org/10.1007/BF02900224
+
 output: [xi_sl, gamma_sl]
+
 N_sl: [12, 12]
+
 n_sl: 20
 a_sl: 2.25
 h_0_sl-sl: 1.0e+9

examples/config/phase/mechanical/plastic/phenopowerlaw_Mg.yaml

@@ -1,4 +1,5 @@
 type: phenopowerlaw
+
 references:
   - F. Wang et al.,
     Acta Materialia 80:77-93, 2014,

examples/config/phase/mechanical/plastic/phenopowerlaw_Sn-beta.yaml

@@ -0,0 +1,33 @@
+type: phenopowerlaw
+
+references:
+  - A. Chakraborty and P. Eisenlohr,
+    Journal of Applied Physics 124:025302, 2018,
+    https://doi.org/10.1063/1.5029933
+
+output: [xi_sl, gamma_sl]
+
+N_sl: [2, 2, 2, 4, 2, 4, 2, 2, 4, 0, 0, 8]
+
+n_sl: 6.0
+a_sl: 2.0
+h_0_sl-sl: 20.0e+6
+xi_0_sl:   [8.5e+6,  4.3e+6, 10.4e+6, 4.5e+6,  5.6e+6,  5.1e+6,  7.4e+6, 15.0e+6, 6.6e+6, 0.0, 0.0, 12.0e+6]
+xi_inf_sl: [11.0e+6, 9.0e+6, 11.0e+6, 9.0e+6, 10.0e+6, 10.0e+6, 10.0e+6, 10.0e+6, 9.0e+6, 0.0, 0.0, 13.0e+6]
+h_sl-sl: [+1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,
+          +1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,
+          +1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,
+          +1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,
+          +1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0, # 50
+          +1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,
+          +1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,
+          +1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,
+          +1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,
+          -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, # 100
+          -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0,
+          -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0,
+          -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0,
+          -1.0, -1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,
+          +1.0, -1.0, -1.0,  1.0,  1.0, -1.0, -1.0,  1.0,  1.0,  1.0, # 150
+          +1.0,  1.0,  1.0,  1.0,  1.0,  1.0]  # unused entries are indicated by -1.0
+dot_gamma_0_sl: 2.6e-8

examples/config/phase/mechanical/plastic/phenopowerlaw_Ti.yaml

@@ -1,4 +1,5 @@
 type: phenopowerlaw
+
 references:
   - C. Zambaldi et al.,
     Journal of Materials Research 27(1):356-367, 2021,
@@ -6,9 +7,11 @@ references:
   - L. Wang et al.,
     Acta Materialia 132:598-610, 2017,
     https://doi.org/10.1016/j.actamat.2017.05.015
+
 output: [gamma_sl]
 
 N_sl: [3, 3, 0, 12]   # basal, prism, -, 1. pyr<c+a>
+
 n_sl: 20
 a_sl: 2.0
 dot_gamma_0_sl: 0.001

examples/config/phase/thermal/AISI304.yaml

@@ -5,5 +5,6 @@ references:
   - R.H. Bogaard et al.
     Thermochimica Acta 218:373-393, 1993
     https://doi.org/10.1016/0040-6031(93)80437-F
+
 C_p: 470.0
 K_11: 14.34

examples/config/phase/thermal/Al.yaml

@@ -1,5 +1,14 @@
 references:
-  - https://www.engineeringtoolbox.com/thermal-conductivity-metals-d_858.html
+  - J.G. Hust and A.B. Lankford,
+    Thermal Conductivity of Aluminum, Copper, Iron, and Tungsten from 1K to the Melting Point,
+    US Department of Commerce, Boulder, Colorado, 1984,
+    fit to Tab. 3.4.1 (RRR=1000, T_min=200K, T_max=900K)
   - https://www.engineeringtoolbox.com/specific-heat-metals-d_152.html
+
+K_11: 2.380e+2
+K_11,T: 2.068e-3
+K_11,T^2: -7.765e-5
+
+T_ref: 293.15
+
 C_p: 910.0
-K_11: 236.0

examples/config/phase/thermal/Au.yaml

@@ -1,4 +1,5 @@
 references:
   - https://de.wikipedia.org/wiki/Gold
+
 C_p: 128.0
 K_11: 320.0

examples/config/phase/thermal/Cu.yaml

@@ -1,4 +1,14 @@
 references:
+  - J.G. Hust and A.B. Lankford,
+    Thermal Conductivity of Aluminum, Copper, Iron, and Tungsten from 1K to the Melting Point,
+    US Department of Commerce, Boulder, Colorado, 1984,
+    fit to Tab. 2.4.1 (RRR=1000, T_min=200K, T_max=1000K)
   - https://www.mit.edu/~6.777/matprops/copper.htm
+
+K_11: 4.039e+2
+K_11,T: -8.119e-2
+K_11,T^2: 1.454e-5
+
+T_ref: 293.15
+
 C_p: 385.0
-K_11: 401.0

examples/config/phase/thermal/Fe.yaml

@@ -0,0 +1,14 @@
+references:
+  - J.G. Hust and A.B. Lankford,
+    Thermal Conductivity of Aluminum, Copper, Iron, and Tungsten from 1K to the Melting Point,
+    US Department of Commerce, Boulder, Colorado, 1984,
+    fit to Tab. 4.4.1 (RRR=300, T_min=200K, T_max=1000K)
+  - https://www.engineeringtoolbox.com/specific-heat-metals-d_152.html
+
+K_11: 8.055e+1
+K_11,T: -1.051e-1
+K_11,T^2: 5.464e-5
+
+T_ref: 293.15
+
+C_p: 450.0

examples/config/phase/thermal/Ni.yaml

@@ -0,0 +1,14 @@
+references:
+  - Y.S. Touloukian et al.,
+    TPRC Data Series Volume 1. Thermal conductivity - metallic elements and alloys,
+    IFI/Plenum, 1970,
+    fit to Tab. 35R (T_min=150K, T_max=500K)
+  - https://www.engineeringtoolbox.com/specific-heat-metals-d_152.html
+
+K_11: 9.132e+1
+K_11,T: -1.525e-1
+K_11,T^2: 3.053e-4
+
+T_ref: 293.15
+
+C_p: 440.0

examples/config/phase/thermal/Sn-beta.yaml

@@ -0,0 +1,18 @@
+references:
+  - Y.S. Touloukian et al.,
+    TPRC Data Series Volume 1. Thermal conductivity - metallic elements and alloys,
+    IFI/Plenum, 1970,
+    fit to Tab. 61R (T_min=100K, T_max=400K)
+  - https://www.engineeringtoolbox.com/specific-heat-metals-d_152.html
+
+K_11: 7.414e+1
+K_11,T: -6.465e-2
+K_11,T^2: 2.066e-4
+
+K_33: 5.147e+1
+K_33,T: -4.506e-2
+K_33,T^2: 1.435e-4
+
+T_ref: 293.15
+
+C_p: 210.0

examples/config/phase/thermal/W.yaml

@@ -1,4 +1,14 @@
 references:
+  - J.G. Hust and A.B. Lankford,
+    Thermal Conductivity of Aluminum, Copper, Iron, and Tungsten from 1K to the Melting Point,
+    US Department of Commerce, Boulder, Colorado, 1984,
+    fit to Tab. 5.4.1 (RRR=300, T_min=200K, T_max=1000K)
   - https://www.mit.edu/~6.777/matprops/tungsten.htm
+
+K_11: 1.758e+2
+K_11,T: -1.605e-1
+K_11,T^2: 1.160e-4
+
+T_ref: 293.15
+
 C_p: 132.51
-K_11: 178.0

examples/config/phase/thermal/source/dissipation_generic.yaml

@@ -1,2 +1,3 @@
 type: dissipation
+
 kappa: .9

examples/config/phase/thermal/source/externalheat_ramp-and-hold.yaml

@@ -1,3 +1,4 @@
 type: externalheat
+
 f_T: [1, 1, 0, 0]
 t_n: [0, 500, 500.001, 1000]

examples/config/phase/thermal/steel-0.5C.yaml

@@ -1,5 +1,6 @@
 references:
   - https://www.engineeringtoolbox.com/thermal-conductivity-metals-d_858.html
   - https://www.engineeringtoolbox.com/specific-heat-metals-d_152.html
+
 C_p: 490.0
 K_11: 54.0

python/damask/VERSION

@@ -1 +1 @@
-v3.0.0-alpha5-608-g3e8d1a60d
+v3.0.0-alpha5-696-g6fce27dee

python/damask/_configmaterial.py

@@ -178,8 +178,9 @@ class ConfigMaterial(Config):
         table : damask.Table
             Table that contains material information.
         **kwargs
-            Keyword arguments where the key is the name and the value specifies
+            Keyword arguments where the key is the property name and
-            the label of the data column in the table.
+            the value specifies either the label of the data column in the table
+            or a constant value.
 
         Returns
         -------
@@ -211,8 +212,23 @@ class ConfigMaterial(Config):
         homogenization: {}
         phase: {}
 
+        >>> cm.from_table(t,O='qu',phase='phase',homogenization='single_crystal')
+        material:
+          - constituents:
+              - O: [0.19, 0.8, 0.24, -0.51]
+                v: 1.0
+                phase: Aluminum
+            homogenization: single_crystal
+          - constituents:
+              - O: [0.8, 0.19, 0.24, -0.51]
+                v: 1.0
+                phase: Steel
+            homogenization: single_crystal
+        homogenization: {}
+        phase: {}
+
         """
-        kwargs_ = {k:table.get(v) for k,v in kwargs.items()}
+        kwargs_ = {k:table.get(v) if v in table.labels else np.atleast_2d([v]*len(table)).T for k,v in kwargs.items()}
 
         _,idx = np.unique(np.hstack(list(kwargs_.values())),return_index=True,axis=0)
         idx = np.sort(idx)

python/damask/_crystal.py

@@ -2,10 +2,11 @@ from typing import Union, Dict, List, Tuple
 
 import numpy as np
 
+from ._typehints import FloatSequence, CrystalFamily, CrystalLattice, CrystalKinematics
 from . import util
 from . import Rotation
 
-lattice_symmetries = {
+lattice_symmetries: Dict[CrystalLattice, CrystalFamily] = {
                 'aP': 'triclinic',
 
                 'mP': 'monoclinic',
@@ -30,9 +31,9 @@ lattice_symmetries = {
 class Crystal():
     """Crystal lattice."""
 
-    def __init__(self,*,
+    def __init__(self, *,
-                 family = None,
+                 family: CrystalFamily = None,
-                 lattice = None,
+                 lattice: CrystalLattice = None,
                  a: float = None, b: float = None, c: float = None,
                  alpha: float = None, beta: float = None, gamma: float = None,
                  degrees: bool = False):
@@ -130,9 +131,8 @@ class Crystal():
             Crystal to check for equality.
 
         """
-        if not isinstance(other, Crystal):
+        return NotImplemented if not isinstance(other, Crystal) else \
-            return NotImplemented
+               self.lattice == other.lattice and \
-        return self.lattice == other.lattice and \
                self.parameters == other.parameters and \
                self.family == other.family
 
@@ -208,7 +208,7 @@ class Crystal():
         ...    }
 
         """
-        _basis  = {
+        _basis: Dict[CrystalFamily, Dict[str, np.ndarray]]  = {
             'cubic':    {'improper':np.array([ [-1.            ,  0.            ,  1. ],
                                                [ np.sqrt(2.)   , -np.sqrt(2.)   ,  0. ],
                                                [ 0.            ,  np.sqrt(3.)   ,  0. ] ]),
@@ -315,19 +315,19 @@ class Crystal():
                                               self.lattice[-1],None),dtype=float)
 
     def to_lattice(self, *,
-                   direction: np.ndarray = None,
+                   direction: FloatSequence = None,
-                   plane: np.ndarray = None) -> np.ndarray:
+                   plane: FloatSequence = None) -> np.ndarray:
         """
         Calculate lattice vector corresponding to crystal frame direction or plane normal.
 
         Parameters
         ----------
-        direction|plane : numpy.ndarray of shape (...,3)
+        direction|plane : numpy.ndarray, shape (...,3)
             Vector along direction or plane normal.
 
         Returns
         -------
-        Miller : numpy.ndarray of shape (...,3)
+        Miller : numpy.ndarray, shape (...,3)
             Lattice vector of direction or plane.
             Use util.scale_to_coprime to convert to (integer) Miller indices.
 
@@ -341,19 +341,19 @@ class Crystal():
 
 
     def to_frame(self, *,
-                 uvw: np.ndarray = None,
+                 uvw: FloatSequence = None,
-                 hkl: np.ndarray = None) -> np.ndarray:
+                 hkl: FloatSequence = None) -> np.ndarray:
         """
         Calculate crystal frame vector along lattice direction [uvw] or plane normal (hkl).
 
         Parameters
         ----------
-        uvw|hkl : numpy.ndarray of shape (...,3)
+        uvw|hkl : numpy.ndarray, shape (...,3)
             Miller indices of crystallographic direction or plane normal.
 
         Returns
         -------
-        vector : numpy.ndarray of shape (...,3)
+        vector : numpy.ndarray, shape (...,3)
             Crystal frame vector along [uvw] direction or (hkl) plane normal.
 
         """
@@ -366,7 +366,7 @@ class Crystal():
 
 
     def kinematics(self,
-                   mode: str) -> Dict[str, List[np.ndarray]]:
+                   mode: CrystalKinematics) -> Dict[str, List[np.ndarray]]:
         """
         Return crystal kinematics systems.
 
@@ -381,7 +381,7 @@ class Crystal():
             Directions and planes of deformation mode families.
 
         """
-        _kinematics = {
+        _kinematics: Dict[CrystalLattice, Dict[CrystalKinematics, List[np.ndarray]]] = {
             'cF': {
                 'slip': [np.array([
                            [+0,+1,-1, +1,+1,+1],
@@ -626,7 +626,7 @@ class Crystal():
 
 
     def relation_operations(self,
-                            model: str) -> Tuple[str, Rotation]:
+                            model: str) -> Tuple[CrystalLattice, Rotation]:
         """
         Crystallographic orientation relationships for phase transformations.
 
@@ -658,7 +658,7 @@ class Crystal():
         https://doi.org/10.1016/j.actamat.2004.11.021
 
         """
-        _orientation_relationships = {
+        _orientation_relationships: Dict[str, Dict[CrystalLattice,np.ndarray]] = {
           'KS': {
             'cF' : np.array([
                 [[-1, 0, 1],[ 1, 1, 1]],

python/damask/_orientation.py

@@ -1,8 +1,10 @@
 import inspect
 import copy
+from typing import Union, Callable, List, Dict, Any, Tuple, TypeVar
 
 import numpy as np
 
+from ._typehints import FloatSequence, IntSequence, CrystalFamily, CrystalLattice
 from . import Rotation
 from . import Crystal
 from . import util
@@ -33,6 +35,7 @@ _parameter_doc = \
 
        """
 
+MyType = TypeVar('MyType', bound='Orientation')
 
 class Orientation(Rotation,Crystal):
     """
@@ -93,12 +96,13 @@ class Orientation(Rotation,Crystal):
 
     @util.extend_docstring(_parameter_doc)
     def __init__(self,
-                 rotation = np.array([1.0,0.0,0.0,0.0]), *,
+                 rotation: Union[FloatSequence, Rotation] = np.array([1.,0.,0.,0.]),
-                 family = None,
+                 *,
-                 lattice = None,
+                 family: CrystalFamily = None,
-                 a = None,b = None,c = None,
+                 lattice: CrystalLattice = None,
-                 alpha = None,beta = None,gamma = None,
+                 a: float = None, b: float = None, c: float = None,
-                 degrees = False):
+                 alpha: float = None, beta: float = None, gamma: float = None,
+                 degrees: bool = False):
         """
         New orientation.
 
@@ -115,13 +119,13 @@ class Orientation(Rotation,Crystal):
                               a=a,b=b,c=c, alpha=alpha,beta=beta,gamma=gamma, degrees=degrees)
 
 
-    def __repr__(self):
+    def __repr__(self) -> str:
         """Represent."""
         return '\n'.join([Crystal.__repr__(self),
                           Rotation.__repr__(self)])
 
-
+    def __copy__(self: MyType,
-    def __copy__(self,rotation=None):
+                 rotation: Union[FloatSequence, Rotation] = None) -> MyType:
         """Create deep copy."""
         dup = copy.deepcopy(self)
         if rotation is not None:
@@ -131,7 +135,9 @@ class Orientation(Rotation,Crystal):
     copy = __copy__
 
 
-    def __eq__(self,other):
+
+    def __eq__(self,
+               other: object) -> bool:
         """
         Equal to other.
 
@@ -141,12 +147,15 @@ class Orientation(Rotation,Crystal):
             Orientation to check for equality.
 
         """
+        if not isinstance(other, Orientation):
+            return NotImplemented
         matching_type = self.family == other.family and \
                         self.lattice == other.lattice and \
                         self.parameters == other.parameters
         return np.logical_and(matching_type,super(self.__class__,self.reduced).__eq__(other.reduced))
 
-    def __ne__(self,other):
+    def __ne__(self,
+               other: object) -> bool:
         """
         Not equal to other.
 
@@ -156,10 +165,14 @@ class Orientation(Rotation,Crystal):
             Orientation to check for equality.
 
         """
-        return np.logical_not(self==other)
+        return np.logical_not(self==other) if isinstance(other, Orientation) else NotImplemented
 
 
-    def isclose(self,other,rtol=1e-5,atol=1e-8,equal_nan=True):
+    def isclose(self: MyType,
+                other: MyType,
+                rtol: float = 1e-5,
+                atol: float = 1e-8,
+                equal_nan: bool = True) -> bool:
         """
         Report where values are approximately equal to corresponding ones of other Orientation.
 
@@ -176,7 +189,7 @@ class Orientation(Rotation,Crystal):
 
         Returns
         -------
-        mask : numpy.ndarray bool
+        mask : numpy.ndarray of bool, shape (self.shape)
             Mask indicating where corresponding orientations are close.
 
         """
@@ -187,7 +200,11 @@ class Orientation(Rotation,Crystal):
 
 
 
-    def allclose(self,other,rtol=1e-5,atol=1e-8,equal_nan=True):
+    def allclose(self: MyType,
+                 other: MyType,
+                 rtol: float = 1e-5,
+                 atol: float = 1e-8,
+                 equal_nan: bool = True) -> bool:
         """
         Test whether all values are approximately equal to corresponding ones of other Orientation.
 
@@ -208,10 +225,11 @@ class Orientation(Rotation,Crystal):
             Whether all values are close between both orientations.
 
         """
-        return np.all(self.isclose(other,rtol,atol,equal_nan))
+        return bool(np.all(self.isclose(other,rtol,atol,equal_nan)))
 
 
-    def __mul__(self,other):
+    def __mul__(self: MyType,
+                other: Union[Rotation, 'Orientation']) -> MyType:
         """
         Compose this orientation with other.
 
@@ -226,14 +244,15 @@ class Orientation(Rotation,Crystal):
             Compound rotation self*other, i.e. first other then self rotation.
 
         """
-        if isinstance(other,Orientation) or isinstance(other,Rotation):
+        if isinstance(other, (Orientation,Rotation)):
-            return self.copy(rotation=Rotation.__mul__(self,Rotation(other.quaternion)))
+            return self.copy(Rotation(self.quaternion)*Rotation(other.quaternion))
         else:
             raise TypeError('use "O@b", i.e. matmul, to apply Orientation "O" to object "b"')
 
 
     @staticmethod
-    def _split_kwargs(kwargs,target):
+    def _split_kwargs(kwargs: Dict[str, Any],
+                      target: Callable) -> Tuple[Dict[str, Any], ...]:
         """
         Separate keyword arguments in 'kwargs' targeted at 'target' from general keyword arguments of Orientation objects.
 
@@ -252,7 +271,7 @@ class Orientation(Rotation,Crystal):
             Valid keyword arguments of Orientation object.
 
         """
-        kws = ()
+        kws: Tuple[Dict[str, Any], ...] = ()
         for t in (target,Orientation.__init__):
             kws += ({key: kwargs[key] for key in set(inspect.signature(t).parameters) & set(kwargs)},)
 
@@ -264,105 +283,108 @@ class Orientation(Rotation,Crystal):
 
 
     @classmethod
-    @util.extended_docstring(Rotation.from_random,_parameter_doc)
+    @util.extended_docstring(Rotation.from_random, _parameter_doc)
-    def from_random(cls,**kwargs):
+    def from_random(cls, **kwargs) -> 'Orientation':
         kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_random)
         return cls(rotation=Rotation.from_random(**kwargs_rot),**kwargs_ori)
 
 
     @classmethod
     @util.extended_docstring(Rotation.from_quaternion,_parameter_doc)
-    def from_quaternion(cls,**kwargs):
+    def from_quaternion(cls, **kwargs) -> 'Orientation':
         kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_quaternion)
         return cls(rotation=Rotation.from_quaternion(**kwargs_rot),**kwargs_ori)
 
 
     @classmethod
     @util.extended_docstring(Rotation.from_Euler_angles,_parameter_doc)
-    def from_Euler_angles(cls,**kwargs):
+    def from_Euler_angles(cls, **kwargs) -> 'Orientation':
         kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_Euler_angles)
         return cls(rotation=Rotation.from_Euler_angles(**kwargs_rot),**kwargs_ori)
 
 
     @classmethod
     @util.extended_docstring(Rotation.from_axis_angle,_parameter_doc)
-    def from_axis_angle(cls,**kwargs):
+    def from_axis_angle(cls, **kwargs) -> 'Orientation':
         kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_axis_angle)
         return cls(rotation=Rotation.from_axis_angle(**kwargs_rot),**kwargs_ori)
 
 
     @classmethod
     @util.extended_docstring(Rotation.from_basis,_parameter_doc)
-    def from_basis(cls,**kwargs):
+    def from_basis(cls, **kwargs) -> 'Orientation':
         kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_basis)
         return cls(rotation=Rotation.from_basis(**kwargs_rot),**kwargs_ori)
 
 
     @classmethod
     @util.extended_docstring(Rotation.from_matrix,_parameter_doc)
-    def from_matrix(cls,**kwargs):
+    def from_matrix(cls, **kwargs) -> 'Orientation':
         kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_matrix)
         return cls(rotation=Rotation.from_matrix(**kwargs_rot),**kwargs_ori)
 
 
     @classmethod
     @util.extended_docstring(Rotation.from_Rodrigues_vector,_parameter_doc)
-    def from_Rodrigues_vector(cls,**kwargs):
+    def from_Rodrigues_vector(cls, **kwargs) -> 'Orientation':
         kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_Rodrigues_vector)
         return cls(rotation=Rotation.from_Rodrigues_vector(**kwargs_rot),**kwargs_ori)
 
 
     @classmethod
     @util.extended_docstring(Rotation.from_homochoric,_parameter_doc)
-    def from_homochoric(cls,**kwargs):
+    def from_homochoric(cls, **kwargs) -> 'Orientation':
         kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_homochoric)
         return cls(rotation=Rotation.from_homochoric(**kwargs_rot),**kwargs_ori)
 
 
     @classmethod
     @util.extended_docstring(Rotation.from_cubochoric,_parameter_doc)
-    def from_cubochoric(cls,**kwargs):
+    def from_cubochoric(cls, **kwargs) -> 'Orientation':
         kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_cubochoric)
         return cls(rotation=Rotation.from_cubochoric(**kwargs_rot),**kwargs_ori)
 
 
     @classmethod
     @util.extended_docstring(Rotation.from_spherical_component,_parameter_doc)
-    def from_spherical_component(cls,**kwargs):
+    def from_spherical_component(cls, **kwargs) -> 'Orientation':
         kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_spherical_component)
         return cls(rotation=Rotation.from_spherical_component(**kwargs_rot),**kwargs_ori)
 
 
     @classmethod
     @util.extended_docstring(Rotation.from_fiber_component,_parameter_doc)
-    def from_fiber_component(cls,**kwargs):
+    def from_fiber_component(cls, **kwargs) -> 'Orientation':
         kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_fiber_component)
         return cls(rotation=Rotation.from_fiber_component(**kwargs_rot),**kwargs_ori)
 
 
     @classmethod
     @util.extend_docstring(_parameter_doc)
-    def from_directions(cls,uvw,hkl,**kwargs):
+    def from_directions(cls,
+                        uvw: FloatSequence,
+                        hkl: FloatSequence,
+                        **kwargs) -> 'Orientation':
         """
         Initialize orientation object from two crystallographic directions.
 
         Parameters
         ----------
-        uvw : list, numpy.ndarray of shape (...,3)
+        uvw : numpy.ndarray, shape (...,3)
-            lattice direction aligned with lab frame x-direction.
+            Lattice direction aligned with lab frame x-direction.
-        hkl : list, numpy.ndarray of shape (...,3)
+        hkl : numpy.ndarray, shape (...,3)
-            lattice plane normal aligned with lab frame z-direction.
+            Lattice plane normal aligned with lab frame z-direction.
 
         """
         o = cls(**kwargs)
         x = o.to_frame(uvw=uvw)
         z = o.to_frame(hkl=hkl)
         om = np.stack([x,np.cross(z,x),z],axis=-2)
-        return o.copy(rotation=Rotation.from_matrix(tensor.transpose(om/np.linalg.norm(om,axis=-1,keepdims=True))))
+        return o.copy(Rotation.from_matrix(tensor.transpose(om/np.linalg.norm(om,axis=-1,keepdims=True))))
 
 
     @property
-    def equivalent(self):
+    def equivalent(self: MyType) -> MyType:
         """
         Orientations that are symmetrically equivalent.
 
@@ -372,11 +394,11 @@ class Orientation(Rotation,Crystal):
         """
         sym_ops = self.symmetry_operations
         o = sym_ops.broadcast_to(sym_ops.shape+self.shape,mode='right')
-        return self.copy(rotation=o*Rotation(self.quaternion).broadcast_to(o.shape,mode='left'))
+        return self.copy(o*Rotation(self.quaternion).broadcast_to(o.shape,mode='left'))
 
 
     @property
-    def reduced(self):
+    def reduced(self: MyType) -> MyType:
         """Select symmetrically equivalent orientation that falls into fundamental zone according to symmetry."""
         eq   = self.equivalent
         ok   = eq.in_FZ
@@ -387,13 +409,13 @@ class Orientation(Rotation,Crystal):
 
 
     @property
-    def in_FZ(self):
+    def in_FZ(self) -> Union[np.bool_, np.ndarray]:
         """
         Check whether orientation falls into fundamental zone of own symmetry.
 
         Returns
         -------
-        in : numpy.ndarray of bool, quaternion.shape
+        in : numpy.ndarray of bool, shape (self.shape)
             Whether Rodrigues-Frank vector falls into fundamental zone.
 
         Notes
@@ -431,13 +453,13 @@ class Orientation(Rotation,Crystal):
 
 
     @property
-    def in_disorientation_FZ(self):
+    def in_disorientation_FZ(self) -> np.ndarray:
         """
         Check whether orientation falls into fundamental zone of disorientations.
 
         Returns
         -------
-        in : numpy.ndarray of bool, quaternion.shape
+        in : numpy.ndarray of bool, shape (self.shape)
             Whether Rodrigues-Frank vector falls into disorientation FZ.
 
         References
@@ -471,8 +493,9 @@ class Orientation(Rotation,Crystal):
             else:
                 return np.ones_like(rho[...,0],dtype=bool)
 
-
+    def disorientation(self,
-    def disorientation(self,other,return_operators=False):
+                       other: 'Orientation',
+                       return_operators: bool = False) -> object:
         """
         Calculate disorientation between myself and given other orientation.
 
@@ -490,7 +513,7 @@ class Orientation(Rotation,Crystal):
         -------
         disorientation : Orientation
             Disorientation between self and other.
-        operators : numpy.ndarray int of shape (...,2), conditional
+        operators : numpy.ndarray of int, shape (...,2), conditional
             Index of symmetrically equivalent orientation that rotated vector to the SST.
 
         Notes
@@ -530,7 +553,7 @@ class Orientation(Rotation,Crystal):
             raise NotImplementedError('disorientation between different crystal families')
 
         blend = util.shapeblender(self.shape,other.shape)
-        s =  self.equivalent
+        s = self.equivalent
         o = other.equivalent
 
         s_ = s.reshape((s.shape[0],1)+ self.shape).broadcast_to((s.shape[0],o.shape[0])+blend,mode='right')
@@ -557,13 +580,15 @@ class Orientation(Rotation,Crystal):
                )
 
 
-    def average(self,weights=None,return_cloud=False):
+    def average(self,
+                weights: FloatSequence = None,
+                return_cloud: bool = False):
         """
         Return orientation average over last dimension.
 
         Parameters
         ----------
-        weights : numpy.ndarray, optional
+        weights : numpy.ndarray, shape (self.shape), optional
             Relative weights of orientations.
         return_cloud : bool, optional
             Return the set of symmetrically equivalent orientations that was used in averaging.
@@ -583,31 +608,30 @@ class Orientation(Rotation,Crystal):
 
         """
         eq = self.equivalent
-        m  = eq.misorientation(self[...,0].reshape((1,)+self.shape[:-1]+(1,))
+        m  = eq.misorientation(self[...,0].reshape((1,)+self.shape[:-1]+(1,))             # type: ignore
-                                          .broadcast_to(eq.shape))\
+                                          .broadcast_to(eq.shape)).as_axis_angle()[...,3] # type: ignore
-               .as_axis_angle()[...,3]
         r = Rotation(np.squeeze(np.take_along_axis(eq.quaternion,
                                                    np.argmin(m,axis=0)[np.newaxis,...,np.newaxis],
                                                    axis=0),
                                 axis=0))
-        return (
+        return ((self.copy(Rotation(r).average(weights)),self.copy(Rotation(r))) if return_cloud else
-                (self.copy(rotation=Rotation(r).average(weights)),
+                self.copy(Rotation(r).average(weights))
-                 self.copy(rotation=Rotation(r)))
-                if return_cloud else
-                self.copy(rotation=Rotation(r).average(weights))
                )
 
 
-    def to_SST(self,vector,proper=False,return_operators=False):
+    def to_SST(self,
+               vector: FloatSequence,
+               proper: bool = False,
+               return_operators: bool = False) -> np.ndarray:
         """
         Rotate vector to ensure it falls into (improper or proper) standard stereographic triangle of crystal symmetry.
 
         Parameters
         ----------
-        vector : numpy.ndarray of shape (...,3)
+        vector : numpy.ndarray, shape (...,3)
             Lab frame vector to align with crystal frame direction.
             Shape of vector blends with shape of own rotation array.
-            For example, a rotation array of shape (3,2) and a (2,4) vector array result in (3,2,4) outputs.
+            For example, a rotation array of shape (3,2) and a vector array of shape (2,4) result in (3,2,4) outputs.
         proper : bool, optional
             Consider only vectors with z >= 0, hence combine two neighboring SSTs.
             Defaults to False.
@@ -617,15 +641,18 @@ class Orientation(Rotation,Crystal):
 
         Returns
         -------
-        vector_SST : numpy.ndarray of shape (...,3)
+        vector_SST : numpy.ndarray, shape (...,3)
             Rotated vector falling into SST.
-        operators : numpy.ndarray int of shape (...), conditional
+        operators : numpy.ndarray of int, shape (...), conditional
             Index of symmetrically equivalent orientation that rotated vector to SST.
 
         """
+        vector_ = np.array(vector,float)
+        if vector_.shape[-1] != 3:
+            raise ValueError('input is not a field of three-dimensional vectors')
         eq  = self.equivalent
-        blend = util.shapeblender(eq.shape,np.array(vector).shape[:-1])
+        blend = util.shapeblender(eq.shape,vector_.shape[:-1])
-        poles = eq.broadcast_to(blend,mode='right') @ np.broadcast_to(np.array(vector),blend+(3,))
+        poles = eq.broadcast_to(blend,mode='right') @ np.broadcast_to(vector_,blend+(3,))
         ok    = self.in_SST(poles,proper=proper)
         ok   &= np.cumsum(ok,axis=0) == 1
         loc   = np.where(ok)
@@ -637,13 +664,15 @@ class Orientation(Rotation,Crystal):
                )
 
 
-    def in_SST(self,vector,proper=False):
+    def in_SST(self,
+               vector: FloatSequence,
+               proper: bool = False) -> Union[np.bool_, np.ndarray]:
         """
         Check whether given crystal frame vector falls into standard stereographic triangle of own symmetry.
 
         Parameters
         ----------
-        vector : numpy.ndarray of shape (...,3)
+        vector : numpy.ndarray, shape (...,3)
             Vector to check.
         proper : bool, optional
             Consider only vectors with z >= 0, hence combine two neighboring SSTs.
@@ -655,39 +684,43 @@ class Orientation(Rotation,Crystal):
             Whether vector falls into SST.
 
         """
-        if not isinstance(vector,np.ndarray) or vector.shape[-1] != 3:
+        vector_ = np.array(vector,float)
+        if vector_.shape[-1] != 3:
             raise ValueError('input is not a field of three-dimensional vectors')
 
         if self.standard_triangle is None:                                                          # direct exit for no symmetry
-            return  np.ones_like(vector[...,0],bool)
+            return np.ones_like(vector_[...,0],bool)
 
         if proper:
             components_proper   = np.around(np.einsum('...ji,...i',
-                                                      np.broadcast_to(self.standard_triangle['proper'], vector.shape+(3,)),
+                                                      np.broadcast_to(self.standard_triangle['proper'],   vector_.shape+(3,)),
-                                                      vector), 12)
+                                                      vector_), 12)
             components_improper = np.around(np.einsum('...ji,...i',
-                                                      np.broadcast_to(self.standard_triangle['improper'], vector.shape+(3,)),
+                                                      np.broadcast_to(self.standard_triangle['improper'], vector_.shape+(3,)),
-                                                      vector), 12)
+                                                      vector_), 12)
             return   np.all(components_proper   >= 0.0,axis=-1) \
                    | np.all(components_improper >= 0.0,axis=-1)
         else:
             components = np.around(np.einsum('...ji,...i',
-                                             np.broadcast_to(self.standard_triangle['improper'], vector.shape+(3,)),
+                                             np.broadcast_to(self.standard_triangle['improper'], vector_.shape+(3,)),
-                                             np.block([vector[...,:2],np.abs(vector[...,2:3])])), 12)
+                                             np.block([vector_[...,:2],np.abs(vector_[...,2:3])])), 12)
 
             return np.all(components >= 0.0,axis=-1)
 
 
-    def IPF_color(self,vector,in_SST=True,proper=False):
+    def IPF_color(self,
+                  vector: FloatSequence,
+                  in_SST: bool = True,
+                  proper: bool = False) -> np.ndarray:
         """
         Map vector to RGB color within standard stereographic triangle of own symmetry.
 
         Parameters
         ----------
-        vector : numpy.ndarray of shape (...,3)
+        vector : numpy.ndarray, shape (...,3)
             Vector to colorize.
             Shape of vector blends with shape of own rotation array.
-            For example, a rotation array of shape (3,2) and a (2,4) vector array result in (3,2,4) outputs.
+            For example, a rotation array of shape (3,2) and a vector array of shape (2,4) result in (3,2,4) outputs.
         in_SST : bool, optional
             Consider symmetrically equivalent orientations such that poles are located in SST.
             Defaults to True.
@@ -697,7 +730,7 @@ class Orientation(Rotation,Crystal):
 
         Returns
         -------
-        rgb : numpy.ndarray of shape (...,3)
+        rgb : numpy.ndarray, shape (...,3)
            RGB array of IPF colors.
 
         Examples
@@ -721,35 +754,35 @@ class Orientation(Rotation,Crystal):
 
         if proper:
             components_proper   = np.around(np.einsum('...ji,...i',
-                                                      np.broadcast_to(self.standard_triangle['proper'], vector_.shape+(3,)),
+                                                      np.broadcast_to(self.standard_triangle['proper'],   vector_.shape+(3,)),
                                                       vector_), 12)
             components_improper = np.around(np.einsum('...ji,...i',
                                                       np.broadcast_to(self.standard_triangle['improper'], vector_.shape+(3,)),
                                                       vector_), 12)
-            in_SST = np.all(components_proper   >= 0.0,axis=-1) \
+            in_SST_ = np.all(components_proper   >= 0.0,axis=-1) \
-                   | np.all(components_improper >= 0.0,axis=-1)
+                    | np.all(components_improper >= 0.0,axis=-1)
-            components = np.where((in_SST & np.all(components_proper   >= 0.0,axis=-1))[...,np.newaxis],
+            components = np.where((in_SST_ & np.all(components_proper   >= 0.0,axis=-1))[...,np.newaxis],
                                   components_proper,components_improper)
         else:
             components = np.around(np.einsum('...ji,...i',
                                              np.broadcast_to(self .standard_triangle['improper'], vector_.shape+(3,)),
                                              np.block([vector_[...,:2],np.abs(vector_[...,2:3])])), 12)
 
-            in_SST = np.all(components >= 0.0,axis=-1)
+            in_SST_ = np.all(components >= 0.0,axis=-1)
 
         with np.errstate(invalid='ignore',divide='ignore'):
             rgb = (components/np.linalg.norm(components,axis=-1,keepdims=True))**0.5                # smoothen color ramps
             rgb = np.clip(rgb,0.,1.)                                                                # clip intensity
             rgb /= np.max(rgb,axis=-1,keepdims=True)                                                # normalize to (HS)V = 1
-        rgb[np.broadcast_to(~in_SST[...,np.newaxis],rgb.shape)] = 0.0
+        rgb[np.broadcast_to(~in_SST_[...,np.newaxis],rgb.shape)] = 0.0
 
         return rgb
 
 
     @property
-    def symmetry_operations(self):
+    def symmetry_operations(self) -> Rotation:
         """Symmetry operations as Rotations."""
-        _symmetry_operations = {
+        _symmetry_operations: Dict[CrystalFamily, List]  = {
             'cubic':         [
                               [ 1.0,            0.0,            0.0,            0.0            ],
                               [ 0.0,            1.0,            0.0,            0.0            ],
@@ -819,22 +852,25 @@ class Orientation(Rotation,Crystal):
 ####################################################################################################
     # functions that require lattice, not just family
 
-    def to_pole(self,*,uvw=None,hkl=None,with_symmetry=False):
+    def to_pole(self, *,
+                uvw: FloatSequence = None,
+                hkl: FloatSequence = None,
+                with_symmetry: bool = False) -> np.ndarray:
         """
         Calculate lab frame vector along lattice direction [uvw] or plane normal (hkl).
 
         Parameters
         ----------
-        uvw|hkl : numpy.ndarray of shape (...,3)
+        uvw|hkl : numpy.ndarray, shape (...,3)
             Miller indices of crystallographic direction or plane normal.
             Shape of vector blends with shape of own rotation array.
-            For example, a rotation array of shape (3,2) and a (2,4) vector array result in (3,2,4) outputs.
+            For example, a rotation array, shape (3,2) and a vector array of shape (2,4) result in (3,2,4) outputs.
         with_symmetry : bool, optional
             Calculate all N symmetrically equivalent vectors.
 
         Returns
         -------
-        vector : numpy.ndarray of shape (...,3) or (...,N,3)
+        vector : numpy.ndarray, shape (...,3) or (...,N,3)
             Lab frame vector (or vectors if with_symmetry) along [uvw] direction or (hkl) plane normal.
 
         """
@@ -846,23 +882,24 @@ class Orientation(Rotation,Crystal):
             blend += sym_ops.shape
             v = sym_ops.broadcast_to(shape) \
               @ np.broadcast_to(v.reshape(util.shapeshifter(v.shape,shape+(3,))),shape+(3,))
-        return ~(self.broadcast_to(blend)) \
+        return ~(self.broadcast_to(blend))@ np.broadcast_to(v,blend+(3,))
-               @ np.broadcast_to(v,blend+(3,))
 
 
-    def Schmid(self,*,N_slip=None,N_twin=None):
+    def Schmid(self, *,
+               N_slip: IntSequence = None,
+               N_twin: IntSequence = None) -> np.ndarray:
         u"""
         Calculate Schmid matrix P = d ⨂ n in the lab frame for selected deformation systems.
 
         Parameters
         ----------
-        N_slip|N_twin : iterable of int
+        N_slip|N_twin : '*' or iterable of int
             Number of deformation systems per family of the deformation system.
             Use '*' to select all.
 
         Returns
         -------
-        P : numpy.ndarray of shape (N,...,3,3)
+        P : numpy.ndarray, shape (N,...,3,3)
             Schmid matrix for each of the N deformation systems.
 
         Examples
@@ -887,6 +924,8 @@ class Orientation(Rotation,Crystal):
                             (self.kinematics('twin'),N_twin)
         if active == '*': active = [len(a) for a in kinematics['direction']]
 
+        if not active:
+            raise ValueError('Schmid matrix not defined')
         d = self.to_frame(uvw=np.vstack([kinematics['direction'][i][:n] for i,n in enumerate(active)]))
         p = self.to_frame(hkl=np.vstack([kinematics['plane'][i][:n] for i,n in enumerate(active)]))
         P = np.einsum('...i,...j',d/np.linalg.norm(d,axis=1,keepdims=True),
@@ -897,7 +936,8 @@ class Orientation(Rotation,Crystal):
                @ np.broadcast_to(P.reshape(util.shapeshifter(P.shape,shape)),shape)
 
 
-    def related(self,model):
+    def related(self: MyType,
+                model: str) -> MyType:
         """
         Orientations derived from the given relationship.
 

python/damask/_table.py

@@ -44,6 +44,13 @@ class Table:
         return '\n'.join(['# '+c for c in self.comments])+'\n'+data_repr
 
 
+    def __eq__(self,
+               other: object) -> bool:
+        """Compare to other Table."""
+        return NotImplemented if not isinstance(other,Table) else \
+               self.shapes == other.shapes and self.data.equals(other.data)
+
+
     def __getitem__(self,
                     item: Union[slice, Tuple[slice, ...]]) -> 'Table':
         """
@@ -75,20 +82,22 @@ class Table:
            colB  colA
         0     1     0
         2     7     6
-        >>> tbl[1:2,'colB']
+        >>> tbl[[True,False,False,True],'colB']
            colB
-        1     4
+        0     1
-        2     7
+        3    10
 
         """
-        item = (item,slice(None,None,None)) if isinstance(item,slice) else \
+        item_ = (item,slice(None,None,None)) if isinstance(item,(slice,np.ndarray)) else \
-               item if isinstance(item[0],slice) else \
+                (np.array(item),slice(None,None,None)) if isinstance(item,list) and np.array(item).dtype == np.bool_ else \
-               (slice(None,None,None),item)
+                (np.array(item[0]),item[1]) if isinstance(item[0],list) else \
-        sliced = self.data.loc[item]
+                item if isinstance(item[0],(slice,np.ndarray)) else \
-        cols = np.array(sliced.columns if isinstance(sliced,pd.core.frame.DataFrame) else [item[1]])
+                (slice(None,None,None),item)
+        sliced = self.data.loc[item_]
+        cols = np.array(sliced.columns if isinstance(sliced,pd.core.frame.DataFrame) else [item_[1]])
         _,idx = np.unique(cols,return_index=True)
         return self.__class__(data=sliced,
-                              shapes = {k:self.shapes[k] for k in cols[np.sort(idx)]},
+                              shapes={k:self.shapes[k] for k in cols[np.sort(idx)]},
                               comments=self.comments)
 
 
@@ -185,7 +194,7 @@ class Table:
 
         Returns
         -------
-        mask : numpy.ndarray bool
+        mask : numpy.ndarray of bool
             Mask indicating where corresponding table values are close.
 
         """
@@ -363,9 +372,9 @@ class Table:
         label : str
             Column label.
         data : numpy.ndarray
-            New data.
+            Replacement data.
         info : str, optional
-            Human-readable information about the new data.
+            Human-readable information about the modified data.
 
         Returns
         -------
@@ -398,9 +407,9 @@ class Table:
         label : str
             Column label.
         data : numpy.ndarray
-            Modified data.
+            New data.
         info : str, optional
-            Human-readable information about the modified data.
+            Human-readable information about the new data.
 
         Returns
         -------
@@ -476,7 +485,7 @@ class Table:
                 labels: Union[str, List[str]],
                 ascending: Union[bool, List[bool]] = True) -> 'Table':
         """
-        Sort table by values of given labels.
+        Sort table by data of given columns.
 
         Parameters
         ----------

python/damask/_typehints.py

@@ -1,6 +1,6 @@
 """Functionality for typehints."""
 
-from typing import Sequence, Union, TextIO
+from typing import Sequence, Union, Literal, TextIO
 from pathlib import Path
 
 import numpy as np
@@ -9,6 +9,9 @@ import numpy as np
 FloatSequence = Union[np.ndarray,Sequence[float]]
 IntSequence = Union[np.ndarray,Sequence[int]]
 FileHandle = Union[TextIO, str, Path]
+CrystalFamily = Union[None,Literal['triclinic', 'monoclinic', 'orthorhombic', 'tetragonal', 'hexagonal', 'cubic']]
+CrystalLattice = Union[None,Literal['aP', 'mP', 'mS', 'oP', 'oS', 'oI', 'oF', 'tP', 'tI', 'hP', 'cP', 'cI', 'cF']]
+CrystalKinematics = Literal['slip', 'twin']
 NumpyRngSeed = Union[int, IntSequence, np.random.SeedSequence, np.random.Generator]
 # BitGenerator does not exists in older numpy versions
 #NumpyRngSeed = Union[int, IntSequence, np.random.SeedSequence, np.random.BitGenerator, np.random.Generator]

python/damask/util.py

@@ -9,7 +9,7 @@ import re
 import fractions
 from collections import abc
 from functools import reduce
-from typing import Union, Tuple, Iterable, Callable, Dict, List, Any, Literal, SupportsIndex, Sequence
+from typing import Union, Tuple, Iterable, Callable, Dict, List, Any, Literal
 from pathlib import Path
 
 import numpy as np
@@ -427,7 +427,7 @@ def hybrid_IA(dist: np.ndarray,
 def shapeshifter(fro: Tuple[int, ...],
                  to: Tuple[int, ...],
                  mode: Literal['left','right'] = 'left',
-                 keep_ones: bool = False) -> Sequence[SupportsIndex]:
+                 keep_ones: bool = False) -> Tuple[int, ...]:
     """
     Return dimensions that reshape 'fro' to become broadcastable to 'to'.
 
@@ -490,7 +490,7 @@ def shapeshifter(fro: Tuple[int, ...],
 
 
 def shapeblender(a: Tuple[int, ...],
-                 b: Tuple[int, ...]) -> Sequence[SupportsIndex]:
+                 b: Tuple[int, ...]) -> Tuple[int, ...]:
     """
     Return a shape that overlaps the rightmost entries of 'a' with the leftmost of 'b'.
 

python/tests/test_ConfigMaterial.py

@@ -96,6 +96,18 @@ class TestConfigMaterial:
         for i,m in enumerate(c['material']):
             assert m['homogenization'] == 1 and (m['constituents'][0]['O'] == [1,0,1,1]).all()
 
+    def test_from_table_with_constant(self):
+        N = np.random.randint(3,10)
+        a = np.vstack((np.hstack((np.arange(N),np.arange(N)[::-1])),
+                       np.ones(N*2),np.zeros(N*2),np.ones(N*2),np.ones(N*2),
+                       np.ones(N*2),
+                      )).T
+        t = Table(a,{'varying':1,'constant':4,'ones':1})
+        c = ConfigMaterial.from_table(t,**{'phase':'varying','O':'constant','homogenization':1})
+        assert len(c['material']) == N
+        for i,m in enumerate(c['material']):
+            assert m['homogenization'] == 1 and (m['constituents'][0]['O'] == [1,0,1,1]).all()
+
     @pytest.mark.parametrize('N,n,kw',[
                                         (1,1,{'phase':'Gold',
                                               'O':[1,0,0,0],

python/tests/test_Table.py

@@ -59,10 +59,14 @@ class TestTable:
 
     @pytest.mark.parametrize('N',[1,3,4])
     def test_slice(self,default,N):
+        mask = np.random.choice([True,False],len(default))
         assert len(default[:N]) == 1+N
         assert len(default[:N,['F','s']]) == 1+N
+        assert len(default[mask,['F','s']]) == np.count_nonzero(mask)
+        assert default[mask,['F','s']] == default[mask][['F','s']] == default[['F','s']][mask]
+        assert default[np.logical_not(mask),['F','s']] != default[mask][['F','s']]
         assert default[N:].get('F').shape == (len(default)-N,3,3)
-        assert (default[:N,['v','s']].data == default['v','s'][:N].data).all().all()
+        assert default[:N,['v','s']].data.equals(default['v','s'][:N].data)
 
     @pytest.mark.parametrize('mode',['str','path'])
     def test_write_read(self,default,tmp_path,mode):

src/CPFEM.f90

@@ -192,11 +192,10 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS
 
     else validCalculation
       if (debugCPFEM%extensive)  print'(a,i8,1x,i2)', '<< CPFEM >> calculation for elFE ip ',elFE,ip
-      call homogenization_mechanical_response(dt,[ip,ip],[elCP,elCP])
+      call homogenization_mechanical_response(dt,(elCP-1)*discretization_nIPs + ip,(elCP-1)*discretization_nIPs + ip)
       if (.not. terminallyIll) &
         call homogenization_mechanical_response2(dt,[ip,ip],[elCP,elCP])
 
-
       terminalIllness: if (terminallyIll) then
 
         call random_number(rnd)

src/discretization.f90

@@ -12,7 +12,8 @@ module discretization
 
   integer,     public, protected :: &
     discretization_nIPs, &
-    discretization_Nelems
+    discretization_Nelems, &
+    discretization_Ncells
 
   integer,     public, protected, dimension(:),   allocatable :: &
     discretization_materialAt                                                                       !ToDo: discretization_ID_material
@@ -53,6 +54,7 @@ subroutine discretization_init(materialAt,&
 
   discretization_Nelems = size(materialAt,1)
   discretization_nIPs   = size(IPcoords0,2)/discretization_Nelems
+  discretization_Ncells = discretization_Nelems*discretization_nIPs
 
   discretization_materialAt = materialAt
 

src/grid/spectral_utilities.f90

@@ -31,7 +31,7 @@ module spectral_utilities
 !--------------------------------------------------------------------------------------------------
 ! grid related information
   real(pReal), protected,  public                :: wgt                                             !< weighting factor 1/Nelems
-  integer,     protected,  public                :: grid1Red                                        !< cells(1)/2
+  integer,     protected,  public                :: cells1Red                                       !< cells(1)/2
   real(pReal), protected,  public,  dimension(3) :: scaledGeomSize                                  !< scaled geometry size for calculation of divergence
 
 !--------------------------------------------------------------------------------------------------
@@ -201,7 +201,7 @@ subroutine spectral_utilities_init
                                 num_grid%get_asString('PETSc_options',defaultVal=''),err_PETSc)
   CHKERRQ(err_PETSc)
 
-  grid1Red = cells(1)/2 + 1
+  cells1Red = cells(1)/2 + 1
   wgt = 1.0/real(product(cells),pReal)
 
   num%memory_efficient      = num_grid%get_asInt('memory_efficient',      defaultVal=1) > 0         ! ToDo: should be logical in YAML file
@@ -265,8 +265,8 @@ subroutine spectral_utilities_init
   gridFFTW = int(cells,C_INTPTR_T)
   alloc_local = fftw_mpi_local_size_3d(gridFFTW(3), gridFFTW(2), gridFFTW(1)/2 +1, &
                                        PETSC_COMM_WORLD, local_K, local_K_offset)
-  allocate (xi1st (3,grid1Red,cells(2),cells3),source = cmplx(0.0_pReal,0.0_pReal,pReal))           ! frequencies for first derivatives, only half the size for first dimension
+  allocate (xi1st (3,cells1Red,cells(2),cells3),source = cmplx(0.0_pReal,0.0_pReal,pReal))          ! frequencies for first derivatives, only half the size for first dimension
-  allocate (xi2nd (3,grid1Red,cells(2),cells3),source = cmplx(0.0_pReal,0.0_pReal,pReal))           ! frequencies for second derivatives, only half the size for first dimension
+  allocate (xi2nd (3,cells1Red,cells(2),cells3),source = cmplx(0.0_pReal,0.0_pReal,pReal))          ! frequencies for second derivatives, only half the size for first dimension
 
   tensorField = fftw_alloc_complex(tensorSize*alloc_local)
   call c_f_pointer(tensorField, tensorField_real,    [3_C_INTPTR_T,3_C_INTPTR_T, &
@@ -333,7 +333,7 @@ subroutine spectral_utilities_init
       do j = 1, cells(2)
         k_s(2) = j - 1
         if (j > cells(2)/2 + 1) k_s(2) = k_s(2) - cells(2)                                          ! running from 0,1,...,N/2,N/2+1,-N/2,-N/2+1,...,-1
-          do i = 1, grid1Red
+          do i = 1, cells1Red
             k_s(1) = i - 1                                                                          ! symmetry, junst running from 0,1,...,N/2,N/2+1
             xi2nd(1:3,i,j,k-cells3Offset) = utilities_getFreqDerivative(k_s)
             where(mod(cells,2)==0 .and. [i,j,k] == cells/2+1 .and. &
@@ -347,7 +347,7 @@ subroutine spectral_utilities_init
   if (num%memory_efficient) then                                                                    ! allocate just single fourth order tensor
     allocate (gamma_hat(3,3,3,3,1,1,1), source = cmplx(0.0_pReal,0.0_pReal,pReal))
   else                                                                                              ! precalculation of gamma_hat field
-    allocate (gamma_hat(3,3,3,3,grid1Red,cells(2),cells3), source = cmplx(0.0_pReal,0.0_pReal,pReal))
+    allocate (gamma_hat(3,3,3,3,cells1Red,cells(2),cells3), source = cmplx(0.0_pReal,0.0_pReal,pReal))
   endif
 
 end subroutine spectral_utilities_init
@@ -362,7 +362,7 @@ end subroutine spectral_utilities_init
 subroutine utilities_updateGamma(C)
 
   real(pReal), intent(in), dimension(3,3,3,3) :: C                                                  !< input stiffness to store as reference stiffness
-  complex(pReal),              dimension(3,3) :: temp33_complex, xiDyad_cmplx
+  complex(pReal),              dimension(3,3) :: temp33_cmplx, xiDyad_cmplx
   real(pReal),                 dimension(6,6) :: A, A_inv
   integer :: &
     i, j, k, &
@@ -372,27 +372,40 @@ subroutine utilities_updateGamma(C)
   C_ref = C
 
   if (.not. num%memory_efficient) then
-    gamma_hat =  cmplx(0.0_pReal,0.0_pReal,pReal)                                                   ! for the singular point and any non invertible A
+    gamma_hat = cmplx(0.0_pReal,0.0_pReal,pReal)                                                    ! for the singular point and any non invertible A
-    do k = cells3Offset+1, cells3Offset+cells3; do j = 1, cells(2); do i = 1, grid1Red
+    !$OMP PARALLEL DO PRIVATE(l,m,n,o,temp33_cmplx,xiDyad_cmplx,A,A_inv,err)
+    do k = cells3Offset+1, cells3Offset+cells3; do j = 1, cells(2); do i = 1, cells1Red
       if (any([i,j,k] /= 1)) then                                                                   ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1
-        do concurrent (l = 1:3, m = 1:3)
+#ifndef __INTEL_COMPILER
+        do concurrent(l = 1:3, m = 1:3)
           xiDyad_cmplx(l,m) = conjg(-xi1st(l,i,j,k-cells3Offset))*xi1st(m,i,j,k-cells3Offset)
         end do
         do concurrent(l = 1:3, m = 1:3)
-          temp33_complex(l,m) = sum(cmplx(C_ref(l,1:3,m,1:3),0.0_pReal)*xiDyad_cmplx)
+          temp33_cmplx(l,m) = sum(cmplx(C_ref(l,1:3,m,1:3),0.0_pReal)*xiDyad_cmplx)
         end do
-        A(1:3,1:3) = temp33_complex%re; A(4:6,4:6) =  temp33_complex%re
+#else
-        A(1:3,4:6) = temp33_complex%im; A(4:6,1:3) = -temp33_complex%im
+        forall(l = 1:3, m = 1:3) &
+          xiDyad_cmplx(l,m) = conjg(-xi1st(l,i,j,k-cells3Offset))*xi1st(m,i,j,k-cells3Offset)
+        forall(l = 1:3, m = 1:3) &
+          temp33_cmplx(l,m) = sum(cmplx(C_ref(l,1:3,m,1:3),0.0_pReal)*xiDyad_cmplx)
+#endif
+        A(1:3,1:3) = temp33_cmplx%re; A(4:6,4:6) =  temp33_cmplx%re
+        A(1:3,4:6) = temp33_cmplx%im; A(4:6,1:3) = -temp33_cmplx%im
         if (abs(math_det33(A(1:3,1:3))) > 1e-16) then
           call math_invert(A_inv, err, A)
-          temp33_complex = cmplx(A_inv(1:3,1:3),A_inv(1:3,4:6),pReal)
+          temp33_cmplx = cmplx(A_inv(1:3,1:3),A_inv(1:3,4:6),pReal)
+#ifndef __INTEL_COMPILER
           do concurrent(l=1:3, m=1:3, n=1:3, o=1:3)
-            gamma_hat(l,m,n,o,i,j,k-cells3Offset) = temp33_complex(l,n)* &
+            gamma_hat(l,m,n,o,i,j,k-cells3Offset) = temp33_cmplx(l,n) * xiDyad_cmplx(o,m)
-                                                    conjg(-xi1st(o,i,j,k-cells3Offset))*xi1st(m,i,j,k-cells3Offset)
           end do
+#else
+          forall(l=1:3, m=1:3, n=1:3, o=1:3) &
+            gamma_hat(l,m,n,o,i,j,k-cells3Offset) = temp33_cmplx(l,n) * xiDyad_cmplx(o,m)
+#endif
         end if
       end if
     end do; end do; end do
+    !$OMP END PARALLEL DO
   endif
 
 end subroutine utilities_updateGamma
@@ -405,7 +418,7 @@ end subroutine utilities_updateGamma
 !--------------------------------------------------------------------------------------------------
 subroutine utilities_FFTtensorForward
 
-  tensorField_real(1:3,1:3,cells(1)+1:grid1Red*2,:,:) = 0.0_pReal
+  tensorField_real(1:3,1:3,cells(1)+1:cells1Red*2,:,:) = 0.0_pReal
   call fftw_mpi_execute_dft_r2c(planTensorForth,tensorField_real,tensorField_fourier)
 
 end subroutine utilities_FFTtensorForward
@@ -429,7 +442,7 @@ end subroutine utilities_FFTtensorBackward
 !--------------------------------------------------------------------------------------------------
 subroutine utilities_FFTscalarForward
 
-  scalarField_real(cells(1)+1:grid1Red*2,:,:) = 0.0_pReal
+  scalarField_real(cells(1)+1:cells1Red*2,:,:) = 0.0_pReal
   call fftw_mpi_execute_dft_r2c(planScalarForth,scalarField_real,scalarField_fourier)
 
 end subroutine utilities_FFTscalarForward
@@ -454,7 +467,7 @@ end subroutine utilities_FFTscalarBackward
 !--------------------------------------------------------------------------------------------------
 subroutine utilities_FFTvectorForward
 
-  vectorField_real(1:3,cells(1)+1:grid1Red*2,:,:) = 0.0_pReal
+  vectorField_real(1:3,cells(1)+1:cells1Red*2,:,:) = 0.0_pReal
   call fftw_mpi_execute_dft_r2c(planVectorForth,vectorField_real,vectorField_fourier)
 
 end subroutine utilities_FFTvectorForward
@@ -478,7 +491,7 @@ end subroutine utilities_FFTvectorBackward
 subroutine utilities_fourierGammaConvolution(fieldAim)
 
   real(pReal), intent(in), dimension(3,3) :: fieldAim                                               !< desired average value of the field after convolution
-  complex(pReal),          dimension(3,3) :: temp33_complex, xiDyad_cmplx
+  complex(pReal),          dimension(3,3) :: temp33_cmplx, xiDyad_cmplx
   real(pReal),             dimension(6,6) :: A, A_inv
 
   integer :: &
@@ -493,38 +506,61 @@ subroutine utilities_fourierGammaConvolution(fieldAim)
 !--------------------------------------------------------------------------------------------------
 ! do the actual spectral method calculation (mechanical equilibrium)
   memoryEfficient: if (num%memory_efficient) then
-    do k = 1, cells3; do j = 1, cells(2); do i = 1, grid1Red
+    !$OMP PARALLEL DO PRIVATE(l,m,n,o,temp33_cmplx,xiDyad_cmplx,A,A_inv,err,gamma_hat)
+    do k = 1, cells3; do j = 1, cells(2); do i = 1, cells1Red
       if (any([i,j,k+cells3Offset] /= 1)) then                                                      ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1
+#ifndef __INTEL_COMPILER
         do concurrent(l = 1:3, m = 1:3)
           xiDyad_cmplx(l,m) = conjg(-xi1st(l,i,j,k))*xi1st(m,i,j,k)
         end do
         do concurrent(l = 1:3, m = 1:3)
-          temp33_complex(l,m) = sum(cmplx(C_ref(l,1:3,m,1:3),0.0_pReal)*xiDyad_cmplx)
+          temp33_cmplx(l,m) = sum(cmplx(C_ref(l,1:3,m,1:3),0.0_pReal)*xiDyad_cmplx)
         end do
-        A(1:3,1:3) = temp33_complex%re; A(4:6,4:6) =  temp33_complex%re
+#else
-        A(1:3,4:6) = temp33_complex%im; A(4:6,1:3) = -temp33_complex%im
+        forall(l = 1:3, m = 1:3) &
+          xiDyad_cmplx(l,m) = conjg(-xi1st(l,i,j,k))*xi1st(m,i,j,k)
+        forall(l = 1:3, m = 1:3) &
+          temp33_cmplx(l,m) = sum(cmplx(C_ref(l,1:3,m,1:3),0.0_pReal)*xiDyad_cmplx)
+#endif
+        A(1:3,1:3) = temp33_cmplx%re; A(4:6,4:6) =  temp33_cmplx%re
+        A(1:3,4:6) = temp33_cmplx%im; A(4:6,1:3) = -temp33_cmplx%im
         if (abs(math_det33(A(1:3,1:3))) > 1e-16) then
           call math_invert(A_inv, err, A)
-          temp33_complex = cmplx(A_inv(1:3,1:3),A_inv(1:3,4:6),pReal)
+          temp33_cmplx = cmplx(A_inv(1:3,1:3),A_inv(1:3,4:6),pReal)
+#ifndef __INTEL_COMPILER
           do concurrent(l=1:3, m=1:3, n=1:3, o=1:3)
-            gamma_hat(l,m,n,o,1,1,1) =  temp33_complex(l,n)*conjg(-xi1st(o,i,j,k))*xi1st(m,i,j,k)
+            gamma_hat(l,m,n,o,1,1,1) = temp33_cmplx(l,n)*xiDyad_cmplx(o,m)
           end do
+          do concurrent(l = 1:3, m = 1:3)
+            temp33_cmplx(l,m) = sum(gamma_hat(l,m,1:3,1:3,1,1,1)*tensorField_fourier(1:3,1:3,i,j,k))
+          end do
+#else
+          forall(l=1:3, m=1:3, n=1:3, o=1:3) &
+            gamma_hat(l,m,n,o,1,1,1) = temp33_cmplx(l,n)*xiDyad_cmplx(o,m)
+          forall(l = 1:3, m = 1:3) &
+            temp33_cmplx(l,m) = sum(gamma_hat(l,m,1:3,1:3,1,1,1)*tensorField_fourier(1:3,1:3,i,j,k))
+#endif
+          tensorField_fourier(1:3,1:3,i,j,k) = temp33_cmplx
         else
-          gamma_hat(1:3,1:3,1:3,1:3,1,1,1) = cmplx(0.0_pReal,0.0_pReal,pReal)
+          tensorField_fourier(1:3,1:3,i,j,k) = cmplx(0.0_pReal,0.0_pReal,pReal)
         end if
-        do concurrent(l = 1:3, m = 1:3)
-          temp33_Complex(l,m) = sum(gamma_hat(l,m,1:3,1:3,1,1,1)*tensorField_fourier(1:3,1:3,i,j,k))
-        end do
-        tensorField_fourier(1:3,1:3,i,j,k) = temp33_Complex
       end if
     end do; end do; end do
+    !$OMP END PARALLEL DO
   else memoryEfficient
-    do k = 1, cells3;  do j = 1, cells(2);  do i = 1,grid1Red
+    !$OMP PARALLEL DO PRIVATE(l,m,temp33_cmplx)
+    do k = 1, cells3;  do j = 1, cells(2);  do i = 1,cells1Red
+#ifndef __INTEL_COMPILER
       do concurrent(l = 1:3, m = 1:3)
-        temp33_Complex(l,m) = sum(gamma_hat(l,m,1:3,1:3,i,j,k) * tensorField_fourier(1:3,1:3,i,j,k))
+        temp33_cmplx(l,m) = sum(gamma_hat(l,m,1:3,1:3,i,j,k)*tensorField_fourier(1:3,1:3,i,j,k))
       end do
-      tensorField_fourier(1:3,1:3,i,j,k) = temp33_Complex
+#else
+      forall(l = 1:3, m = 1:3) &
+        temp33_cmplx(l,m) = sum(gamma_hat(l,m,1:3,1:3,i,j,k)*tensorField_fourier(1:3,1:3,i,j,k))
+#endif
+      tensorField_fourier(1:3,1:3,i,j,k) = temp33_cmplx
     end do; end do; end do
+    !$OMP END PARALLEL DO
   end if memoryEfficient
 
   if (cells3Offset == 0) tensorField_fourier(1:3,1:3,1,1,1) = cmplx(fieldAim/wgt,0.0_pReal,pReal)
@@ -544,12 +580,14 @@ subroutine utilities_fourierGreenConvolution(D_ref, mu_ref, Delta_t)
 
 !--------------------------------------------------------------------------------------------------
 ! do the actual spectral method calculation
-  do k = 1, cells3; do j = 1, cells(2) ;do i = 1, grid1Red
+  !$OMP PARALLEL DO PRIVATE(GreenOp_hat)
+  do k = 1, cells3; do j = 1, cells(2) ;do i = 1, cells1Red
     GreenOp_hat = cmplx(1.0_pReal,0.0_pReal,pReal) &
                 / (cmplx(mu_ref,0.0_pReal,pReal) + cmplx(Delta_t,0.0_pReal) &
                    * sum(conjg(xi1st(1:3,i,j,k))* matmul(cmplx(D_ref,0.0_pReal),xi1st(1:3,i,j,k))))
     scalarField_fourier(i,j,k) = scalarField_fourier(i,j,k)*GreenOp_hat
   enddo; enddo; enddo
+  !$OMP END PARALLEL DO
 
 end subroutine utilities_fourierGreenConvolution
 
@@ -572,7 +610,7 @@ real(pReal) function utilities_divergenceRMS()
 ! calculating RMS divergence criterion in Fourier space
   utilities_divergenceRMS = 0.0_pReal
   do k = 1, cells3; do j = 1, cells(2)
-    do i = 2, grid1Red -1                                                                           ! Has somewhere a conj. complex counterpart. Therefore count it twice.
+    do i = 2, cells1Red -1                                                                          ! Has somewhere a conj. complex counterpart. Therefore count it twice.
       utilities_divergenceRMS = utilities_divergenceRMS &
             + 2.0_pReal*(sum (real(matmul(tensorField_fourier(1:3,1:3,i,j,k), &                     ! (sqrt(real(a)**2 + aimag(a)**2))**2 = real(a)**2 + aimag(a)**2, i.e. do not take square root and square again
                                           conjg(-xi1st(1:3,i,j,k))*rescaledGeom))**2) &             ! --> sum squared L_2 norm of vector
@@ -584,10 +622,10 @@ real(pReal) function utilities_divergenceRMS()
                                   conjg(-xi1st(1:3,1,j,k))*rescaledGeom))**2) &
                + sum(aimag(matmul(tensorField_fourier(1:3,1:3,1       ,j,k), &
                                   conjg(-xi1st(1:3,1,j,k))*rescaledGeom))**2) &
-               + sum( real(matmul(tensorField_fourier(1:3,1:3,grid1Red,j,k), &
+               + sum( real(matmul(tensorField_fourier(1:3,1:3,cells1Red,j,k), &
-                                  conjg(-xi1st(1:3,grid1Red,j,k))*rescaledGeom))**2) &
+                                  conjg(-xi1st(1:3,cells1Red,j,k))*rescaledGeom))**2) &
-               + sum(aimag(matmul(tensorField_fourier(1:3,1:3,grid1Red,j,k), &
+               + sum(aimag(matmul(tensorField_fourier(1:3,1:3,cells1Red,j,k), &
-                                  conjg(-xi1st(1:3,grid1Red,j,k))*rescaledGeom))**2)
+                                  conjg(-xi1st(1:3,cells1Red,j,k))*rescaledGeom))**2)
   enddo; enddo
   if (cells(1) == 1) utilities_divergenceRMS = utilities_divergenceRMS * 0.5_pReal                   ! counted twice in case of cells(1) == 1
   call MPI_Allreduce(MPI_IN_PLACE,utilities_divergenceRMS,1_MPI_INTEGER_KIND,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,err_MPI)
@@ -617,7 +655,7 @@ real(pReal) function utilities_curlRMS()
   utilities_curlRMS = 0.0_pReal
 
   do k = 1, cells3; do j = 1, cells(2);
-    do i = 2, grid1Red - 1
+    do i = 2, cells1Red - 1
       do l = 1, 3
         curl_fourier(l,1) = (+tensorField_fourier(l,3,i,j,k)*xi1st(2,i,j,k)*rescaledGeom(2) &
                              -tensorField_fourier(l,2,i,j,k)*xi1st(3,i,j,k)*rescaledGeom(3))
@@ -640,12 +678,12 @@ real(pReal) function utilities_curlRMS()
     utilities_curlRMS = utilities_curlRMS &
                       + sum(curl_fourier%re**2 + curl_fourier%im**2)                                ! this layer (DC) does not have a conjugate complex counterpart (if cells(1) /= 1)
     do l = 1, 3
-      curl_fourier = (+tensorField_fourier(l,3,grid1Red,j,k)*xi1st(2,grid1Red,j,k)*rescaledGeom(2) &
+      curl_fourier = (+tensorField_fourier(l,3,cells1Red,j,k)*xi1st(2,cells1Red,j,k)*rescaledGeom(2) &
-                      -tensorField_fourier(l,2,grid1Red,j,k)*xi1st(3,grid1Red,j,k)*rescaledGeom(3))
+                      -tensorField_fourier(l,2,cells1Red,j,k)*xi1st(3,cells1Red,j,k)*rescaledGeom(3))
-      curl_fourier = (+tensorField_fourier(l,1,grid1Red,j,k)*xi1st(3,grid1Red,j,k)*rescaledGeom(3) &
+      curl_fourier = (+tensorField_fourier(l,1,cells1Red,j,k)*xi1st(3,cells1Red,j,k)*rescaledGeom(3) &
-                      -tensorField_fourier(l,3,grid1Red,j,k)*xi1st(1,grid1Red,j,k)*rescaledGeom(1))
+                      -tensorField_fourier(l,3,cells1Red,j,k)*xi1st(1,cells1Red,j,k)*rescaledGeom(1))
-      curl_fourier = (+tensorField_fourier(l,2,grid1Red,j,k)*xi1st(1,grid1Red,j,k)*rescaledGeom(1) &
+      curl_fourier = (+tensorField_fourier(l,2,cells1Red,j,k)*xi1st(1,cells1Red,j,k)*rescaledGeom(1) &
-                      -tensorField_fourier(l,1,grid1Red,j,k)*xi1st(2,grid1Red,j,k)*rescaledGeom(2))
+                      -tensorField_fourier(l,1,cells1Red,j,k)*xi1st(2,cells1Red,j,k)*rescaledGeom(2))
     enddo
     utilities_curlRMS = utilities_curlRMS &
                       + sum(curl_fourier%re**2 + curl_fourier%im**2)                                ! this layer (Nyquist) does not have a conjugate complex counterpart (if cells(1) /= 1)
@@ -736,9 +774,10 @@ subroutine utilities_fourierScalarGradient()
 
   integer :: i, j, k
 
-  do k = 1, cells3;  do j = 1, cells(2);  do i = 1,grid1Red
+
+  do k = 1, cells3;  do j = 1, cells(2);  do i = 1,cells1Red
     vectorField_fourier(1:3,i,j,k) = scalarField_fourier(i,j,k)*xi1st(1:3,i,j,k)                    ! ToDo: no -conjg?
-  enddo; enddo; enddo
+  end do; end do; end do
 
 end subroutine utilities_fourierScalarGradient
 
@@ -748,11 +787,9 @@ end subroutine utilities_fourierScalarGradient
 !--------------------------------------------------------------------------------------------------
 subroutine utilities_fourierVectorDivergence()
 
-  integer :: i, j, k
 
-  do k = 1, cells3;  do j = 1, cells(2);  do i = 1,grid1Red
+  scalarField_fourier(1:cells1Red,1:cells(2),1:cells3) = sum(vectorField_fourier(1:3,1:cells1Red,1:cells(2),1:cells3) &
-    scalarField_fourier(i,j,k) = sum(vectorField_fourier(1:3,i,j,k)*conjg(-xi1st(1:3,i,j,k)))
+                                                             *conjg(-xi1st),1)
-  enddo; enddo; enddo
 
 end subroutine utilities_fourierVectorDivergence
 
@@ -764,11 +801,12 @@ subroutine utilities_fourierVectorGradient()
 
   integer :: i, j, k, m, n
 
-  do k = 1, cells3;  do j = 1, cells(2);  do i = 1,grid1Red
+
+  do k = 1, cells3;  do j = 1, cells(2);  do i = 1,cells1Red
     do m = 1, 3; do n = 1, 3
       tensorField_fourier(m,n,i,j,k) = vectorField_fourier(m,i,j,k)*xi1st(n,i,j,k)
-    enddo; enddo
+    end do; end do
-  enddo; enddo; enddo
+  end do; end do; end do
 
 end subroutine utilities_fourierVectorGradient
 
@@ -780,9 +818,10 @@ subroutine utilities_fourierTensorDivergence()
 
   integer :: i, j, k
 
-  do k = 1, cells3;  do j = 1, cells(2);  do i = 1,grid1Red
+
+  do k = 1, cells3;  do j = 1, cells(2);  do i = 1,cells1Red
     vectorField_fourier(:,i,j,k) = matmul(tensorField_fourier(:,:,i,j,k),conjg(-xi1st(:,i,j,k)))
-  enddo; enddo; enddo
+  end do; end do; end do
 
 end subroutine utilities_fourierTensorDivergence
 
@@ -812,9 +851,9 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,&
 
   homogenization_F  = reshape(F,[3,3,product(cells(1:2))*cells3])                                   ! set materialpoint target F to estimated field
 
-  call homogenization_mechanical_response(Delta_t,[1,1],[1,product(cells(1:2))*cells3])             ! calculate P field
+  call homogenization_mechanical_response(Delta_t,1,product(cells(1:2))*cells3)                     ! calculate P field
   if (.not. terminallyIll) &
-    call homogenization_thermal_response(Delta_t,[1,1],[1,product(cells(1:2))*cells3])
+    call homogenization_thermal_response(Delta_t,1,product(cells(1:2))*cells3)
   if (.not. terminallyIll) &
     call homogenization_mechanical_response2(Delta_t,[1,1],[1,product(cells(1:2))*cells3])
 
@@ -884,11 +923,10 @@ pure function utilities_calculateRate(heterogeneous,field0,field,dt,avRate)
   real(pReal),             dimension(3,3,cells(1),cells(2),cells3) :: &
     utilities_calculateRate
 
-  if (heterogeneous) then
+
-    utilities_calculateRate = (field-field0) / dt
+  utilities_calculateRate = merge((field-field0) / dt, &
-  else
+                                  spread(spread(spread(avRate,3,cells(1)),4,cells(2)),5,cells3), &
-    utilities_calculateRate = spread(spread(spread(avRate,3,cells(1)),4,cells(2)),5,cells3)
+                                  heterogeneous)
-  endif
 
 end function utilities_calculateRate
 
@@ -980,6 +1018,7 @@ end function utilities_getFreqDerivative
 subroutine utilities_updateCoords(F)
 
   real(pReal),   dimension(3,3,cells(1),cells(2),cells3), intent(in) :: F
+
   real(pReal),   dimension(3,  cells(1),cells(2),cells3)             :: IPcoords
   real(pReal),   dimension(3,  cells(1),cells(2),cells3+2)           :: IPfluct_padded              ! Fluctuations of cell center displacement (padded along z for MPI)
   real(pReal),   dimension(3,  cells(1)+1,cells(2)+1,cells3+1)       :: nodeCoords
@@ -1010,20 +1049,23 @@ subroutine utilities_updateCoords(F)
                         1, 1, 1, &
                         0, 1, 1  ], [3,8])
 
+
   step = geomSize/real(cells, pReal)
  !--------------------------------------------------------------------------------------------------
  ! integration in Fourier space to get fluctuations of cell center discplacements
   tensorField_real(1:3,1:3,1:cells(1),1:cells(2),1:cells3) = F
   call utilities_FFTtensorForward()
 
-  do k = 1, cells3; do j = 1, cells(2); do i = 1, grid1Red
+  !$OMP PARALLEL DO
+  do k = 1, cells3; do j = 1, cells(2); do i = 1, cells1Red
     if (any([i,j,k+cells3Offset] /= 1)) then
       vectorField_fourier(1:3,i,j,k) = matmul(tensorField_fourier(1:3,1:3,i,j,k),xi2nd(1:3,i,j,k)) &
                                      / sum(conjg(-xi2nd(1:3,i,j,k))*xi2nd(1:3,i,j,k)) * cmplx(wgt,0.0,pReal)
     else
       vectorField_fourier(1:3,i,j,k) = cmplx(0.0,0.0,pReal)
-    endif
+    end if
-  enddo; enddo; enddo
+  end do; end do; end do
+  !$OMP END PARALLEL DO
 
   call fftw_mpi_execute_dft_c2r(planVectorBack,vectorField_fourier,vectorField_real)
 
@@ -1041,7 +1083,7 @@ subroutine utilities_updateCoords(F)
   rank_b = modulo(worldrank-1_MPI_INTEGER_KIND,worldsize)
 
   ! send bottom layer to process below
-  call MPI_Isend(IPfluct_padded(:,:,:,2),      c,MPI_DOUBLE,rank_b,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,request(1),err_MPI)
+  call MPI_Isend(IPfluct_padded(:,:,:,2),       c,MPI_DOUBLE,rank_b,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,request(1),err_MPI)
   if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
   call MPI_Irecv(IPfluct_padded(:,:,:,cells3+2),c,MPI_DOUBLE,rank_t,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,request(2),err_MPI)
   if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
@@ -1049,7 +1091,7 @@ subroutine utilities_updateCoords(F)
   ! send top layer to process above
   call MPI_Isend(IPfluct_padded(:,:,:,cells3+1),c,MPI_DOUBLE,rank_t,1_MPI_INTEGER_KIND,MPI_COMM_WORLD,request(3),err_MPI)
   if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
-  call MPI_Irecv(IPfluct_padded(:,:,:,1),      c,MPI_DOUBLE,rank_b,1_MPI_INTEGER_KIND,MPI_COMM_WORLD,request(4),err_MPI)
+  call MPI_Irecv(IPfluct_padded(:,:,:,1),       c,MPI_DOUBLE,rank_b,1_MPI_INTEGER_KIND,MPI_COMM_WORLD,request(4),err_MPI)
   if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
 
   call MPI_Waitall(4,request,status,err_MPI)

src/homogenization.f90

@@ -222,14 +222,13 @@ end subroutine homogenization_init
 !--------------------------------------------------------------------------------------------------
 !> @brief
 !--------------------------------------------------------------------------------------------------
-subroutine homogenization_mechanical_response(Delta_t,FEsolving_execIP,FEsolving_execElem)
+subroutine homogenization_mechanical_response(Delta_t,cell_start,cell_end)
 
   real(pReal), intent(in) :: Delta_t                                                                !< time increment
-  integer, dimension(2), intent(in) :: FEsolving_execElem, FEsolving_execIP
+  integer, intent(in) :: &
+    cell_start, cell_end
   integer :: &
     NiterationMPstate, &
-    ip, &                                                                                            !< integration point number
-    el, &                                                                                            !< element number
     co, ce, ho, en
   logical :: &
     converged
@@ -237,45 +236,42 @@ subroutine homogenization_mechanical_response(Delta_t,FEsolving_execIP,FEsolving
     doneAndHappy
 
 
-  !$OMP PARALLEL DO PRIVATE(ce,en,ho,NiterationMPstate,converged,doneAndHappy)
+  !$OMP PARALLEL DO PRIVATE(en,ho,co,NiterationMPstate,converged,doneAndHappy)
-  do el = FEsolving_execElem(1),FEsolving_execElem(2)
+  do ce = cell_start, cell_end
 
-    do ip = FEsolving_execIP(1),FEsolving_execIP(2)
+    en = material_homogenizationEntry(ce)
-      ce = (el-1)*discretization_nIPs + ip
+    ho = material_homogenizationID(ce)
-      en = material_homogenizationEntry(ce)
-      ho = material_homogenizationID(ce)
 
-      call phase_restore(ce,.false.) ! wrong name (is more a forward function)
+    call phase_restore(ce,.false.) ! wrong name (is more a forward function)
 
-      if(homogState(ho)%sizeState > 0)  homogState(ho)%state(:,en) = homogState(ho)%state0(:,en)
+    if(homogState(ho)%sizeState > 0)  homogState(ho)%state(:,en) = homogState(ho)%state0(:,en)
-      if(damageState_h(ho)%sizeState > 0) damageState_h(ho)%state(:,en) = damageState_h(ho)%state0(:,en)
+    if(damageState_h(ho)%sizeState > 0) damageState_h(ho)%state(:,en) = damageState_h(ho)%state0(:,en)
-      call damage_partition(ce)
+    call damage_partition(ce)
 
-      doneAndHappy = [.false.,.true.]
+    doneAndHappy = [.false.,.true.]
 
-      NiterationMPstate = 0
+    NiterationMPstate = 0
-      convergenceLooping: do while (.not. (terminallyIll .or. doneAndHappy(1)) &
+    convergenceLooping: do while (.not. (terminallyIll .or. doneAndHappy(1)) &
-                                    .and. NiterationMPstate < num%nMPstate)
+                                  .and. NiterationMPstate < num%nMPstate)
-        NiterationMPstate = NiterationMPstate + 1
+      NiterationMPstate = NiterationMPstate + 1
 
-        call mechanical_partition(homogenization_F(1:3,1:3,ce),ce)
+      call mechanical_partition(homogenization_F(1:3,1:3,ce),ce)
-        converged = all([(phase_mechanical_constitutive(Delta_t,co,ip,el),co=1,homogenization_Nconstituents(ho))])
+      converged = all([(phase_mechanical_constitutive(Delta_t,co,ce),co=1,homogenization_Nconstituents(ho))])
-        if (converged) then
+      if (converged) then
-          doneAndHappy = mechanical_updateState(Delta_t,homogenization_F(1:3,1:3,ce),ce)
+        doneAndHappy = mechanical_updateState(Delta_t,homogenization_F(1:3,1:3,ce),ce)
-          converged = all(doneAndHappy)
+        converged = all(doneAndHappy)
-        else
+      else
-          doneAndHappy = [.true.,.false.]
+        doneAndHappy = [.true.,.false.]
-        endif
+      end if
-      enddo convergenceLooping
+    end do convergenceLooping
 
-      converged = converged .and. all([(phase_damage_constitutive(Delta_t,co,ip,el),co=1,homogenization_Nconstituents(ho))])
+    converged = converged .and. all([(phase_damage_constitutive(Delta_t,co,ce),co=1,homogenization_Nconstituents(ho))])
 
-      if (.not. converged) then
+    if (.not. converged) then
-        if (.not. terminallyIll) print*, ' Cell ', ce, ' terminally ill'
+      if (.not. terminallyIll) print*, ' Cell ', ce, ' terminally ill'
-        terminallyIll = .true.
+      terminallyIll = .true.
-      endif
+    end if
-    enddo
+  end do
-  enddo
   !$OMP END PARALLEL DO
 
 end subroutine homogenization_mechanical_response
@@ -284,31 +280,27 @@ end subroutine homogenization_mechanical_response
 !--------------------------------------------------------------------------------------------------
 !> @brief
 !--------------------------------------------------------------------------------------------------
-subroutine homogenization_thermal_response(Delta_t,FEsolving_execIP,FEsolving_execElem)
+subroutine homogenization_thermal_response(Delta_t,cell_start,cell_end)
 
   real(pReal), intent(in) :: Delta_t                                                                !< time increment
-  integer, dimension(2), intent(in) :: FEsolving_execElem, FEsolving_execIP
+  integer, intent(in) :: &
+    cell_start, cell_end
   integer :: &
-    ip, &                                                                                           !< integration point number
-    el, &                                                                                           !< element number
     co, ce, ho
 
 
-  !$OMP PARALLEL DO PRIVATE(ho,ce)
+  !$OMP PARALLEL DO PRIVATE(ho)
-  do el = FEsolving_execElem(1),FEsolving_execElem(2)
+  do ce = cell_start, cell_end
     if (terminallyIll) continue
-    do ip = FEsolving_execIP(1),FEsolving_execIP(2)
+    ho = material_homogenizationID(ce)
-      ce = (el-1)*discretization_nIPs + ip
+    call thermal_partition(ce)
-      ho = material_homogenizationID(ce)
+    do co = 1, homogenization_Nconstituents(ho)
-      call thermal_partition(ce)
+      if (.not. phase_thermal_constitutive(Delta_t,material_phaseID(co,ce),material_phaseEntry(co,ce))) then
-      do co = 1, homogenization_Nconstituents(ho)
+        if (.not. terminallyIll) print*, ' Cell ', ce, ' terminally ill'
-        if (.not. phase_thermal_constitutive(Delta_t,material_phaseID(co,ce),material_phaseEntry(co,ce))) then
+        terminallyIll = .true.
-          if (.not. terminallyIll) print*, ' Cell ', ce, ' terminally ill'
+      end if
-          terminallyIll = .true.
+    end do
-        endif
+  end do
-      enddo
-    enddo
-  enddo
   !$OMP END PARALLEL DO
 
 end subroutine homogenization_thermal_response
@@ -331,11 +323,11 @@ subroutine homogenization_mechanical_response2(Delta_t,FEsolving_execIP,FEsolvin
   elementLooping3: do el = FEsolving_execElem(1),FEsolving_execElem(2)
     IpLooping3: do ip = FEsolving_execIP(1),FEsolving_execIP(2)
       ce = (el-1)*discretization_nIPs + ip
-      ho = material_homogenizationID(ce)
+    ho = material_homogenizationID(ce)
-      do co = 1, homogenization_Nconstituents(ho)
+    do co = 1, homogenization_Nconstituents(ho)
-        call crystallite_orientations(co,ip,el)
+      call crystallite_orientations(co,ip,el)
       enddo
-      call mechanical_homogenize(Delta_t,ce)
+    call mechanical_homogenize(Delta_t,ce)
     enddo IpLooping3
   enddo elementLooping3
   !$OMP END PARALLEL DO

src/homogenization_mechanical.f90

@@ -131,20 +131,15 @@ module subroutine mechanical_homogenize(Delta_t,ce)
   integer :: co
 
 
-  homogenization_P(1:3,1:3,ce)            = phase_P(1,ce)
+  homogenization_P(1:3,1:3,ce)            = phase_P(1,ce)*material_v(1,ce)
-  homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = phase_mechanical_dPdF(Delta_t,1,ce)
+  homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = phase_mechanical_dPdF(Delta_t,1,ce)*material_v(1,ce)
   do co = 2, homogenization_Nconstituents(material_homogenizationID(ce))
     homogenization_P(1:3,1:3,ce)            = homogenization_P(1:3,1:3,ce) &
-                                            + phase_P(co,ce)
+                                            + phase_P(co,ce)*material_v(co,ce)
     homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = homogenization_dPdF(1:3,1:3,1:3,1:3,ce) &
-                                            + phase_mechanical_dPdF(Delta_t,co,ce)
+                                            + phase_mechanical_dPdF(Delta_t,co,ce)*material_v(co,ce)
   end do
 
-  homogenization_P(1:3,1:3,ce)            = homogenization_P(1:3,1:3,ce) &
-                                          / real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal)
-  homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = homogenization_dPdF(1:3,1:3,1:3,1:3,ce) &
-                                          / real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal)
-
 end subroutine mechanical_homogenize
 
 

src/homogenization_thermal.f90

@@ -105,13 +105,11 @@ module function homogenization_mu_T(ce) result(mu)
   integer :: co
 
 
-  mu = phase_mu_T(1,ce)
+  mu = phase_mu_T(1,ce)*material_v(1,ce)
   do co = 2, homogenization_Nconstituents(material_homogenizationID(ce))
-    mu = mu + phase_mu_T(co,ce)
+    mu = mu + phase_mu_T(co,ce)*material_v(co,ce)
   end do
 
-  mu = mu / real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal)
-
 end function homogenization_mu_T
 
 
@@ -126,13 +124,11 @@ module function homogenization_K_T(ce) result(K)
   integer :: co
 
 
-  K = phase_K_T(1,ce)
+  K = phase_K_T(1,ce)*material_v(1,ce)
   do co = 2, homogenization_Nconstituents(material_homogenizationID(ce))
-    K = K + phase_K_T(co,ce)
+    K = K + phase_K_T(co,ce)*material_v(co,ce)
   end do
 
-  K = K / real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal)
-
 end function homogenization_K_T
 
 
@@ -147,13 +143,11 @@ module function homogenization_f_T(ce) result(f)
   integer :: co
 
 
-  f = phase_f_T(material_phaseID(1,ce),material_phaseEntry(1,ce))
+  f = phase_f_T(material_phaseID(1,ce),material_phaseEntry(1,ce))*material_v(1,ce)
   do co = 2, homogenization_Nconstituents(material_homogenizationID(ce))
-    f = f + phase_f_T(material_phaseID(co,ce),material_phaseEntry(co,ce))
+    f = f + phase_f_T(material_phaseID(co,ce),material_phaseEntry(co,ce))*material_v(co,ce)
   end do
 
-  f = f/real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal)
-
 end function homogenization_f_T
 
 

src/material.f90

@@ -17,16 +17,17 @@ module material
   implicit none
   private
 
-  type :: tRotationContainer
+  type, public :: tRotationContainer
-    type(tRotation), dimension(:),  allocatable :: data
+    type(tRotation), dimension(:), allocatable :: data
-  end type
+  end type tRotationContainer
-  type :: tTensorContainer
+
+  type, public :: tTensorContainer
     real(pReal), dimension(:,:,:), allocatable :: data
-  end type
+  end type tTensorContainer
 
 
-  type(tRotationContainer), dimension(:), allocatable :: material_O_0
+  type(tRotationContainer), dimension(:), allocatable, public, protected :: material_O_0
-  type(tTensorContainer),   dimension(:), allocatable :: material_F_i_0
+  type(tTensorContainer),   dimension(:), allocatable, public, protected :: material_F_i_0
 
   integer, dimension(:), allocatable, public, protected :: &
     homogenization_Nconstituents                                                                    !< number of grains in each homogenization
@@ -37,20 +38,17 @@ module material
     material_name_phase, &                                                                          !< name of each phase
     material_name_homogenization                                                                    !< name of each homogenization
 
-  integer, dimension(:),     allocatable, public, protected :: &                                    ! (elem)
+  integer, dimension(:),   allocatable, public, protected :: &                                      ! (cell)
-    material_homogenizationID, &                                                                    !< per cell TODO: material_ID_homogenization
+    material_homogenizationID, &                                                                    ! TODO: rename to material_ID_homogenization
-    material_homogenizationEntry                                                                    !< per cell TODO: material_entry_homogenization
+    material_homogenizationEntry                                                                    ! TODO: rename to material_entry_homogenization
-  integer, dimension(:,:),   allocatable, public, protected :: &                                    ! (constituent,elem)
+  integer, dimension(:,:), allocatable, public, protected :: &                                      ! (constituent,cell)
-    material_phaseAt, &                                                                             !< phase ID of each element TODO: remove
+    material_phaseID, &                                                                             ! TODO: rename to material_ID_phase
-    material_phaseID, &                                                                             !< per (constituent,cell) TODO: material_ID_phase
+    material_phaseEntry                                                                             ! TODO: rename to material_entry_phase
-    material_phaseEntry                                                                             !< per (constituent,cell) TODO: material_entry_phase
+
-  integer, dimension(:,:,:), allocatable, public, protected :: &                                    ! (constituent,IP,elem)
+  real(pReal), dimension(:,:), allocatable, public, protected :: &
-    material_phaseMemberAt                                                                          !TODO: remove
+    material_v                                                                                      ! fraction
+
   public :: &
-    tTensorContainer, &
-    tRotationContainer, &
-    material_F_i_0, &
-    material_O_0, &
     material_init
 
 contains
@@ -97,11 +95,12 @@ subroutine parse()
     counterPhase, &
     counterHomogenization
 
-  real(pReal) :: &
+  real(pReal) :: v
-    frac
   integer :: &
-    el, ip, co, ma, &
+    el, ip, &
-    h, ce
+    ho, ph, &
+    co, ce, &
+    ma
 
   materials       => config_material%get('material')
   phases          => config_material%get('phase')
@@ -118,51 +117,52 @@ subroutine parse()
 #endif
 
   allocate(homogenization_Nconstituents(homogenizations%length))
-  do h=1, homogenizations%length
+  do ho=1, homogenizations%length
-    homogenization => homogenizations%get(h)
+    homogenization => homogenizations%get(ho)
-    homogenization_Nconstituents(h) = homogenization%get_asInt('N_constituents')
+    homogenization_Nconstituents(ho) = homogenization%get_asInt('N_constituents')
   end do
   homogenization_maxNconstituents = maxval(homogenization_Nconstituents)
 
   allocate(counterPhase(phases%length),source=0)
   allocate(counterHomogenization(homogenizations%length),source=0)
 
-  allocate(material_phaseAt(homogenization_maxNconstituents,discretization_Nelems),source=0)
+  allocate(material_homogenizationID(discretization_Ncells),source=0)
-  allocate(material_phaseMemberAt(homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems),source=0)
+  allocate(material_homogenizationEntry(discretization_Ncells),source=0)
 
+  allocate(material_phaseID(homogenization_maxNconstituents,discretization_Ncells),source=0)
+  allocate(material_phaseEntry(homogenization_maxNconstituents,discretization_Ncells),source=0)
 
-  allocate(material_homogenizationID(discretization_nIPs*discretization_Nelems),source=0)
+  allocate(material_v(homogenization_maxNconstituents,discretization_Ncells),source=0.0_pReal)
-  allocate(material_homogenizationEntry(discretization_nIPs*discretization_Nelems),source=0)
-  allocate(material_phaseID(homogenization_maxNconstituents,discretization_nIPs*discretization_Nelems),source=0)
-  allocate(material_phaseEntry(homogenization_maxNconstituents,discretization_nIPs*discretization_Nelems),source=0)
 
   do el = 1, discretization_Nelems
-    material     => materials%get(discretization_materialAt(el))
+    material => materials%get(discretization_materialAt(el))
-    constituents => material%get('constituents')
 
+    ho = homogenizations%getIndex(material%get_asString('homogenization'))
     do ip = 1, discretization_nIPs
       ce = (el-1)*discretization_nIPs + ip
-      material_homogenizationID(ce) = homogenizations%getIndex(material%get_asString('homogenization'))
+      material_homogenizationID(ce) = ho
-      counterHomogenization(material_homogenizationID(ce)) = counterHomogenization(material_homogenizationID(ce)) + 1
+      counterHomogenization(ho) = counterHomogenization(ho) + 1
-      material_homogenizationEntry(ce) = counterHomogenization(material_homogenizationID(ce))
+      material_homogenizationEntry(ce) = counterHomogenization(ho)
     end do
 
-    frac = 0.0_pReal
+    constituents => material%get('constituents')
     do co = 1, constituents%length
       constituent => constituents%get(co)
-      frac = frac + constituent%get_asFloat('v')
 
-      material_phaseAt(co,el) = phases%getIndex(constituent%get_asString('phase'))
+      v = constituent%get_asFloat('v')
+
+      ph = phases%getIndex(constituent%get_asString('phase'))
       do ip = 1, discretization_nIPs
         ce = (el-1)*discretization_nIPs + ip
-        counterPhase(material_phaseAt(co,el)) = counterPhase(material_phaseAt(co,el)) + 1
+        material_phaseID(co,ce) = ph
-        material_phaseMemberAt(co,ip,el)      = counterPhase(material_phaseAt(co,el))
+        counterPhase(ph) = counterPhase(ph) + 1
-        material_phaseEntry(co,ce) = counterPhase(material_phaseAt(co,el))
+        material_phaseEntry(co,ce) = counterPhase(ph)
-        material_phaseID(co,ce) = material_phaseAt(co,el)
+        material_v(co,ce) = v
       end do
 
     end do
-    if (dNeq(frac,1.0_pReal,1.e-12_pReal)) call IO_error(153,ext_msg='constituent')
+    if (dNeq(sum(material_v(1:constituents%length,ce)),1.0_pReal,1.e-9_pReal)) &
+      call IO_error(153,ext_msg='constituent')
 
   end do
 

src/math.f90

@@ -262,9 +262,8 @@ pure function math_identity4th()
     math_identity4th(i,j,k,l) = 0.5_pReal*(math_I3(i,k)*math_I3(j,l)+math_I3(i,l)*math_I3(j,k))
   enddo
 #else
-  do i=1,3; do j=1,3; do k=1,3; do l=1,3
+  forall(i=1:3, j=1:3, k=1:3, l=1:3) &
     math_identity4th(i,j,k,l) = 0.5_pReal*(math_I3(i,k)*math_I3(j,l)+math_I3(i,l)*math_I3(j,k))
-  enddo; enddo; enddo; enddo
 #endif
 
 end function math_identity4th
@@ -338,9 +337,7 @@ pure function math_outer(A,B)
     math_outer(i,j) = A(i)*B(j)
   enddo
 #else
-  do i=1,size(A,1); do j=1,size(B,1)
+  forall(i=1:size(A,1), j=1:size(B,1)) math_outer(i,j) = A(i)*B(j)
-    math_outer(i,j) = A(i)*B(j)
-  enddo; enddo
 #endif
 
 end function math_outer
@@ -387,9 +384,7 @@ pure function math_mul3333xx33(A,B)
     math_mul3333xx33(i,j) = sum(A(i,j,1:3,1:3)*B(1:3,1:3))
   enddo
 #else
-  do i=1,3; do j=1,3
+  forall (i=1:3, j=1:3) math_mul3333xx33(i,j) = sum(A(i,j,1:3,1:3)*B(1:3,1:3))
-    math_mul3333xx33(i,j) = sum(A(i,j,1:3,1:3)*B(1:3,1:3))
-  enddo; enddo
 #endif
 
 end function math_mul3333xx33
@@ -411,9 +406,7 @@ pure function math_mul3333xx3333(A,B)
     math_mul3333xx3333(i,j,k,l) = sum(A(i,j,1:3,1:3)*B(1:3,1:3,k,l))
   enddo
 #else
-  do i=1,3; do j=1,3; do k=1,3; do l=1,3
+  forall(i=1:3, j=1:3, k=1:3, l=1:3) math_mul3333xx3333(i,j,k,l) = sum(A(i,j,1:3,1:3)*B(1:3,1:3,k,l))
-    math_mul3333xx3333(i,j,k,l) = sum(A(i,j,1:3,1:3)*B(1:3,1:3,k,l))
-  enddo; enddo; enddo; enddo
 #endif
 
 end function math_mul3333xx3333
@@ -752,9 +745,7 @@ pure function math_3333to99(m3333)
     math_3333to99(i,j) = m3333(MAPPLAIN(1,i),MAPPLAIN(2,i),MAPPLAIN(1,j),MAPPLAIN(2,j))
   enddo
 #else
-  do i=1,9; do j=1,9
+  forall(i=1:9, j=1:9) math_3333to99(i,j) = m3333(MAPPLAIN(1,i),MAPPLAIN(2,i),MAPPLAIN(1,j),MAPPLAIN(2,j))
-    math_3333to99(i,j) = m3333(MAPPLAIN(1,i),MAPPLAIN(2,i),MAPPLAIN(1,j),MAPPLAIN(2,j))
-  enddo; enddo
 #endif
 
 end function math_3333to99
@@ -775,9 +766,7 @@ pure function math_99to3333(m99)
     math_99to3333(MAPPLAIN(1,i),MAPPLAIN(2,i),MAPPLAIN(1,j),MAPPLAIN(2,j)) = m99(i,j)
   enddo
 #else
-  do i=1,9; do j=1,9
+  forall(i=1:9, j=1:9) math_99to3333(MAPPLAIN(1,i),MAPPLAIN(2,i),MAPPLAIN(1,j),MAPPLAIN(2,j)) = m99(i,j)
-    math_99to3333(MAPPLAIN(1,i),MAPPLAIN(2,i),MAPPLAIN(1,j),MAPPLAIN(2,j)) = m99(i,j)
-  enddo; enddo
 #endif
 
 end function math_99to3333
@@ -810,9 +799,7 @@ pure function math_sym3333to66(m3333,weighted)
     math_sym3333to66(i,j) = w(i)*w(j)*m3333(MAPNYE(1,i),MAPNYE(2,i),MAPNYE(1,j),MAPNYE(2,j))
   enddo
 #else
-  do i=1,6; do j=1,6
+  forall(i=1:6, j=1:6) math_sym3333to66(i,j) = w(i)*w(j)*m3333(MAPNYE(1,i),MAPNYE(2,i),MAPNYE(1,j),MAPNYE(2,j))
-    math_sym3333to66(i,j) = w(i)*w(j)*m3333(MAPNYE(1,i),MAPNYE(2,i),MAPNYE(1,j),MAPNYE(2,j))
-  enddo; enddo
 #endif
 
 end function math_sym3333to66
@@ -950,9 +937,7 @@ pure function math_3333toVoigt66_stiffness(C) result(C_tilde)
     C_tilde(i,j) = C(MAPVOIGT(1,i),MAPVOIGT(2,i),MAPVOIGT(1,j),MAPVOIGT(2,j))
   end do
 #else
-  do i=1,6; do j=1,6
+  forall(i=1:6, j=1:6) C_tilde(i,j) = C(MAPVOIGT(1,i),MAPVOIGT(2,i),MAPVOIGT(1,j),MAPVOIGT(2,j))
-    C_tilde(i,j) = C(MAPVOIGT(1,i),MAPVOIGT(2,i),MAPVOIGT(1,j),MAPVOIGT(2,j))
-  end do; end do
 #endif
 
 end function math_3333toVoigt66_stiffness

src/mesh/FEM_utilities.f90

@@ -150,7 +150,7 @@ subroutine utilities_constitutiveResponse(timeinc,P_av,forwardData)
 
   print'(/,1x,a)', '... evaluating constitutive response ......................................'
 
-  call homogenization_mechanical_response(timeinc,[1,mesh_maxNips],[1,mesh_NcpElems])                ! calculate P field
+  call homogenization_mechanical_response(timeinc,1,mesh_maxNips*mesh_NcpElems)                     ! calculate P field
   if (.not. terminallyIll) &
     call homogenization_mechanical_response2(timeinc,[1,mesh_maxNips],[1,mesh_NcpElems])
   cutBack = .false.

src/phase.f90

@@ -228,15 +228,15 @@ module phase
 
     end function phase_thermal_constitutive
 
-    module function phase_damage_constitutive(Delta_t,co,ip,el) result(converged_)
+    module function phase_damage_constitutive(Delta_t,co,ce) result(converged_)
       real(pReal), intent(in) :: Delta_t
-      integer, intent(in) :: co, ip, el
+      integer, intent(in) :: co, ce
       logical :: converged_
     end function phase_damage_constitutive
 
-    module function phase_mechanical_constitutive(Delta_t,co,ip,el) result(converged_)
+    module function phase_mechanical_constitutive(Delta_t,co,ce) result(converged_)
       real(pReal), intent(in) :: Delta_t
-      integer, intent(in) :: co, ip, el
+      integer, intent(in) :: co, ce
       logical :: converged_
     end function phase_mechanical_constitutive
 
@@ -264,19 +264,18 @@ module phase
       real(pReal) :: f
     end function phase_f_T
 
-    module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e)
+    module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,ip,el)
       integer, intent(in) :: &
         ph, &
-        i, &
+        ip, &
-        e
+        el
         type(tRotationContainer), dimension(:), intent(in) :: orientation
     end subroutine plastic_nonlocal_updateCompatibility
 
-    module subroutine plastic_dependentState(co,ip,el)
+    module subroutine plastic_dependentState(ph,en)
       integer, intent(in) :: &
-        co, &                                                                                       !< component-ID of integration point
+        ph, &
-        ip, &                                                                                       !< integration point
+        en
-        el                                                                                          !< element
     end subroutine plastic_dependentState
 
 
@@ -503,8 +502,8 @@ subroutine crystallite_init()
     ce, &
     co, &                                                                                           !< counter in integration point component loop
     ip, &                                                                                           !< counter in integration point loop
-    el                                                                                              !< counter in element loop
+    el, &                                                                                           !< counter in element loop
-
+    en, ph
   class(tNode), pointer :: &
     num_crystallite, &
     phases
@@ -543,13 +542,15 @@ subroutine crystallite_init()
 
   phases => config_material%get('phase')
 
-  !$OMP PARALLEL DO PRIVATE(ce)
+  !$OMP PARALLEL DO PRIVATE(ce,ph,en)
   do el = 1, discretization_Nelems
     do ip = 1, discretization_nIPs
       ce = (el-1)*discretization_nIPs + ip
       do co = 1,homogenization_Nconstituents(material_homogenizationID(ce))
+        en = material_phaseEntry(co,ce)
+        ph = material_phaseID(co,ce)
         call crystallite_orientations(co,ip,el)
-        call plastic_dependentState(co,ip,el)                                                       ! update dependent state variables to be consistent with basic states
+        call plastic_dependentState(ph,en)                                                          ! update dependent state variables to be consistent with basic states
      end do
     end do
   end do
@@ -577,8 +578,8 @@ subroutine crystallite_orientations(co,ip,el)
 
   call phase_O(ph)%data(en)%fromMatrix(transpose(math_rotationalPart(mechanical_F_e(ph,en))))
 
-  if (plasticState(material_phaseAt(1,el))%nonlocal) &
+  if (plasticState(material_phaseID(1,(el-1)*discretization_nIPs + ip))%nonlocal) &
-    call plastic_nonlocal_updateCompatibility(phase_O,material_phaseAt(1,el),ip,el)
+    call plastic_nonlocal_updateCompatibility(phase_O,material_phaseID(1,(el-1)*discretization_nIPs + ip),ip,el)
 
 
 end subroutine crystallite_orientations

src/phase_damage.f90

@@ -127,21 +127,20 @@ end subroutine damage_init
 !--------------------------------------------------------------------------------------------------
 !> @brief calculate stress (P)
 !--------------------------------------------------------------------------------------------------
-module function phase_damage_constitutive(Delta_t,co,ip,el) result(converged_)
+module function phase_damage_constitutive(Delta_t,co,ce) result(converged_)
 
   real(pReal), intent(in) :: Delta_t
   integer, intent(in) :: &
     co, &
-    ip, &
+    ce
-    el
   logical :: converged_
 
   integer :: &
     ph, en
 
 
-  ph = material_phaseID(co,(el-1)*discretization_nIPs + ip)
+  ph = material_phaseID(co,ce)
-  en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip)
+  en = material_phaseEntry(co,ce)
 
   converged_ = .not. integrateDamageState(Delta_t,ph,en)
 

src/phase_mechanical.f90

@@ -79,15 +79,12 @@ submodule(phase) mechanical
         en
     end subroutine plastic_isotropic_LiAndItsTangent
 
-    module function plastic_dotState(subdt,co,ip,el,ph,en) result(dotState)
+    module function plastic_dotState(subdt,ph,en) result(dotState)
       integer, intent(in) :: &
-        co, &                                                                                       !< constituent
-        ip, &                                                                                       !< integration point
-        el, &                                                                                       !< element
         ph, &
         en
       real(pReal),  intent(in) :: &
-        subdt                                                                                       !< timestep
+        subdt                                                                                           !< timestep
       real(pReal), dimension(plasticState(ph)%sizeDotState) :: &
         dotState
     end function plastic_dotState
@@ -365,13 +362,11 @@ end subroutine mechanical_results
 !> @brief calculation of stress (P) with time integration based on a residuum in Lp and
 !> intermediate acceleration of the Newton-Raphson correction
 !--------------------------------------------------------------------------------------------------
-function integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) result(broken)
+function integrateStress(F,subFp0,subFi0,Delta_t,ph,en) result(broken)
 
   real(pReal), dimension(3,3), intent(in) :: F,subFp0,subFi0
   real(pReal),                 intent(in) :: Delta_t
-  integer, intent(in)::         el, &                                                               ! element index
+  integer, intent(in) :: ph, en
-                                      ip, &                                                         ! integration point index
-                                      co                                                            ! grain index
 
   real(pReal), dimension(3,3)::       Fp_new, &                                                     ! plastic deformation gradient at end of timestep
                                       invFp_new, &                                                  ! inverse of Fp_new
@@ -419,19 +414,13 @@ function integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) result(broken)
                                       ierr, &                                                       ! error indicator for LAPACK
                                       o, &
                                       p, &
-                                      ph, &
-                                      en, &
                                       jacoCounterLp, &
                                       jacoCounterLi                                                 ! counters to check for Jacobian update
   logical :: error,broken
 
 
   broken = .true.
-
+  call plastic_dependentState(ph,en)
-  ph = material_phaseID(co,(el-1)*discretization_nIPs + ip)
-  en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip)
-
-  call plastic_dependentState(co,ip,el)
 
   Lpguess = phase_mechanical_Lp(ph)%data(1:3,1:3,en)                                              ! take as first guess
   Liguess = phase_mechanical_Li(ph)%data(1:3,1:3,en)                                              ! take as first guess
@@ -579,37 +568,32 @@ end function integrateStress
 !> @brief integrate stress, state with adaptive 1st order explicit Euler method
 !> using Fixed Point Iteration to adapt the stepsize
 !--------------------------------------------------------------------------------------------------
-function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) result(broken)
+function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en) result(broken)
 
   real(pReal), intent(in),dimension(3,3) :: F_0,F,subFp0,subFi0
   real(pReal), intent(in),dimension(:)   :: subState0
   real(pReal), intent(in) :: Delta_t
   integer, intent(in) :: &
-    el, &                                                                                           !< element index in element loop
+    ph, &
-    ip, &                                                                                           !< integration point index in ip loop
+    en
-    co                                                                                              !< grain index in grain loop
   logical :: &
     broken
 
   integer :: &
     NiterationState, &                                                                              !< number of iterations in state loop
-    ph, &
-    en, &
     sizeDotState
   real(pReal) :: &
     zeta
-  real(pReal), dimension(plasticState(material_phaseID(co,(el-1)*discretization_nIPs+ip))%sizeDotState) :: &
+  real(pReal), dimension(plasticState(ph)%sizeDotState) :: &
     r, &                                                                                            ! state residuum
     dotState
-  real(pReal), dimension(plasticState(material_phaseID(co,(el-1)*discretization_nIPs+ip))%sizeDotState,2) :: &
+  real(pReal), dimension(plasticState(ph)%sizeDotState,2) :: &
     dotState_last
 
 
-  ph = material_phaseID(co,(el-1)*discretization_nIPs + ip)
-  en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip)
   broken = .true.
 
-  dotState = plastic_dotState(Delta_t, co,ip,el,ph,en)
+  dotState = plastic_dotState(Delta_t,ph,en)
   if (any(IEEE_is_NaN(dotState))) return
 
   sizeDotState = plasticState(ph)%sizeDotState
@@ -620,10 +604,10 @@ function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) resul
     dotState_last(1:sizeDotState,2) = merge(dotState_last(1:sizeDotState,1),0.0, nIterationState > 1)
     dotState_last(1:sizeDotState,1) = dotState
 
-    broken = integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el)
+    broken = integrateStress(F,subFp0,subFi0,Delta_t,ph,en)
     if(broken) exit iteration
 
-    dotState = plastic_dotState(Delta_t, co,ip,el,ph,en)
+    dotState = plastic_dotState(Delta_t,ph,en)
     if (any(IEEE_is_NaN(dotState))) exit iteration
 
     zeta = damper(dotState,dotState_last(1:sizeDotState,1),dotState_last(1:sizeDotState,2))
@@ -672,31 +656,26 @@ end function integrateStateFPI
 !--------------------------------------------------------------------------------------------------
 !> @brief integrate state with 1st order explicit Euler method
 !--------------------------------------------------------------------------------------------------
-function integrateStateEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) result(broken)
+function integrateStateEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en) result(broken)
 
   real(pReal), intent(in),dimension(3,3) :: F_0,F,subFp0,subFi0
   real(pReal), intent(in),dimension(:)   :: subState0
   real(pReal), intent(in) :: Delta_t
   integer, intent(in) :: &
-    el, &                                                                                            !< element index in element loop
+    ph, &
-    ip, &                                                                                            !< integration point index in ip loop
+    en                                                                                               !< grain index in grain loop
-    co                                                                                               !< grain index in grain loop
   logical :: &
     broken
 
-  real(pReal), dimension(plasticState(material_phaseID(co,(el-1)*discretization_nIPs+ip))%sizeDotState) :: &
+  real(pReal), dimension(plasticState(ph)%sizeDotState) :: &
     dotState
   integer :: &
-    ph, &
-    en, &
     sizeDotState
 
 
-  ph = material_phaseID(co,(el-1)*discretization_nIPs + ip)
-  en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip)
   broken = .true.
 
-  dotState = plastic_dotState(Delta_t, co,ip,el,ph,en)
+  dotState = plastic_dotState(Delta_t,ph,en)
   if (any(IEEE_is_NaN(dotState))) return
 
   sizeDotState = plasticState(ph)%sizeDotState
@@ -706,7 +685,7 @@ function integrateStateEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) res
   broken = plastic_deltaState(ph,en)
   if(broken) return
 
-  broken = integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el)
+  broken = integrateStress(F,subFp0,subFi0,Delta_t,ph,en)
 
 end function integrateStateEuler
 
@@ -714,32 +693,27 @@ end function integrateStateEuler
 !--------------------------------------------------------------------------------------------------
 !> @brief integrate stress, state with 1st order Euler method with adaptive step size
 !--------------------------------------------------------------------------------------------------
-function integrateStateAdaptiveEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) result(broken)
+function integrateStateAdaptiveEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en) result(broken)
 
   real(pReal), intent(in),dimension(3,3) :: F_0,F,subFp0,subFi0
   real(pReal), intent(in),dimension(:)   :: subState0
   real(pReal), intent(in) :: Delta_t
   integer, intent(in) :: &
-    el, &                                                                                            !< element index in element loop
+    ph, &
-    ip, &                                                                                            !< integration point index in ip loop
+    en
-    co                                                                                               !< grain index in grain loop
   logical :: &
     broken
 
   integer :: &
-    ph, &
-    en, &
     sizeDotState
-  real(pReal), dimension(plasticState(material_phaseID(co,(el-1)*discretization_nIPs+ip))%sizeDotState) :: &
+  real(pReal), dimension(plasticState(ph)%sizeDotState) :: &
     r, &
     dotState
 
 
-  ph = material_phaseID(co,(el-1)*discretization_nIPs + ip)
-  en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip)
   broken = .true.
 
-  dotState = plastic_dotState(Delta_t, co,ip,el,ph,en)
+  dotState = plastic_dotState(Delta_t,ph,en)
   if (any(IEEE_is_NaN(dotState))) return
 
   sizeDotState = plasticState(ph)%sizeDotState
@@ -751,10 +725,10 @@ function integrateStateAdaptiveEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip
   broken = plastic_deltaState(ph,en)
   if(broken) return
 
-  broken = integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el)
+  broken = integrateStress(F,subFp0,subFi0,Delta_t,ph,en)
   if(broken) return
 
-  dotState = plastic_dotState(Delta_t, co,ip,el,ph,en)
+  dotState = plastic_dotState(Delta_t,ph,en)
   if (any(IEEE_is_NaN(dotState))) return
 
   broken = .not. converged(r + 0.5_pReal * dotState * Delta_t, &
@@ -767,12 +741,12 @@ end function integrateStateAdaptiveEuler
 !---------------------------------------------------------------------------------------------------
 !> @brief Integrate state (including stress integration) with the classic Runge Kutta method
 !---------------------------------------------------------------------------------------------------
-function integrateStateRK4(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) result(broken)
+function integrateStateRK4(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en) result(broken)
 
   real(pReal), intent(in),dimension(3,3) :: F_0,F,subFp0,subFi0
   real(pReal), intent(in),dimension(:)   :: subState0
   real(pReal), intent(in) :: Delta_t
-  integer, intent(in) :: co,ip,el
+  integer, intent(in) :: ph, en
   logical :: broken
 
   real(pReal), dimension(3,3), parameter :: &
@@ -787,7 +761,7 @@ function integrateStateRK4(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) resul
     B = [1.0_pReal/6.0_pReal, 1.0_pReal/3.0_pReal, 1.0_pReal/3.0_pReal, 1.0_pReal/6.0_pReal]
 
 
-  broken = integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el,A,B,C)
+  broken = integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en,A,B,C)
 
 end function integrateStateRK4
 
@@ -795,12 +769,12 @@ end function integrateStateRK4
 !---------------------------------------------------------------------------------------------------
 !> @brief Integrate state (including stress integration) with the Cash-Carp method
 !---------------------------------------------------------------------------------------------------
-function integrateStateRKCK45(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) result(broken)
+function integrateStateRKCK45(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en) result(broken)
 
   real(pReal), intent(in),dimension(3,3) :: F_0,F,subFp0,subFi0
   real(pReal), intent(in),dimension(:)   :: subState0
   real(pReal), intent(in) :: Delta_t
-  integer, intent(in) :: co,ip,el
+  integer, intent(in) :: ph, en
   logical :: broken
 
   real(pReal), dimension(5,5), parameter :: &
@@ -822,7 +796,7 @@ function integrateStateRKCK45(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) re
       13525.0_pReal/55296.0_pReal, 277.0_pReal/14336.0_pReal,  1._pReal/4._pReal]
 
 
-  broken = integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el,A,B,C,DB)
+  broken = integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en,A,B,C,DB)
 
 end function integrateStateRKCK45
 
@@ -831,7 +805,7 @@ end function integrateStateRKCK45
 !> @brief Integrate state (including stress integration) with an explicit Runge-Kutta method or an
 !! embedded explicit Runge-Kutta method
 !--------------------------------------------------------------------------------------------------
-function integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el,A,B,C,DB) result(broken)
+function integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en,A,B,C,DB) result(broken)
 
   real(pReal), intent(in),dimension(3,3) :: F_0,F,subFp0,subFi0
   real(pReal), intent(in),dimension(:)   :: subState0
@@ -840,28 +814,23 @@ function integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el,A,B,C,D
   real(pReal), dimension(:),   intent(in) :: B, C
   real(pReal), dimension(:),   intent(in), optional :: DB
   integer, intent(in) :: &
-    el, &                                                                                            !< element index in element loop
+    ph, &
-    ip, &                                                                                            !< integration point index in ip loop
+    en
-    co                                                                                               !< grain index in grain loop
   logical :: broken
 
   integer :: &
     stage, &                                                                                        ! stage index in integration stage loop
     n, &
-    ph, &
-    en, &
     sizeDotState
-  real(pReal), dimension(plasticState(material_phaseID(co,(el-1)*discretization_nIPs+ip))%sizeDotState) :: &
+  real(pReal), dimension(plasticState(ph)%sizeDotState) :: &
     dotState
-  real(pReal), dimension(plasticState(material_phaseID(co,(el-1)*discretization_nIPs+ip))%sizeDotState,size(B)) :: &
+  real(pReal), dimension(plasticState(ph)%sizeDotState,size(B)) :: &
     plastic_RKdotState
 
 
-  ph = material_phaseID(co,(el-1)*discretization_nIPs + ip)
-  en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip)
   broken = .true.
 
-  dotState = plastic_dotState(Delta_t, co,ip,el,ph,en)
+  dotState = plastic_dotState(Delta_t,ph,en)
   if (any(IEEE_is_NaN(dotState))) return
 
   sizeDotState = plasticState(ph)%sizeDotState
@@ -879,10 +848,10 @@ function integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el,A,B,C,D
     plasticState(ph)%state(1:sizeDotState,en) = subState0 &
                                               + dotState * Delta_t
 
-    broken = integrateStress(F_0 + (F-F_0) * Delta_t*C(stage),subFp0,subFi0,Delta_t*C(stage),co,ip,el)
+    broken = integrateStress(F_0+(F-F_0)*Delta_t*C(stage),subFp0,subFi0,Delta_t*C(stage), ph,en)
     if(broken) exit
 
-    dotState = plastic_dotState(Delta_t*C(stage), co,ip,el,ph,en)
+    dotState = plastic_dotState(Delta_t*C(stage), ph,en)
     if (any(IEEE_is_NaN(dotState))) exit
 
   enddo
@@ -904,7 +873,7 @@ function integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el,A,B,C,D
   broken = plastic_deltaState(ph,en)
   if(broken) return
 
-  broken = integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el)
+  broken = integrateStress(F,subFp0,subFi0,Delta_t,ph,en)
 
 end function integrateStateRK
 
@@ -1006,13 +975,12 @@ end subroutine mechanical_forward
 !--------------------------------------------------------------------------------------------------
 !> @brief calculate stress (P)
 !--------------------------------------------------------------------------------------------------
-module function phase_mechanical_constitutive(Delta_t,co,ip,el) result(converged_)
+module function phase_mechanical_constitutive(Delta_t,co,ce) result(converged_)
 
   real(pReal), intent(in) :: Delta_t
   integer, intent(in) :: &
     co, &
-    ip, &
+    ce
-    el
   logical :: converged_
 
   real(pReal) :: &
@@ -1028,16 +996,15 @@ module function phase_mechanical_constitutive(Delta_t,co,ip,el) result(converged
     subLi0, &
     subF0, &
     subF
-  real(pReal), dimension(:), allocatable :: subState0
+  real(pReal), dimension(plasticState(material_phaseID(co,ce))%sizeState) :: subState0
 
 
-  ph = material_phaseID(co,(el-1)*discretization_nIPs + ip)
+  ph = material_phaseID(co,ce)
-  en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip)
+  en = material_phaseEntry(co,ce)
-  sizeDotState = plasticState(ph)%sizeDotState
 
+  subState0 = plasticState(ph)%state0(:,en)
   subLi0 = phase_mechanical_Li0(ph)%data(1:3,1:3,en)
   subLp0 = phase_mechanical_Lp0(ph)%data(1:3,1:3,en)
-  allocate(subState0,source=plasticState(ph)%State0(:,en))
   subFp0 = phase_mechanical_Fp0(ph)%data(1:3,1:3,en)
   subFi0 = phase_mechanical_Fi0(ph)%data(1:3,1:3,en)
   subF0  = phase_mechanical_F0(ph)%data(1:3,1:3,en)
@@ -1070,7 +1037,7 @@ module function phase_mechanical_constitutive(Delta_t,co,ip,el) result(converged
       subStep       = num%subStepSizeCryst * subStep
       phase_mechanical_Fp(ph)%data(1:3,1:3,en) = subFp0
       phase_mechanical_Fi(ph)%data(1:3,1:3,en) = subFi0
-      phase_mechanical_S(ph)%data(1:3,1:3,en) = phase_mechanical_S0(ph)%data(1:3,1:3,en)            ! why no subS0 ? is S0 of any use?
+      phase_mechanical_S(ph)%data(1:3,1:3,en) = phase_mechanical_S0(ph)%data(1:3,1:3,en)
       if (subStep < 1.0_pReal) then                                                                 ! actual (not initial) cutback
         phase_mechanical_Lp(ph)%data(1:3,1:3,en) = subLp0
         phase_mechanical_Li(ph)%data(1:3,1:3,en) = subLi0
@@ -1082,9 +1049,10 @@ module function phase_mechanical_constitutive(Delta_t,co,ip,el) result(converged
 !--------------------------------------------------------------------------------------------------
 !  prepare for integration
     if (todo) then
+      sizeDotState = plasticState(ph)%sizeDotState
       subF = subF0 &
            + subStep * (phase_mechanical_F(ph)%data(1:3,1:3,en) - phase_mechanical_F0(ph)%data(1:3,1:3,en))
-      converged_ = .not. integrateState(subF0,subF,subFp0,subFi0,subState0(1:sizeDotState),subStep * Delta_t,co,ip,el)
+      converged_ = .not. integrateState(subF0,subF,subFp0,subFi0,subState0(1:sizeDotState),subStep * Delta_t,ph,en)
     endif
 
   enddo cutbackLooping

src/phase_mechanical_plastic.f90

@@ -160,16 +160,14 @@ submodule(phase:mechanical) plastic
         dotState
     end function dislotungsten_dotState
 
-    module subroutine nonlocal_dotState(Mp,timestep,ph,en,ip,el)
+    module subroutine nonlocal_dotState(Mp,timestep,ph,en)
       real(pReal), dimension(3,3), intent(in) :: &
         Mp                                                                                          !< MandelStress
       real(pReal), intent(in) :: &
         timestep                                                                                    !< substepped crystallite time increment
       integer, intent(in) :: &
         ph, &
-        en, &
+        en
-        ip, &                                                                                       !< current integration point
-        el                                                                                          !< current element number
     end subroutine nonlocal_dotState
 
     module subroutine dislotwin_dependentState(ph,en)
@@ -184,12 +182,10 @@ submodule(phase:mechanical) plastic
         en
     end subroutine dislotungsten_dependentState
 
-    module subroutine nonlocal_dependentState(ph, en, ip, el)
+    module subroutine nonlocal_dependentState(ph,en)
       integer, intent(in) :: &
         ph, &
-        en, &
+        en
-        ip, &                                                                                       !< current integration point
-        el                                                                                          !< current element number
     end subroutine nonlocal_dependentState
 
     module subroutine plastic_kinehardening_deltaState(Mp,ph,en)
@@ -299,12 +295,9 @@ end subroutine plastic_LpAndItsTangents
 !--------------------------------------------------------------------------------------------------
 !> @brief contains the constitutive equation for calculating the rate of change of microstructure
 !--------------------------------------------------------------------------------------------------
-module function plastic_dotState(subdt,co,ip,el,ph,en) result(dotState)
+module function plastic_dotState(subdt,ph,en) result(dotState)
 
   integer, intent(in) :: &
-    co, &                                                                                           !< component-ID of integration point
-    ip, &                                                                                           !< integration point
-    el, &                                                                                           !< element
     ph, &
     en
   real(pReal),  intent(in) :: &
@@ -337,7 +330,7 @@ module function plastic_dotState(subdt,co,ip,el,ph,en) result(dotState)
         dotState = dislotungsten_dotState(Mp,ph,en)
 
       case (PLASTIC_NONLOCAL_ID) plasticType
-        call nonlocal_dotState(Mp,subdt,ph,en,ip,el)
+        call nonlocal_dotState(Mp,subdt,ph,en)
         dotState = plasticState(ph)%dotState(:,en)
 
     end select plasticType
@@ -349,21 +342,13 @@ end function plastic_dotState
 !--------------------------------------------------------------------------------------------------
 !> @brief calls microstructure function of the different plasticity constitutive models
 !--------------------------------------------------------------------------------------------------
-module subroutine plastic_dependentState(co, ip, el)
+module subroutine plastic_dependentState(ph,en)
 
   integer, intent(in) :: &
-    co, &                                                                                           !< component-ID of integration point
-    ip, &                                                                                           !< integration point
-    el                                                                                              !< element
-
-  integer :: &
     ph, &
     en
 
 
-  ph = material_phaseID(co,(el-1)*discretization_nIPs + ip)
-  en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip)
-
   plasticType: select case (phase_plasticity(ph))
 
     case (PLASTIC_DISLOTWIN_ID) plasticType
@@ -373,7 +358,7 @@ module subroutine plastic_dependentState(co, ip, el)
       call dislotungsten_dependentState(ph,en)
 
     case (PLASTIC_NONLOCAL_ID) plasticType
-      call nonlocal_dependentState(ph,en,ip,el)
+      call nonlocal_dependentState(ph,en)
 
   end select plasticType
 
@@ -391,10 +376,9 @@ module function plastic_deltaState(ph, en) result(broken)
     en
   logical :: broken
 
-  real(pReal),               dimension(3,3) :: &
+  real(pReal), dimension(3,3) :: &
     Mp
   integer :: &
-    myOffset, &
     mySize
 
 
@@ -419,7 +403,7 @@ module function plastic_deltaState(ph, en) result(broken)
 
       broken = any(IEEE_is_NaN(plasticState(ph)%deltaState(:,en)))
       if (.not. broken) then
-        mySize = plasticState(ph)%sizeDeltaState
+        mySize   = plasticState(ph)%sizeDeltaState
         plasticState(ph)%deltaState2(1:mySize,en) = plasticState(ph)%deltaState2(1:mySize,en) &
                                                   + plasticState(ph)%deltaState(1:mySize,en)
       end if

src/phase_mechanical_plastic_nonlocal.f90

@@ -15,6 +15,9 @@ submodule(phase:plastic) nonlocal
 
   type :: tGeometry
     real(pReal), dimension(:), allocatable :: V_0
+    integer, dimension(:,:,:), allocatable :: IPneighborhood
+    real(pReal), dimension(:,:), allocatable :: IParea, IPcoordinates
+    real(pReal), dimension(:,:,:), allocatable :: IPareaNormal
   end type tGeometry
 
   type(tGeometry), dimension(:), allocatable :: geom
@@ -123,8 +126,10 @@ submodule(phase:plastic) nonlocal
 
   type :: tNonlocalDependentState
     real(pReal), allocatable, dimension(:,:) :: &
-     tau_pass, &
+      tau_pass, &
-     tau_Back
+      tau_Back
+    real(pReal), allocatable, dimension(:,:,:,:,:) :: &
+      compatibility
   end type tNonlocalDependentState
 
   type :: tNonlocalState
@@ -160,7 +165,6 @@ submodule(phase:plastic) nonlocal
     state0
 
   type(tParameters), dimension(:), allocatable :: param                                             !< containers of constitutive parameters
-
   type(tNonlocalDependentState), dimension(:), allocatable :: dependentState
 
 contains
@@ -406,6 +410,10 @@ module function plastic_nonlocal_init() result(myPlasticity)
     call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,sizeDeltaState,0)     ! ToDo: state structure does not follow convention
 
     allocate(geom(ph)%V_0(Nmembers))
+    allocate(geom(ph)%IPneighborhood(3,nIPneighbors,Nmembers))
+    allocate(geom(ph)%IPareaNormal(3,nIPneighbors,Nmembers))
+    allocate(geom(ph)%IParea(nIPneighbors,Nmembers))
+    allocate(geom(ph)%IPcoordinates(3,Nmembers))
     call storeGeometry(ph)
 
     if(plasticState(ph)%nonlocal .and. .not. allocated(IPneighborhood)) &
@@ -489,6 +497,7 @@ module function plastic_nonlocal_init() result(myPlasticity)
 
     allocate(dst%tau_pass(prm%sum_N_sl,Nmembers),source=0.0_pReal)
     allocate(dst%tau_back(prm%sum_N_sl,Nmembers),source=0.0_pReal)
+    allocate(dst%compatibility(2,maxval(param%sum_N_sl),maxval(param%sum_N_sl),nIPneighbors,Nmembers),source=0.0_pReal)
     end associate
 
     if (Nmembers > 0) call stateInit(ini,ph,Nmembers)
@@ -543,13 +552,11 @@ end function plastic_nonlocal_init
 !--------------------------------------------------------------------------------------------------
 !> @brief calculates quantities characterizing the microstructure
 !--------------------------------------------------------------------------------------------------
-module subroutine nonlocal_dependentState(ph, en, ip, el)
+module subroutine nonlocal_dependentState(ph, en)
 
   integer, intent(in) :: &
     ph, &
-    en, &
+    en
-    ip, &
-    el
 
   integer :: &
     no, &                                                                                           !< neighbor offset
@@ -650,12 +657,12 @@ module subroutine nonlocal_dependentState(ph, en, ip, el)
     nRealNeighbors = 0.0_pReal
     neighbor_rhoTotal = 0.0_pReal
     do n = 1,nIPneighbors
-      neighbor_el = IPneighborhood(1,n,ip,el)
+      neighbor_el = geom(ph)%IPneighborhood(1,n,en)
-      neighbor_ip = IPneighborhood(2,n,ip,el)
+      neighbor_ip = geom(ph)%IPneighborhood(2,n,en)
-      no = material_phasememberAt(1,neighbor_ip,neighbor_el)
-      if (neighbor_el > 0 .and. neighbor_ip > 0) then
-        if (material_phaseAt(1,neighbor_el) == ph) then
 
+      if (neighbor_el > 0 .and. neighbor_ip > 0) then
+        if (material_phaseID(1,(neighbor_el-1)*discretization_nIPs + neighbor_ip) == ph) then
+            no = material_phaseEntry(1,(neighbor_el-1)*discretization_nIPs + neighbor_ip)
             nRealNeighbors = nRealNeighbors + 1.0_pReal
             rho_neighbor0 = getRho0(ph,no)
 
@@ -665,11 +672,11 @@ module subroutine nonlocal_dependentState(ph, en, ip, el)
             neighbor_rhoTotal(1,:,n) = sum(abs(rho_neighbor0(:,edg)),2)
             neighbor_rhoTotal(2,:,n) = sum(abs(rho_neighbor0(:,scr)),2)
 
-            connection_latticeConf(1:3,n) = matmul(invFe, discretization_IPcoords(1:3,neighbor_el+neighbor_ip-1) &
+            connection_latticeConf(1:3,n) = matmul(invFe, geom(ph)%IPcoordinates(1:3,no) &
-                                          - discretization_IPcoords(1:3,el+neighbor_ip-1))
+                                          - geom(ph)%IPcoordinates(1:3,en))
-            normal_latticeConf = matmul(transpose(invFp), IPareaNormal(1:3,n,ip,el))
+            normal_latticeConf = matmul(transpose(invFp), geom(ph)%IPareaNormal(1:3,n,en))
             if (math_inner(normal_latticeConf,connection_latticeConf(1:3,n)) < 0.0_pReal) &         ! neighboring connection points in opposite direction to face normal: must be periodic image
-              connection_latticeConf(1:3,n) = normal_latticeConf * IPvolume(ip,el)/IParea(n,ip,el)  ! instead take the surface normal scaled with the diameter of the cell
+              connection_latticeConf(1:3,n) = normal_latticeConf * geom(ph)%V_0(en)/geom(ph)%IParea(n,en)  ! instead take the surface normal scaled with the diameter of the cell
         else
           ! local neighbor or different lattice structure or different constitution instance -> use central values instead
           connection_latticeConf(1:3,n) = 0.0_pReal
@@ -939,7 +946,7 @@ end subroutine plastic_nonlocal_deltaState
 !> @brief calculates the rate of change of microstructure
 !---------------------------------------------------------------------------------------------------
 module subroutine nonlocal_dotState(Mp,timestep, &
-                                    ph,en,ip,el)
+                                    ph,en)
 
   real(pReal), dimension(3,3), intent(in) :: &
     Mp                                                                                              !< MandelStress
@@ -947,9 +954,7 @@ module subroutine nonlocal_dotState(Mp,timestep, &
     timestep                                                                                        !< substepped crystallite time increment
   integer, intent(in) :: &
     ph, &
-    en, &
+    en
-    ip, &                                                                                           !< current integration point
-    el                                                                                              !< current element number
 
   integer ::  &
     c, &                                                                                            !< character of dislocation
@@ -1099,7 +1104,7 @@ module subroutine nonlocal_dotState(Mp,timestep, &
                                          - rhoDip(s,1) / timestep - rhoDotAthermalAnnihilation(s,9) &
                                                                   - rhoDotSingle2DipoleGlide(s,9))  ! make sure that we do not annihilate more dipoles than we have
 
-  rhoDot = rhoDotFlux(timestep, ph,en,ip,el) &
+  rhoDot = rhoDotFlux(timestep, ph,en) &
          + rhoDotMultiplication &
          + rhoDotSingle2DipoleGlide &
          + rhoDotAthermalAnnihilation &
@@ -1110,7 +1115,7 @@ module subroutine nonlocal_dotState(Mp,timestep, &
      .or. any(rho(:,dip) + rhoDot(:,9:10) * timestep < -prm%atol_rho)) then
 #ifdef DEBUG
     if (debugConstitutive%extensive) then
-      print'(a,i5,a,i2)', '<< CONST >> evolution rate leads to negative density at el ',el,' ip ',ip
+      print'(a,i5,a,i2)', '<< CONST >> evolution rate leads to negative density at ph ',ph,' en ',en
       print'(a)', '<< CONST >> enforcing cutback !!!'
     end if
 #endif
@@ -1129,20 +1134,17 @@ end subroutine nonlocal_dotState
 !> @brief calculates the rate of change of microstructure
 !---------------------------------------------------------------------------------------------------
 #if __INTEL_COMPILER >= 2020
-non_recursive function rhoDotFlux(timestep,ph,en,ip,el)
+non_recursive function rhoDotFlux(timestep,ph,en)
 #else
-function rhoDotFlux(timestep,ph,en,ip,el)
+function rhoDotFlux(timestep,ph,en)
 #endif
   real(pReal), intent(in) :: &
     timestep                                                                                        !< substepped crystallite time increment
   integer, intent(in) :: &
     ph, &
-    en, &
+    en
-    ip, &                                                                                           !< current integration point
-    el                                                                                              !< current element number
 
   integer ::  &
-    neighbor_ph, &                                                                                  !< phase of my neighbor's plasticity
     ns, &                                                                                           !< short notation for the total number of active slip systems
     c, &                                                                                            !< character of dislocation
     n, &                                                                                            !< index of my current neighbor
@@ -1215,14 +1217,14 @@ function rhoDotFlux(timestep,ph,en,ip,el)
     !*** check CFL (Courant-Friedrichs-Lewy) condition for flux
     if (any( abs(dot_gamma) > 0.0_pReal &                                                           ! any active slip system ...
             .and. prm%C_CFL * abs(v0) * timestep &
-                > IPvolume(ip,el) / maxval(IParea(:,ip,el)))) then                                  ! ...with velocity above critical value (we use the reference volume and area for simplicity here)
+                > geom(ph)%V_0(en)/ maxval(geom(ph)%IParea(:,en)))) then                            ! ...with velocity above critical value (we use the reference volume and area for simplicity here)
 #ifdef DEBUG
     if (debugConstitutive%extensive) then
-      print'(a,i5,a,i2)', '<< CONST >> CFL condition not fullfilled at el ',el,' ip ',ip
+      print'(a,i5,a,i2)', '<< CONST >> CFL condition not fullfilled at ph ',ph,' en ',en
       print'(a,e10.3,a,e10.3)', '<< CONST >> velocity is at  ', &
         maxval(abs(v0), abs(dot_gamma) > 0.0_pReal &
                        .and.  prm%C_CFL * abs(v0) * timestep &
-                             > IPvolume(ip,el) / maxval(IParea(:,ip,el))), &
+                             > geom(ph)%V_0(en) / maxval(geom(ph)%IParea(:,en))), &
         ' at a timestep of ',timestep
       print*, '<< CONST >> enforcing cutback !!!'
     end if
@@ -1245,19 +1247,18 @@ function rhoDotFlux(timestep,ph,en,ip,el)
 
     neighbors: do n = 1,nIPneighbors
 
-      neighbor_el = IPneighborhood(1,n,ip,el)
+      neighbor_el = geom(ph)%IPneighborhood(1,n,en)
-      neighbor_ip = IPneighborhood(2,n,ip,el)
+      neighbor_ip = geom(ph)%IPneighborhood(2,n,en)
-      neighbor_n  = IPneighborhood(3,n,ip,el)
+      neighbor_n  = geom(ph)%IPneighborhood(3,n,en)
-      np = material_phaseAt(1,neighbor_el)
+      np = material_phaseID(1,(neighbor_el-1)*discretization_nIPs + neighbor_ip)
-      no = material_phasememberAt(1,neighbor_ip,neighbor_el)
+      no = material_phaseEntry(1,(neighbor_el-1)*discretization_nIPs + neighbor_ip)
 
       opposite_neighbor = n + mod(n,2) - mod(n+1,2)
-      opposite_el = IPneighborhood(1,opposite_neighbor,ip,el)
+      opposite_el = geom(ph)%IPneighborhood(1,opposite_neighbor,en)
-      opposite_ip = IPneighborhood(2,opposite_neighbor,ip,el)
+      opposite_ip = geom(ph)%IPneighborhood(2,opposite_neighbor,en)
-      opposite_n  = IPneighborhood(3,opposite_neighbor,ip,el)
+      opposite_n  = geom(ph)%IPneighborhood(3,opposite_neighbor,en)
 
       if (neighbor_n > 0) then                                                                      ! if neighbor exists, average deformation gradient
-        neighbor_ph = material_phaseAt(1,neighbor_el)
         neighbor_F = phase_mechanical_F(np)%data(1:3,1:3,no)
         neighbor_Fe = matmul(neighbor_F, math_inv33(phase_mechanical_Fp(np)%data(1:3,1:3,no)))
         Favg = 0.5_pReal * (my_F + neighbor_F)
@@ -1276,12 +1277,12 @@ function rhoDotFlux(timestep,ph,en,ip,el)
       !* The entering flux from my neighbor will be distributed on my slip systems according to the
       !* compatibility
       if (neighbor_n > 0) then
-      if (phase_plasticity(material_phaseAt(1,neighbor_el)) == PLASTIC_NONLOCAL_ID .and. &
+      if (phase_plasticity(np) == PLASTIC_NONLOCAL_ID .and. &
-          any(compatibility(:,:,:,n,ip,el) > 0.0_pReal)) then
+          any(dependentState(ph)%compatibility(:,:,:,n,en) > 0.0_pReal)) then
 
         forall (s = 1:ns, t = 1:4)
-          neighbor_v0(s,t) =          plasticState(np)%state0(iV   (s,t,neighbor_ph),no)
+          neighbor_v0(s,t) =          plasticState(np)%state0(iV   (s,t,np),no)
-          neighbor_rhoSgl0(s,t) = max(plasticState(np)%state0(iRhoU(s,t,neighbor_ph),no),0.0_pReal)
+          neighbor_rhoSgl0(s,t) = max(plasticState(np)%state0(iRhoU(s,t,np),no),0.0_pReal)
         endforall
 
         where (neighbor_rhoSgl0 * IPvolume(neighbor_ip,neighbor_el) ** 0.667_pReal < prm%rho_min &
@@ -1301,12 +1302,12 @@ function rhoDotFlux(timestep,ph,en,ip,el)
                 .and. v0(s,t) * neighbor_v0(s,t) >= 0.0_pReal ) then                                ! ... only if no sign change in flux density
               lineLength = neighbor_rhoSgl0(s,t) * neighbor_v0(s,t) &
                          * math_inner(m(1:3,s,t), normal_neighbor2me) * area                        ! positive line length that wants to enter through this interface
-              where (compatibility(c,:,s,n,ip,el) > 0.0_pReal) &
+              where (dependentState(ph)%compatibility(c,:,s,n,en) > 0.0_pReal) &
                 rhoDotFlux(:,t) = rhoDotFlux(1:ns,t) &
-                                + lineLength/IPvolume(ip,el)*compatibility(c,:,s,n,ip,el)**2        ! transferring to equally signed mobile dislocation type
+                                + lineLength/geom(ph)%V_0(en)*dependentState(ph)%compatibility(c,:,s,n,en)**2        ! transferring to equally signed mobile dislocation type
-              where (compatibility(c,:,s,n,ip,el) < 0.0_pReal) &
+              where (dependentState(ph)%compatibility(c,:,s,n,en) < 0.0_pReal) &
                 rhoDotFlux(:,topp) = rhoDotFlux(:,topp) &
-                                   + lineLength/IPvolume(ip,el)*compatibility(c,:,s,n,ip,el)**2     ! transferring to opposite signed mobile dislocation type
+                                   + lineLength/geom(ph)%V_0(en)*dependentState(ph)%compatibility(c,:,s,n,en)**2     ! transferring to opposite signed mobile dislocation type
 
             end if
           end do
@@ -1322,28 +1323,28 @@ function rhoDotFlux(timestep,ph,en,ip,el)
       !* In case of reduced transmissivity, part of the leaving flux is stored as dead dislocation density.
       !* That means for an interface of zero transmissivity the leaving flux is fully converted to dead dislocations.
       if (opposite_n > 0) then
-      if (phase_plasticity(material_phaseAt(1,opposite_el)) == PLASTIC_NONLOCAL_ID) then
+      if (phase_plasticity(np) == PLASTIC_NONLOCAL_ID) then
 
         normal_me2neighbor_defConf = math_det33(Favg) &
-                                   * matmul(math_inv33(transpose(Favg)),IPareaNormal(1:3,n,ip,el))  ! normal of the interface in (average) deformed configuration (pointing en => neighbor)
+                                   * matmul(math_inv33(transpose(Favg)),geom(ph)%IPareaNormal(1:3,n,en))  ! normal of the interface in (average) deformed configuration (pointing en => neighbor)
         normal_me2neighbor = matmul(transpose(my_Fe), normal_me2neighbor_defConf) &
                            / math_det33(my_Fe)                                                      ! interface normal in my lattice configuration
-        area = IParea(n,ip,el) * norm2(normal_me2neighbor)
+        area = geom(ph)%IParea(n,en) * norm2(normal_me2neighbor)
         normal_me2neighbor = normal_me2neighbor / norm2(normal_me2neighbor)                         ! normalize the surface normal to unit length
         do s = 1,ns
           do t = 1,4
             c = (t + 1) / 2
             if (v0(s,t) * math_inner(m(1:3,s,t), normal_me2neighbor) > 0.0_pReal ) then             ! flux from en to my neighbor == leaving flux for en (might also be a pure flux from my mobile density to dead density if interface not at all transmissive)
               if (v0(s,t) * neighbor_v0(s,t) >= 0.0_pReal) then                                     ! no sign change in flux density
-                transmissivity = sum(compatibility(c,:,s,n,ip,el)**2)                               ! overall transmissivity from this slip system to my neighbor
+                transmissivity = sum(dependentState(ph)%compatibility(c,:,s,n,en)**2)               ! overall transmissivity from this slip system to my neighbor
               else                                                                                  ! sign change in flux density means sign change in stress which does not allow for dislocations to arive at the neighbor
                 transmissivity = 0.0_pReal
               end if
               lineLength = my_rhoSgl0(s,t) * v0(s,t) &
                          * math_inner(m(1:3,s,t), normal_me2neighbor) * area                        ! positive line length of mobiles that wants to leave through this interface
-              rhoDotFlux(s,t) = rhoDotFlux(s,t) - lineLength / IPvolume(ip,el)                      ! subtract dislocation flux from current type
+              rhoDotFlux(s,t) = rhoDotFlux(s,t) - lineLength / geom(ph)%V_0(en)                     ! subtract dislocation flux from current type
               rhoDotFlux(s,t+4) = rhoDotFlux(s,t+4) &
-                                + lineLength / IPvolume(ip,el) * (1.0_pReal - transmissivity) &
+                                + lineLength / geom(ph)%V_0(en) * (1.0_pReal - transmissivity) &
                                 * sign(1.0_pReal, v0(s,t))                                          ! dislocation flux that is not able to leave through interface (because of low transmissivity) will remain as immobile single density at the material point
             end if
           end do
@@ -1364,14 +1365,14 @@ end function rhoDotFlux
 ! plane normals and signed cosine of the angle between the slip directions. Only the largest values
 ! that sum up to a total of 1 are considered, all others are set to zero.
 !--------------------------------------------------------------------------------------------------
-module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e)
+module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,ip,el)
 
   type(tRotationContainer), dimension(:), intent(in) :: &
     orientation                                                                                     ! crystal orientation
   integer, intent(in) :: &
     ph, &
-    i, &
+    ip, &
-    e
+    el
 
   integer :: &
     n, &                                                                                            ! neighbor index
@@ -1397,16 +1398,16 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e)
   associate(prm => param(ph))
     ns = prm%sum_N_sl
 
-    en = material_phaseMemberAt(1,i,e)
+    en = material_phaseEntry(1,(el-1)*discretization_nIPs + ip)
     !*** start out fully compatible
     my_compatibility = 0.0_pReal
     forall(s1 = 1:ns) my_compatibility(:,s1,s1,:) = 1.0_pReal
 
     neighbors: do n = 1,nIPneighbors
-      neighbor_e = IPneighborhood(1,n,i,e)
+      neighbor_e = IPneighborhood(1,n,ip,el)
-      neighbor_i = IPneighborhood(2,n,i,e)
+      neighbor_i = IPneighborhood(2,n,ip,el)
-      neighbor_me = material_phaseMemberAt(1,neighbor_i,neighbor_e)
+      neighbor_me = material_phaseEntry(1,(neighbor_e-1)*discretization_nIPs + neighbor_i)
-      neighbor_phase = material_phaseAt(1,neighbor_e)
+      neighbor_phase = material_phaseID(1,(neighbor_e-1)*discretization_nIPs + neighbor_i)
 
       if (neighbor_e <= 0 .or. neighbor_i <= 0) then
         !* FREE SURFACE
@@ -1465,7 +1466,7 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e)
 
     end do neighbors
 
-    compatibility(:,:,:,:,i,e) = my_compatibility
+    dependentState(ph)%compatibility(:,:,:,:,material_phaseEntry(1,(el-1)*discretization_nIPs + ip)) = my_compatibility
 
   end associate
 
@@ -1756,14 +1757,29 @@ subroutine storeGeometry(ph)
 
   integer, intent(in) :: ph
 
-  integer :: ce, co
+  integer :: ce, co, nCell
   real(pReal), dimension(:), allocatable :: V
+  integer, dimension(:,:,:), allocatable :: neighborhood
+  real(pReal), dimension(:,:), allocatable :: area, coords
+  real(pReal), dimension(:,:,:), allocatable :: areaNormal
 
+  nCell = product(shape(IPvolume))
+
+  V = reshape(IPvolume,[nCell])
+  neighborhood = reshape(IPneighborhood,[3,nIPneighbors,nCell])
+  area = reshape(IParea,[nIPneighbors,nCell])
+  areaNormal = reshape(IPareaNormal,[3,nIPneighbors,nCell])
+  coords = reshape(discretization_IPcoords,[3,nCell])
 
-  V = reshape(IPvolume,[product(shape(IPvolume))])
   do ce = 1, size(material_homogenizationEntry,1)
     do co = 1, homogenization_maxNconstituents
-      if (material_phaseID(co,ce) == ph) geom(ph)%V_0(material_phaseEntry(co,ce)) = V(ce)
+      if (material_phaseID(co,ce) == ph) then
+        geom(ph)%V_0(material_phaseEntry(co,ce)) = V(ce)
+        geom(ph)%IPneighborhood(:,:,material_phaseEntry(co,ce)) = neighborhood(:,:,ce)
+        geom(ph)%IParea(:,material_phaseEntry(co,ce)) = area(:,ce)
+        geom(ph)%IPareaNormal(:,:,material_phaseEntry(co,ce)) = areaNormal(:,:,ce)
+        geom(ph)%IPcoordinates(:,material_phaseEntry(co,ce)) = coords(:,ce)
+      end if
     end do
   end do
 

src/phase_thermal.f90

@@ -96,17 +96,24 @@ module subroutine thermal_init(phases)
 
   do ph = 1, phases%length
     Nmembers = count(material_phaseID == ph)
-    allocate(current(ph)%T(Nmembers),source=300.0_pReal)
+    allocate(current(ph)%T(Nmembers),source=T_ROOM)
     allocate(current(ph)%dot_T(Nmembers),source=0.0_pReal)
     phase => phases%get(ph)
     thermal => phase%get('thermal',defaultVal=emptyDict)
-    param(ph)%C_p = thermal%get_asFloat('C_p',defaultVal=0.0_pReal)                                 ! ToDo: make mandatory?
-    param(ph)%K(1,1) = thermal%get_asFloat('K_11',defaultVal=0.0_pReal)                             ! ToDo: make mandatory?
-    param(ph)%K(3,3) = thermal%get_asFloat('K_33',defaultVal=0.0_pReal)                             ! ToDo: depends on symmtery
-    param(ph)%K = lattice_symmetrize_33(param(ph)%K,phase_lattice(ph))
 
-    sources => thermal%get('source',defaultVal=emptyList)
+    ! ToDo: temperature dependency of K and C_p
-    thermal_Nsources(ph) = sources%length
+    if (thermal%length > 0) then
+      param(ph)%C_p = thermal%get_asFloat('C_p')
+      param(ph)%K(1,1) = thermal%get_asFloat('K_11')
+      if (any(phase_lattice(ph) == ['hP','tI'])) param(ph)%K(3,3) = thermal%get_asFloat('K_33')
+      param(ph)%K = lattice_symmetrize_33(param(ph)%K,phase_lattice(ph))
+
+      sources => thermal%get('source',defaultVal=emptyList)
+      thermal_Nsources(ph) = sources%length
+    else
+      thermal_Nsources(ph) = 0
+    end if
+
     allocate(thermalstate(ph)%p(thermal_Nsources(ph)))
 
   enddo