Merge remote-tracking branch 'origin/development' into 56-parallel-hdf5

.gitlab-ci.yml

@@ -79,6 +79,7 @@ variables:
   Abaqus:                          "$Abaqus2017"
   MSC:                             "$MSC2018_1"
   IntelMarc:                       "$IntelCompiler17_0"
+  IntelAbaqus:                     "$IntelCompiler16_4"
   # ++++++++++++ Documentation +++++++++++++++++++++++++++++++++++++++++
   Doxygen1_8_13:                   "Documentation/Doxygen/1.8.13"
   # ------------ Defaults ----------------------------------------------
@@ -150,6 +151,13 @@ Spectral_geometryPacking:
     - release
 
 ###################################################################################################
+Post_AverageDown:
+  stage: postprocessing
+  script: averageDown/test.py
+  except:
+    - master
+    - release
+
 Post_General:
   stage: postprocessing
   script: PostProcessing/test.py
@@ -192,6 +200,13 @@ Post_ParaviewRelated:
     - master
     - release
 
+Post_OrientationConversion:
+  stage: postprocessing
+  script: OrientationConversion/test.py
+  except:
+    - master
+    - release
+
 ###################################################################################################
 Compile_Spectral_Intel:
   stage: compilePETScIntel
@@ -242,6 +257,13 @@ Compile_Intel_Prepare:
     - release
 
 ###################################################################################################
+Thermal:
+  stage: spectral
+  script: Thermal/test.py
+  except:
+    - master
+    - release
+
 Spectral_PackedGeometry:
   stage: spectral
   script: Spectral_PackedGeometry/test.py
@@ -351,13 +373,6 @@ Phenopowerlaw_singleSlip:
     - master
     - release
 
-HybridIA:
-  stage: spectral
-  script: HybridIA/test.py
-  except:
-    - master
-    - release
-
 TextureComponents:
   stage: spectral
   script: TextureComponents/test.py
@@ -431,14 +446,6 @@ SpectralExample:
   only:
     - development
 
-AbaqusExample:
-  stage: example
-  script:
-    - module load $IntelCompiler16_4 $Abaqus
-    - Abaqus_example/test.py
-  only:
-    - development
-
 ###################################################################################################
 SpectralRuntime:
   stage: performance

CMakeLists.txt

@@ -47,7 +47,7 @@ linker:
 
 # CMake will execute each target in the ${petsc_config_makefile}
 # to acquire corresponding PETSc Variables.
-find_program (MAKE_EXECUTABLE NAMES make gmake)
+find_program (MAKE_EXECUTABLE NAMES gmake make)
 # Find the PETSc includes directory settings
 execute_process (COMMAND ${MAKE_EXECUTABLE} --no-print-directory -f ${petsc_config_makefile} "includes"
                  RESULT_VARIABLE PETSC_INCLUDES_RETURN
@@ -280,7 +280,7 @@ if (CMAKE_Fortran_COMPILER_ID STREQUAL "Intel")
   #   ... for uninitialized variables.
   set (DEBUG_FLAGS "${DEBUG_FLAGS} -ftrapuv")
   #   ... initializes stack local variables to an unusual value to aid error detection
-  set (DEBUG_FLAGS "${DEBUG_FLAGS} -fpe-all0")
+  set (DEBUG_FLAGS "${DEBUG_FLAGS} -fpe-all=0")
   #   ... capture all floating-point exceptions, sets -ftz automatically
 
   set (DEBUG_FLAGS "${DEBUG_FLAGS} -warn")

CONFIG

@@ -8,6 +8,6 @@ set DAMASK_NUM_THREADS    = 4
 set MSC_ROOT              = /opt/msc
 set MARC_VERSION          = 2018.1
 
-set ABAQUS_VERSION        = 2018.1
+set ABAQUS_VERSION        = 2017
 
 set DAMASK_HDF5           = OFF

PRIVATE

@@ -1 +1 @@
-Subproject commit d3bc62220544da0a3198c521e0f73fa07898d357
+Subproject commit e9f93abaecafbfbf11072ae70bca213a7201ed38

VERSION

@@ -1 +1 @@
-v2.0.2-889-g47650e94
+v2.0.2-1107-g17716b4f

env/DAMASK.csh

@@ -44,7 +44,6 @@ if ( $?prompt ) then
   echo "DAMASK             $DAMASK_ROOT"
   echo "Spectral Solver    $SOLVER" 
   echo "Post Processing    $PROCESSING"
-  echo "Multithreading     DAMASK_NUM_THREADS=$DAMASK_NUM_THREADS"
   if ( $?PETSC_DIR) then
     echo "PETSc location     $PETSC_DIR"
   endif
@@ -52,8 +51,10 @@ if ( $?prompt ) then
     echo "MSC.Marc/Mentat    $MSC_ROOT"
   endif
   echo
+  echo "Multithreading     DAMASK_NUM_THREADS=$DAMASK_NUM_THREADS"
   echo `limit datasize`
   echo `limit stacksize`
+  echo
 endif
 
 setenv DAMASK_NUM_THREADS $DAMASK_NUM_THREADS

env/DAMASK.sh

@@ -88,6 +88,7 @@ if [ ! -z "$PS1" ]; then
            size=$(( 1024*$(ulimit -s) )); \
            print('{:.4g} {}'.format(size / (1 << ((int(math.log(size,2) / 10) if size else 0) * 10)), \
            ['bytes','KiB','MiB','GiB','TiB','EiB','ZiB'][int(math.log(size,2) / 10) if size else 0]))")
+  echo
 fi
 
 export DAMASK_NUM_THREADS

env/DAMASK.zsh

@@ -81,6 +81,7 @@ if [ ! -z "$PS1" ]; then
            size=$(( 1024*$(ulimit -s) )); \
            print('{:.4g} {}'.format(size / (1 << ((int(math.log(size,2) / 10) if size else 0) * 10)), \
            ['bytes','KiB','MiB','GiB','TiB','EiB','ZiB'][int(math.log(size,2) / 10) if size else 0]))")
+  echo
 fi
 
 export DAMASK_NUM_THREADS

examples/AbaqusStandard/SX_PX_compression.jnl

@@ -1,23 +0,0 @@
-# Save by m.diehl on 2017_12_06-18.38.26; build 2017 2016_09_27-23.54.59 126836
-from abaqus import *
-upgradeMdb(
-    '/nethome/storage/raid4/m.diehl/DAMASK/examples/AbaqusStandard/SX_PX_compression-6.9-1.cae'
-    , 
-    '/nethome/storage/raid4/m.diehl/DAMASK/examples/AbaqusStandard/SX_PX_compression.cae')
-# Save by m.diehl on 2017_12_06-18.38.26; build 2017 2016_09_27-23.54.59 126836
-from part import *
-from material import *
-from section import *
-from assembly import *
-from step import *
-from interaction import *
-from load import *
-from mesh import *
-from optimization import *
-from job import *
-from sketch import *
-from visualization import *
-from connectorBehavior import *
-mdb.jobs['Job_sx-px'].setValues(description='compression', userSubroutine=
-    '$HOME/DAMASK/src/DAMASK_abaqus_std.f')
-# Save by m.diehl on 2017_12_06-18.39.44; build 2017 2016_09_27-23.54.59 126836

examples/AbaqusStandard/material.config

@@ -1,93 +0,0 @@
-#-------------------#
-<homogenization>
-#-------------------#
-
-[dummy]
-mech	none
-
-[poly]
-mech	isostrain
-Nconstituents	10
-
-
-
-#-------------------#
-<microstructure>
-#-------------------#
-
-
-[Aluminum_001]
-crystallite	1
-(constituent)	phase 1	texture	1	fraction	1.0
-
-[Aluminum_10]
-crystallite	1
-(constituent)	phase 1	texture	2	fraction	1.0
-
-#-------------------#
-<crystallite>
-#-------------------#
-
-[orientation]
-(output) eulerangles
-(output) grainrotation  # deviation from initial orientation as axis (1-3) and angle in degree (4)
-
-#-------------------#
-<phase>
-#-------------------#
-
-
-[Aluminum_phenopowerlaw]
-# slip only
-elasticity              hooke
-plasticity              phenopowerlaw
-
-(output)                resistance_slip
-(output)                shearrate_slip
-(output)                resolvedstress_slip
-(output)                totalshear
-(output)                resistance_twin
-(output)                shearrate_twin
-(output)                resolvedstress_twin
-(output)                totalvolfrac
-
-lattice_structure       fcc
-Nslip                   12
-Ntwin                    0
-
-c11                     106.75e9
-c12                     60.41e9
-c44                     28.34e9
-
-gdot0_slip              0.001
-n_slip                  20
-tau0_slip                 31e6                # per family
-tausat_slip               63e6               # per family
-a_slip                 2.25
-gdot0_twin              0.001
-n_twin                  20
-tau0_twin                 31e6               # per family
-s_pr                    0                        # push-up factor for slip saturation due to twinning
-twin_b                  0
-twin_c                  0
-twin_d                  0
-twin_e                  0
-h0_slipslip             75e6
-h0_sliptwin             0
-h0_twinslip             0
-h0_twintwin             0
-interaction_slipslip    1 1 1.4 1.4 1.4 1.4
-interaction_sliptwin    1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
-interaction_twinslip    1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
-interaction_twintwin    1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 
-atol_resistance         1
-
-
-#-------------------#
-<texture>
-#-------------------#
-
-[001]
-(gauss)		phi1	0.000	Phi	45.000	phi2	0.000	scatter	0.000	fraction	1.000
-
-[random]

examples/AbaqusStandard/numerics.config

@@ -1 +0,0 @@
-fixed_seed 1697667030

lib/damask/orientation.py

@@ -27,15 +27,22 @@ class Rodrigues:
 
 # ******************************************************************************************
 class Quaternion:
-    """
+    u"""
     Orientation represented as unit quaternion.
     
-    All methods and naming conventions based on http://www.euclideanspace.com/maths/algebra/realNormedAlgebra/quaternions.
+    All methods and naming conventions based on Rowenhorst_etal2015
+    Convention 1: coordinate frames are right-handed
+    Convention 2: a rotation angle ω is taken to be positive for a counterclockwise rotation
+                  when viewing from the end point of the rotation axis unit vector towards the origin
+    Convention 3: rotations will be interpreted in the passive sense
+    Convention 4: Euler angle triplets are implemented using the Bunge convention,
+                  with the angular ranges as [0, 2π],[0, π],[0, 2π]
+    Convention 5: the rotation angle ω is limited to the interval [0, π]
 
     w is the real part, (x, y, z) are the imaginary parts.
-    Representation of rotation is in ACTIVE form!
+
-    (Derived directly or through angleAxis, Euler angles, or active matrix)
+    Vector "a" (defined in coordinate system "A") is passively rotated
-    Vector "a" (defined in coordinate system "A") is actively rotated to new coordinates "b".
+               resulting in new coordinates "b" when expressed in system "B".
     b = Q * a
     b = np.dot(Q.asMatrix(),a)
     """
@@ -309,10 +316,12 @@ class Quaternion:
       return np.outer([i for i in self],[i for i in self])
 
     def asMatrix(self):
-      return np.array(
+      qbarhalf = 0.5*(self.w**2 - self.x**2 - self.y**2 - self.z**2)
-        [[1.0-2.0*(self.y*self.y+self.z*self.z),     2.0*(self.x*self.y-self.z*self.w),     2.0*(self.x*self.z+self.y*self.w)],
+      return 2.0*np.array(
-         [    2.0*(self.x*self.y+self.z*self.w), 1.0-2.0*(self.x*self.x+self.z*self.z),     2.0*(self.y*self.z-self.x*self.w)],
+        [[      qbarhalf + self.x**2    , self.x*self.y - self.w*self.z, self.x*self.z + self.w*self.y],
-         [    2.0*(self.x*self.z-self.y*self.w),     2.0*(self.x*self.w+self.y*self.z), 1.0-2.0*(self.x*self.x+self.y*self.y)]])
+         [ self.x*self.y + self.w*self.z,      qbarhalf + self.y**2    , self.y*self.z - self.w*self.x],
+         [ self.x*self.z - self.w*self.y, self.y*self.z + self.w*self.x,      qbarhalf + self.z**2    ],
+        ])
 
     def asAngleAxis(self,
                     degrees = False):
@@ -335,52 +344,23 @@ class Quaternion:
       return np.inf*np.ones(3) if self.w == 0.0 else np.array([self.x, self.y, self.z])/self.w
 
     def asEulers(self,
-                 type = "bunge",
+                 degrees = False):
-                 degrees = False,
+      """Orientation as Bunge-Euler angles."""
-                 standardRange = False):
+      q03 = self.w**2+self.z**2
-      """
+      q12 = self.x**2+self.y**2
-      Orientation as Bunge-Euler angles.
+      chi = np.sqrt(q03*q12)
+      
+      if abs(chi) < 1e-10 and abs(q12) < 1e-10:
+        eulers = np.array([math.atan2(-2*self.w*self.z,self.w**2-self.z**2),0,0])
+      elif abs(chi) < 1e-10 and abs(q03) < 1e-10:
+        eulers = np.array([math.atan2( 2*self.x*self.y,self.x**2-self.y**2),np.pi,0])
+      else:
+        eulers = np.array([math.atan2((self.x*self.z-self.w*self.y)/chi,(-self.w*self.x-self.y*self.z)/chi),
+                           math.atan2(2*chi,q03-q12),
+                           math.atan2((self.w*self.y+self.x*self.z)/chi,( self.y*self.z-self.w*self.x)/chi),
+                          ])
 
-      Conversion of ACTIVE rotation to Euler angles taken from:
+      return np.degrees(eulers) if degrees else eulers
-      Melcher, A.; Unser, A.; Reichhardt, M.; Nestler, B.; Poetschke, M.; Selzer, M.
-      Conversion of EBSD data by a quaternion based algorithm to be used for grain structure simulations
-      Technische Mechanik 30 (2010) pp 401--413.
-      """
-      angles = [0.0,0.0,0.0]
-
-      if type.lower() == 'bunge' or type.lower() == 'zxz':
-        if   abs(self.x) < 1e-4 and abs(self.y) < 1e-4:
-          x = self.w**2 - self.z**2
-          y = 2.*self.w*self.z
-          angles[0] = math.atan2(y,x)
-        elif abs(self.w) < 1e-4 and abs(self.z) < 1e-4:
-          x = self.x**2 - self.y**2
-          y = 2.*self.x*self.y
-          angles[0] = math.atan2(y,x)
-          angles[1] = math.pi
-        else:
-          chi = math.sqrt((self.w**2 + self.z**2)*(self.x**2 + self.y**2))
-
-          x = (self.w * self.x - self.y * self.z)/2./chi
-          y = (self.w * self.y + self.x * self.z)/2./chi
-          angles[0] = math.atan2(y,x)
-
-          x = self.w**2 + self.z**2 - (self.x**2 + self.y**2)
-          y = 2.*chi
-          angles[1] = math.atan2(y,x)
-
-          x = (self.w * self.x + self.y * self.z)/2./chi
-          y = (self.z * self.x - self.y * self.w)/2./chi
-          angles[2] = math.atan2(y,x)
-
-        if standardRange:
-          angles[0] %= 2*math.pi
-          if angles[1] < 0.0:
-            angles[1] += math.pi
-            angles[2] *= -1.0
-          angles[2] %= 2*math.pi
-
-      return np.degrees(angles) if degrees else angles
 
 
 #    # Static constructors
@@ -408,7 +388,7 @@ class Quaternion:
       halfangle = math.atan(np.linalg.norm(rodrigues))
       c = math.cos(halfangle)
       w = c
-      x,y,z = c*rodrigues
+      x,y,z = rodrigues/c
       return cls([w,x,y,z])
 
 
@@ -431,24 +411,19 @@ class Quaternion:
     @classmethod
     def fromEulers(cls,
                    eulers,
-                   type = 'Bunge',
                    degrees = False):
       if not isinstance(eulers, np.ndarray): eulers = np.array(eulers,dtype='d')
       eulers = np.radians(eulers) if degrees else eulers
 
-      c = np.cos(0.5 * eulers)
+      sigma = 0.5*(eulers[0]+eulers[2])
-      s = np.sin(0.5 * eulers)
+      delta = 0.5*(eulers[0]-eulers[2])
+      c = np.cos(0.5*eulers[1])
+      s = np.sin(0.5*eulers[1])
 
-      if type.lower() == 'bunge' or type.lower() == 'zxz':
+      w =  c * np.cos(sigma) 
-        w =   c[0] * c[1] * c[2] - s[0] * c[1] * s[2]
+      x = -s * np.cos(delta) 
-        x =   c[0] * s[1] * c[2] + s[0] * s[1] * s[2]
+      y = -s * np.sin(delta)
-        y = - c[0] * s[1] * s[2] + s[0] * s[1] * c[2]
+      z = -c * np.sin(sigma) 
-        z =   c[0] * c[1] * s[2] + s[0] * c[1] * c[2]
-      else:
-        w = c[0] * c[1] * c[2] - s[0] * s[1] * s[2]
-        x = s[0] * s[1] * c[2] + c[0] * c[1] * s[2]
-        y = s[0] * c[1] * c[2] + c[0] * s[1] * s[2]
-        z = c[0] * s[1] * c[2] - s[0] * c[1] * s[2]
       return cls([w,x,y,z])
 
 
@@ -460,49 +435,16 @@ class Quaternion:
       if m.shape != (3,3) and np.prod(m.shape) == 9:
         m = m.reshape(3,3)
 
-      tr = np.trace(m)
+      w = 0.5*math.sqrt(1.+m[0,0]+m[1,1]+m[2,2])
-      if tr > 1e-8:
+      x = 0.5*math.sqrt(1.+m[0,0]-m[1,1]-m[2,2])
-        s = math.sqrt(tr + 1.0)*2.0
+      y = 0.5*math.sqrt(1.-m[0,0]+m[1,1]-m[2,2])
+      z = 0.5*math.sqrt(1.-m[0,0]-m[1,1]+m[2,2])
 
-        return cls(
+      x *= -1 if m[2,1] < m[1,2] else 1
-          [ s*0.25,
+      y *= -1 if m[0,2] < m[2,0] else 1
-            (m[2,1] - m[1,2])/s,
+      z *= -1 if m[1,0] < m[0,1] else 1
-            (m[0,2] - m[2,0])/s,
-            (m[1,0] - m[0,1])/s,
-          ])
 
-      elif m[0,0] > m[1,1] and m[0,0] > m[2,2]:
+      return cls( np.array([w,x,y,z])/math.sqrt(w**2 + x**2 + y**2 + z**2))
-        t = m[0,0] - m[1,1] - m[2,2] + 1.0
-        s = 2.0*math.sqrt(t)
-
-        return cls(
-          [ (m[2,1] - m[1,2])/s,
-            s*0.25,
-            (m[0,1] + m[1,0])/s,
-            (m[2,0] + m[0,2])/s,
-          ])
-
-      elif m[1,1] > m[2,2]:
-        t = -m[0,0] + m[1,1] - m[2,2] + 1.0
-        s = 2.0*math.sqrt(t)
-
-        return cls(
-          [ (m[0,2] - m[2,0])/s,
-            (m[0,1] + m[1,0])/s,
-            s*0.25,
-            (m[1,2] + m[2,1])/s,
-          ])
-
-      else:
-        t = -m[0,0] - m[1,1] + m[2,2] + 1.0
-        s = 2.0*math.sqrt(t)
-
-        return cls(
-          [ (m[1,0] - m[0,1])/s,
-            (m[2,0] + m[0,2])/s,
-            (m[1,2] + m[2,1])/s,
-            s*0.25,
-          ])
 
 
     @classmethod
@@ -663,7 +605,7 @@ class Symmetry:
                             quaternion,
                             who = []):
     """List of symmetrically equivalent quaternions based on own symmetry."""
-    return [quaternion*q for q in self.symmetryQuats(who)]
+    return [q*quaternion for q in self.symmetryQuats(who)]
 
 
   def inFZ(self,R):
@@ -829,7 +771,7 @@ class Orientation:
       else:
         self.quaternion = Quaternion.fromRandom(randomSeed=random)
     elif isinstance(Eulers, np.ndarray) and Eulers.shape == (3,):                                   # based on given Euler angles
-      self.quaternion = Quaternion.fromEulers(Eulers,type='bunge',degrees=degrees)
+      self.quaternion = Quaternion.fromEulers(Eulers,degrees=degrees)
     elif isinstance(matrix, np.ndarray) :                                                           # based on given rotation matrix
       self.quaternion = Quaternion.fromMatrix(matrix)
     elif isinstance(angleAxis, np.ndarray) and angleAxis.shape == (4,):                             # based on given angle and rotation axis
@@ -855,16 +797,15 @@ class Orientation:
     return 'Symmetry: %s\n' % (self.symmetry) + \
            'Quaternion: %s\n' % (self.quaternion) + \
            'Matrix:\n%s\n' % ( '\n'.join(['\t'.join(map(str,self.asMatrix()[i,:])) for i in range(3)]) ) + \
-           'Bunge Eulers / deg: %s' % ('\t'.join(map(str,self.asEulers('bunge',degrees=True))) )
+           'Bunge Eulers / deg: %s' % ('\t'.join(map(str,self.asEulers(degrees=True))) )
 
   def asQuaternion(self):
     return self.quaternion.asList()
 
   def asEulers(self,
-               type = 'bunge',
                degrees = False,
-               standardRange = False):
+              ):
-    return self.quaternion.asEulers(type, degrees, standardRange)
+    return self.quaternion.asEulers(degrees)
   eulers = property(asEulers)
 
   def asRodrigues(self):
@@ -912,13 +853,13 @@ class Orientation:
     """
     if self.symmetry != other.symmetry: raise TypeError('disorientation between different symmetry classes not supported yet.')
 
-    misQ = self.quaternion.conjugated()*other.quaternion
+    misQ = other.quaternion*self.quaternion.conjugated()
     mySymQs    =  self.symmetry.symmetryQuats() if SST else self.symmetry.symmetryQuats()[:1]       # take all or only first sym operation
     otherSymQs = other.symmetry.symmetryQuats()
     
     for i,sA in enumerate(mySymQs):
       for j,sB in enumerate(otherSymQs):
-        theQ = sA.conjugated()*misQ*sB
+        theQ = sB*misQ*sA.conjugated()
         for k in range(2):
           theQ.conjugate()
           breaker = self.symmetry.inFZ(theQ) \
@@ -929,7 +870,7 @@ class Orientation:
 
 # disorientation, own sym, other sym, self-->other: True, self<--other: False
     return (Orientation(quaternion = theQ,symmetry = self.symmetry.lattice),
-            i,j,k == 1)                                                                             
+            i,j, k == 1)                                                                             
 
 
   def inversePole(self,
@@ -939,10 +880,10 @@ class Orientation:
     """Axis rotated according to orientation (using crystal symmetry to ensure location falls into SST)"""
     if SST:                                                                                         # pole requested to be within SST
       for i,q in enumerate(self.symmetry.equivalentQuaternions(self.quaternion)):                   # test all symmetric equivalent quaternions
-        pole = q.conjugated()*axis                                                                  # align crystal direction to axis
+        pole = q*axis                                                                               # align crystal direction to axis
         if self.symmetry.inSST(pole,proper): break                                                  # found SST version
     else:
-      pole = self.quaternion.conjugated()*axis                                                      # align crystal direction to axis
+      pole = self.quaternion*axis                                                                   # align crystal direction to axis
 
     return (pole,i if SST else 0)
 
@@ -951,7 +892,7 @@ class Orientation:
     color = np.zeros(3,'d')
 
     for q in self.symmetry.equivalentQuaternions(self.quaternion):
-      pole = q.conjugated()*axis                                                                    # align crystal direction to axis
+      pole = q*axis                                                                                 # align crystal direction to axis
       inSST,color = self.symmetry.inSST(pole,color=True)
       if inSST: break
 

processing/post/addMises.py

@@ -1,4 +1,4 @@
-#!/usr/bin/env python3
+#!/usr/bin/env python2.7
 # -*- coding: UTF-8 no BOM -*-
 
 import os,sys,math

processing/post/addOrientations.py

@@ -169,7 +169,7 @@ for name in filenames:
     for output in options.output:
       if   output == 'quaternion': table.data_append(o.asQuaternion())
       elif output == 'rodrigues':  table.data_append(o.asRodrigues())
-      elif output == 'eulers':     table.data_append(o.asEulers('Bunge', degrees=options.degrees))
+      elif output == 'eulers':     table.data_append(o.asEulers(degrees=options.degrees))
     outputAlive = table.data_write()                                                                # output processed line
 
 # ------------------------------------------ output finalization -----------------------------------  

processing/post/addSchmidfactors.py

@@ -190,27 +190,27 @@ if np.sum(input) != 1: parser.error('needs exactly one input format.')
                          (options.quaternion,4,'quaternion'),
                         ][np.where(input)[0][0]]                                                    # select input label that was requested
 
-c_direction = np.zeros((len(slipSystems[options.lattice]),3),'f')
+slip_direction = np.zeros((len(slipSystems[options.lattice]),3),'f')
-c_normal    = np.zeros_like(c_direction)
+slip_normal    = np.zeros_like(slip_direction)
 
 
 if options.lattice in latticeChoices[:2]:
-  c_direction = slipSystems[options.lattice][:,:3]
+  slip_direction = slipSystems[options.lattice][:,:3]
-  c_normal    = slipSystems[options.lattice][:,3:]
+  slip_normal    = slipSystems[options.lattice][:,3:]
 elif options.lattice == latticeChoices[2]:
   # convert 4 Miller index notation of hex to orthogonal 3 Miller index notation
-  for i in range(len(c_direction)):
+  for i in range(len(slip_direction)):
-    c_direction[i] = np.array([slipSystems['hex'][i,0]*1.5,
+    slip_direction[i] = np.array([slipSystems['hex'][i,0]*1.5,
                               (slipSystems['hex'][i,0] + 2.*slipSystems['hex'][i,1])*0.5*np.sqrt(3),
                                slipSystems['hex'][i,3]*options.CoverA,
                               ])
-    c_normal[i]    = np.array([slipSystems['hex'][i,4],
+    slip_normal[i]    = np.array([slipSystems['hex'][i,4],
                               (slipSystems['hex'][i,4] + 2.*slipSystems['hex'][i,5])/np.sqrt(3),
                                slipSystems['hex'][i,7]/options.CoverA,
                               ])
 
-c_direction /= np.tile(np.linalg.norm(c_direction,axis=1),(3,1)).T
+slip_direction /= np.tile(np.linalg.norm(slip_direction,axis=1),(3,1)).T
-c_normal    /= np.tile(np.linalg.norm(c_normal   ,axis=1),(3,1)).T
+slip_normal    /= np.tile(np.linalg.norm(slip_normal   ,axis=1),(3,1)).T
 
 # --- loop over input files ------------------------------------------------------------------------
 
@@ -244,7 +244,7 @@ for name in filenames:
                        .format(       id = i+1,
                                   normal = theNormal,
                                direction = theDirection,
-                              ) for i,(theNormal,theDirection) in enumerate(zip(c_normal,c_direction))])
+                              ) for i,(theNormal,theDirection) in enumerate(zip(slip_normal,slip_direction))])
   table.head_write()
 
 # ------------------------------------------ process data ------------------------------------------
@@ -262,9 +262,9 @@ for name in filenames:
     elif inputtype == 'quaternion':
       o = damask.Orientation(quaternion = np.array(list(map(float,table.data[column:column+4]))),)
 
-    rotForce  = o.quaternion.conjugated() * force
+
-    rotNormal = o.quaternion.conjugated() * normal
+    table.data_append(  np.abs(  np.sum(slip_direction * (o.quaternion * force) ,axis=1) \
-    table.data_append(np.abs(np.sum(c_direction*rotForce,axis=1) * np.sum(c_normal*rotNormal,axis=1)))
+                               * np.sum(slip_normal    * (o.quaternion * normal),axis=1)))
     outputAlive = table.data_write()                                                                # output processed line
 
 # ------------------------------------------ output finalization -----------------------------------  

processing/pre/seeds_fromRandom.py

@@ -127,7 +127,7 @@ options.fraction = np.array(options.fraction)
 options.grid = np.array(options.grid)
 gridSize = options.grid.prod()
 
-if options.randomSeed is None: options.randomSeed = int(os.urandom(4).encode('hex'), 16)
+if options.randomSeed is None: options.randomSeed = int(os.urandom(4).hex(), 16)
 np.random.seed(options.randomSeed)                                                                  # init random generators
 random.seed(options.randomSeed)
 

src/IO.f90

@@ -33,7 +33,6 @@ module IO
    IO_write_jobIntFile, &
    IO_read_realFile, &
    IO_read_intFile, &
-   IO_hybridIA, &
    IO_isBlank, &
    IO_getTag, &
    IO_stringPos, &
@@ -583,223 +582,6 @@ logical function IO_abaqus_hasNoPart(fileUnit)
 620 end function IO_abaqus_hasNoPart
 #endif
 
-!--------------------------------------------------------------------------------------------------
-!> @brief hybrid IA sampling of ODFfile
-!--------------------------------------------------------------------------------------------------
-function IO_hybridIA(Nast,ODFfileName)
- use prec, only: &
-   tol_math_check
-
- implicit none
- integer(pInt),                 intent(in)   :: Nast                                                !< number of samples?
- real(pReal), dimension(3,Nast)              :: IO_hybridIA
- character(len=*),              intent(in)   :: ODFfileName                                         !< name of ODF file including total path
-
-!--------------------------------------------------------------------------------------------------
-! math module is not available
- real(pReal),      parameter  :: PI = 3.141592653589793_pReal
- real(pReal),      parameter  :: INRAD = PI/180.0_pReal
-
- integer(pInt) :: i,j,bin,NnonZero,Nset,Nreps,reps,phi1,Phi,phi2
- integer(pInt), allocatable, dimension(:)   :: chunkPos
- integer(pInt), dimension(3)                :: steps                                                !< number of steps in phi1, Phi, and phi2 direction
- integer(pInt), dimension(4)                :: columns                                              !< columns in linearODF file where eulerangles and density are located
- integer(pInt), dimension(:), allocatable   :: binSet
- real(pReal) :: center,sum_dV_V,prob,dg_0,C,lowerC,upperC,rnd
- real(pReal), dimension(2,3)                :: limits                                               !< starting and end values for eulerangles
- real(pReal), dimension(3)                  :: deltas, &                                            !< angular step size in phi1, Phi, and phi2 direction
-                                               eulers                                               !< euler angles when reading from file
- real(pReal), dimension(:,:,:), allocatable :: dV_V
- character(len=65536) :: line, keyword
- integer(pInt)                                :: headerLength
- integer(pInt),               parameter       :: FILEUNIT = 999_pInt
-
- IO_hybridIA = 0.0_pReal                                                                           ! initialize return value for case of error
- write(6,'(/,a,/)',advance='no') ' Using linear ODF file: '//trim(ODFfileName)
- write(6,'(/,a)') ' Eisenlohr et al., Computational Materials Science, 42(4):670–678, 2008'
- write(6,'(a)')   ' https://doi.org/10.1016/j.commatsci.2007.09.015'
-
-
-!--------------------------------------------------------------------------------------------------
-! parse header of ODF file
- call IO_open_file(FILEUNIT,ODFfileName)
- headerLength = 0_pInt
- line=IO_read(FILEUNIT)
- chunkPos = IO_stringPos(line)
- keyword = IO_lc(IO_StringValue(line,chunkPos,2_pInt,.true.))
- if (keyword(1:4) == 'head') then
-   headerLength = IO_intValue(line,chunkPos,1_pInt) + 1_pInt
- else
-   call IO_error(error_ID=156_pInt, ext_msg='no header found')
- endif
-
-!--------------------------------------------------------------------------------------------------
-! figure out columns containing data
- do i = 1_pInt, headerLength-1_pInt
-   line=IO_read(FILEUNIT)
- enddo
- columns = 0_pInt
- chunkPos = IO_stringPos(line)
- do i = 1_pInt, chunkPos(1)
-   select case ( IO_lc(IO_StringValue(line,chunkPos,i,.true.)) )
-     case ('phi1')
-      columns(1) = i
-     case ('phi')
-      columns(2) = i
-     case ('phi2')
-      columns(3) = i
-     case ('intensity')
-      columns(4) = i
-   end select
- enddo
-
- if (any(columns<1)) call IO_error(error_ID = 156_pInt, ext_msg='could not find expected header')
-
-!--------------------------------------------------------------------------------------------------
-! determine limits, number of steps and step size
- limits(1,1:3) = 721.0_pReal
- limits(2,1:3) =  -1.0_pReal
- steps  = 0_pInt
-
- line=IO_read(FILEUNIT)
- do while (trim(line) /= IO_EOF)
-   chunkPos = IO_stringPos(line)
-   eulers=[IO_floatValue(line,chunkPos,columns(1)),&
-           IO_floatValue(line,chunkPos,columns(2)),&
-           IO_floatValue(line,chunkPos,columns(3))]
-   steps = steps + merge(1,0,eulers>limits(2,1:3))
-   limits(1,1:3) = min(limits(1,1:3),eulers)
-   limits(2,1:3) = max(limits(2,1:3),eulers)
-   line=IO_read(FILEUNIT)
- enddo
-
- deltas = (limits(2,1:3)-limits(1,1:3))/real(steps-1_pInt,pReal)
-
- write(6,'(/,a,/,3(2x,f12.4,1x))',advance='no') ' Starting angles / ° = ',limits(1,1:3)
- write(6,'(/,a,/,3(2x,f12.4,1x))',advance='no') ' Ending angles / ° =   ',limits(2,1:3)
- write(6,'(/,a,/,3(2x,f12.4,1x))',advance='no') ' Angular steps / ° =   ',deltas
-
- if (all(abs(limits(1,1:3)) < tol_math_check)) then
-   write(6,'(/,a,/)',advance='no') ' assuming vertex centered data'
-   center = 0.0_pReal                                                                               ! no need to shift
-   if (any(mod(int(limits(2,1:3),pInt),90)==0)) &
-     call IO_error(error_ID = 156_pInt, ext_msg='linear ODF data repeated at right boundary')
- else
-   write(6,'(/,a,/)',advance='no') ' assuming cell centered data'
-   center = 0.5_pReal                                                                               ! shift data by half of a bin
- endif
-
- limits = limits*INRAD
- deltas = deltas*INRAD
-
-!--------------------------------------------------------------------------------------------------
-! read in data
- allocate(dV_V(steps(3),steps(2),steps(1)),source=0.0_pReal)
- sum_dV_V = 0.0_pReal
- dg_0 = deltas(1)*deltas(3)*2.0_pReal*sin(deltas(2)/2.0_pReal)
- NnonZero = 0_pInt
-
- call IO_checkAndRewind(FILEUNIT)                                                                   ! forward
- do i = 1_pInt, headerLength
-   line=IO_read(FILEUNIT)
- enddo
-
- do phi1=1_pInt,steps(1); do Phi=1_pInt,steps(2); do phi2=1_pInt,steps(3)
-   line=IO_read(FILEUNIT)
-   chunkPos = IO_stringPos(line)
-   eulers=[IO_floatValue(line,chunkPos,columns(1)),&                                                ! read in again for consistency check only
-           IO_floatValue(line,chunkPos,columns(2)),&
-           IO_floatValue(line,chunkPos,columns(3))]*INRAD
-   if (any(abs((real([phi1,phi,phi2],pReal) -1.0_pReal + center)*deltas-eulers)>tol_math_check)) &  ! check if data is in expected order (phi2 fast) and correct for Fortran starting at 1
-     call IO_error(error_ID = 156_pInt, ext_msg='linear ODF data not in expected order')
-
-   prob = IO_floatValue(line,chunkPos,columns(4))
-   if (prob > 0.0_pReal) then
-      NnonZero = NnonZero+1_pInt
-      sum_dV_V = sum_dV_V+prob
-    else
-      prob = 0.0_pReal
-    endif
-    dV_V(phi2,Phi,phi1) = prob*dg_0*sin((real(Phi-1_pInt,pReal)+center)*deltas(2))
- enddo; enddo; enddo
- close(FILEUNIT)
- dV_V = dV_V/sum_dV_V                                                                               ! normalize to 1
-
-!--------------------------------------------------------------------------------------------------
-! now fix bounds
- Nset = max(Nast,NnonZero)                                                                          ! if less than non-zero voxel count requested, sample at least that much
- lowerC = 0.0_pReal
- upperC = real(Nset, pReal)
-
- do while (hybridIA_reps(dV_V,steps,upperC) < Nset)
-   lowerC = upperC
-   upperC = upperC*2.0_pReal
- enddo
-
-!--------------------------------------------------------------------------------------------------
-! binary search for best C
- do
-   C = (upperC+lowerC)/2.0_pReal
-   Nreps = hybridIA_reps(dV_V,steps,C)
-   if (abs(upperC-lowerC) < upperC*1.0e-14_pReal) then
-     C = upperC
-     Nreps = hybridIA_reps(dV_V,steps,C)
-     exit
-   elseif (Nreps < Nset) then
-     lowerC = C
-   elseif (Nreps > Nset) then
-     upperC = C
-   else
-     exit
-   endif
- enddo
-
- allocate(binSet(Nreps))
- bin = 0_pInt                                                                                       ! bin counter
- i = 1_pInt                                                                                         ! set counter
- do phi1=1_pInt,steps(1); do Phi=1_pInt,steps(2) ;do phi2=1_pInt,steps(3)
-   reps = nint(C*dV_V(phi2,Phi,phi1), pInt)
-   binSet(i:i+reps-1) = bin
-   bin = bin+1_pInt                                                                                 ! advance bin
-   i = i+reps                                                                                       ! advance set
- enddo; enddo; enddo
-
- do i=1_pInt,Nast
-   if (i < Nast) then
-     call random_number(rnd)
-     j = nint(rnd*real(Nreps-i,pReal)+real(i,pReal)+0.5_pReal,pInt)
-   else
-     j = i
-   endif
-   bin = binSet(j)
-   IO_hybridIA(1,i) = deltas(1)*(real(mod(bin/(steps(3)*steps(2)),steps(1)),pReal)+center)          ! phi1
-   IO_hybridIA(2,i) = deltas(2)*(real(mod(bin/ steps(3)          ,steps(2)),pReal)+center)          ! Phi
-   IO_hybridIA(3,i) = deltas(3)*(real(mod(bin                    ,steps(3)),pReal)+center)          ! phi2
-   binSet(j) = binSet(i)
- enddo
-
- contains
- !--------------------------------------------------------------------------------------------------
- !> @brief counts hybrid IA repetitions
- !--------------------------------------------------------------------------------------------------
- integer(pInt) pure function hybridIA_reps(dV_V,steps,C)
-
-  implicit none
-  integer(pInt), intent(in), dimension(3)                          :: steps                          !< number of bins in Euler space
-  real(pReal),   intent(in), dimension(steps(3),steps(2),steps(1)) :: dV_V                           !< needs description
-  real(pReal),   intent(in)                                        :: C                              !< needs description
-
-  integer(pInt) :: phi1,Phi,phi2
-
-  hybridIA_reps = 0_pInt
-  do phi1=1_pInt,steps(1); do Phi =1_pInt,steps(2); do phi2=1_pInt,steps(3)
-    hybridIA_reps = hybridIA_reps+nint(C*dV_V(phi2,Phi,phi1), pInt)
-  enddo; enddo; enddo
-
- end function hybridIA_reps
-
-end function IO_hybridIA
-
 
 !--------------------------------------------------------------------------------------------------
 !> @brief identifies strings without content
@@ -1758,7 +1540,6 @@ integer(pInt) function IO_verifyIntValue (string,validChars,myName)
                                  validChars, &                                                        !< valid characters in string
                                  myName                                                               !< name of caller function (for debugging)
  integer(pInt)                :: readStatus, invalidWhere
- !character(len=len(trim(string))) :: trimmed does not work with ifort 14.0.1
 
  IO_verifyIntValue = 0_pInt
 
@@ -1788,7 +1569,6 @@ real(pReal) function IO_verifyFloatValue (string,validChars,myName)
                                  myName                                                               !< name of caller function (for debugging)
 
  integer(pInt)                :: readStatus, invalidWhere
- !character(len=len(trim(string))) :: trimmed does not work with ifort 14.0.1
 
  IO_verifyFloatValue = 0.0_pReal
 

src/constitutive.f90

@@ -513,7 +513,7 @@ subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, S6, Fi, ipc, ip, e
    case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType
      of = phasememberAt(ipc,ip,el)
      instance = phase_plasticityInstance(material_phase(ipc,ip,el))
-     call plastic_phenopowerlaw_LpAndItsTangent   (Lp,dLp_dMp, Mp,instance,of)
+     call plastic_phenopowerlaw_LpAndItsTangent   (Lp,dLp_dMp,Mp,instance,of)
 
    case (PLASTICITY_KINEHARDENING_ID) plasticityType
      call plastic_kinehardening_LpAndItsTangent   (Lp,dLp_dMp99, math_Mandel33to6(Mp),ipc,ip,el)
@@ -525,9 +525,9 @@ subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, S6, Fi, ipc, ip, e
      dLp_dMp = math_Plain99to3333(dLp_dMp99)                                                        ! ToDo: We revert here the last statement in plastic_xx_LpAndItsTanget
 
    case (PLASTICITY_DISLOTWIN_ID) plasticityType
-     call plastic_dislotwin_LpAndItsTangent       (Lp,dLp_dMp99, math_Mandel33to6(Mp), &
+     of = phasememberAt(ipc,ip,el)
-                                                   temperature(ho)%p(tme),ipc,ip,el)
+     instance = phase_plasticityInstance(material_phase(ipc,ip,el))
-     dLp_dMp = math_Plain99to3333(dLp_dMp99)                                                        ! ToDo: We revert here the last statement in plastic_xx_LpAndItsTanget
+     call plastic_dislotwin_LpAndItsTangent       (Lp,dLp_dMp,Mp,temperature(ho)%p(tme),instance,of)
 
    case (PLASTICITY_DISLOUCLA_ID) plasticityType
      call plastic_disloucla_LpAndItsTangent       (Lp,dLp_dMp99, math_Mandel33to6(Mp), &
@@ -922,8 +922,9 @@ subroutine constitutive_collectDotState(S6, FeArray, Fi, FpArray, subdt, subfrac
      call plastic_kinehardening_dotState(math_Mandel33to6(Mp),ipc,ip,el)
 
    case (PLASTICITY_DISLOTWIN_ID) plasticityType
-     call plastic_dislotwin_dotState    (math_Mandel33to6(Mp),temperature(ho)%p(tme), &
+     of = phasememberAt(ipc,ip,el)
-                                         ipc,ip,el)
+     instance = phase_plasticityInstance(material_phase(ipc,ip,el))
+     call plastic_dislotwin_dotState    (Mp,temperature(ho)%p(tme),instance,of)
 
    case (PLASTICITY_DISLOUCLA_ID) plasticityType
      call plastic_disloucla_dotState    (math_Mandel33to6(Mp),temperature(ho)%p(tme), &
@@ -1135,20 +1136,27 @@ function constitutive_postResults(S6, Fi, FeArray, ipc, ip, el)
    case (PLASTICITY_ISOTROPIC_ID) plasticityType
      constitutive_postResults(startPos:endPos) = &
        plastic_isotropic_postResults(S6,ipc,ip,el)
+
    case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType
      of = phasememberAt(ipc,ip,el)
      instance = phase_plasticityInstance(material_phase(ipc,ip,el))
      constitutive_postResults(startPos:endPos) = &
        plastic_phenopowerlaw_postResults(Mp,instance,of)
+
    case (PLASTICITY_KINEHARDENING_ID) plasticityType
      constitutive_postResults(startPos:endPos) = &
        plastic_kinehardening_postResults(S6,ipc,ip,el)
+
    case (PLASTICITY_DISLOTWIN_ID) plasticityType
+     of = phasememberAt(ipc,ip,el)
+     instance = phase_plasticityInstance(material_phase(ipc,ip,el))
      constitutive_postResults(startPos:endPos) = &
-       plastic_dislotwin_postResults(S6,temperature(ho)%p(tme),ipc,ip,el)
+       plastic_dislotwin_postResults(Mp,temperature(ho)%p(tme),instance,of)
+
    case (PLASTICITY_DISLOUCLA_ID) plasticityType
      constitutive_postResults(startPos:endPos) = &
        plastic_disloucla_postResults(S6,temperature(ho)%p(tme),ipc,ip,el)
+
    case (PLASTICITY_NONLOCAL_ID) plasticityType
      constitutive_postResults(startPos:endPos) = &
        plastic_nonlocal_postResults (S6,FeArray,ip,el)

src/lattice.f90

@@ -72,7 +72,7 @@ module lattice
 !--------------------------------------------------------------------------------------------------
 ! face centered cubic
  integer(pInt), dimension(LATTICE_maxNslipFamily), parameter, public :: &
-   LATTICE_fcc_NslipSystem = int([12, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],pInt)                     !< # of slip systems per family for fcc
+   LATTICE_fcc_NslipSystem = int([12, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],pInt)                     !< # of slip systems per family for fcc
 
  integer(pInt), dimension(LATTICE_maxNtwinFamily), parameter, public :: &
    LATTICE_fcc_NtwinSystem = int([12, 0, 0, 0],pInt)                                                !< # of twin systems per family for fcc
@@ -104,11 +104,19 @@ module lattice
       1, 1, 0,     1,-1,-1, &                                                                       ! A6
       0, 1, 1,    -1, 1,-1, &                                                                       ! D1
       1, 0,-1,    -1, 1,-1, &                                                                       ! D4
-     -1,-1, 0,    -1, 1,-1  &                                                                       ! D6
+     -1,-1, 0,    -1, 1,-1, &                                                                       ! D6
+     ! Slip system <110>{110}
+      1, 1, 0,     1,-1, 0, &
+      1,-1, 0,     1, 1, 0, &
+      1, 0, 1,     1, 0,-1, &
+      1, 0,-1,     1, 0, 1, &
+      0, 1, 1,     0, 1,-1, &
+      0, 1,-1,     0, 1, 1  &
      ],pReal),shape(LATTICE_FCC_SYSTEMSLIP))                                                        !< Slip system <110>{111} directions. Sorted according to Eisenlohr & Hantcherli
 
- character(len=*), dimension(1), parameter, public :: LATTICE_FCC_SLIPFAMILY_NAME = &
+ character(len=*), dimension(2), parameter, public :: LATTICE_FCC_SLIPFAMILY_NAME = &
-   ['<0 1 -1>{1 1 1}']
+   ['<0 1 -1>{1 1 1}', &
+    '<0 1 -1>{0 1 1}']
 
  real(pReal), dimension(3+3,LATTICE_fcc_Ntwin), parameter, private :: &
    LATTICE_fcc_systemTwin = reshape(real( [&
@@ -166,25 +174,38 @@ module lattice
 
  integer(pInt), dimension(LATTICE_fcc_Nslip,lattice_fcc_Nslip), parameter, public :: &
    LATTICE_fcc_interactionSlipSlip = reshape(int( [&
-     1,2,2,4,6,5,3,5,5,4,5,6, &  ! ---> slip
+      1, 2, 2, 4, 6, 5, 3, 5, 5, 4, 5, 6,  9,10, 9,10,11,12, &  ! ---> slip
-     2,1,2,6,4,5,5,4,6,5,3,5, &  ! |
+      2, 1, 2, 6, 4, 5, 5, 4, 6, 5, 3, 5,  9,10,11,12, 9,10, &  ! |
-     2,2,1,5,5,3,5,6,4,6,5,4, &  ! |
+      2, 2, 1, 5, 5, 3, 5, 6, 4, 6, 5, 4, 11,12, 9,10, 9,10, &  ! |
-     4,6,5,1,2,2,4,5,6,3,5,5, &  ! v slip
+      4, 6, 5, 1, 2, 2, 4, 5, 6, 3, 5, 5,  9,10,10, 9,12,11, &  ! v slip
-     6,4,5,2,1,2,5,3,5,5,4,6, &
+      6, 4, 5, 2, 1, 2, 5, 3, 5, 5, 4, 6,  9,10,12,11,10, 9, &
-     5,5,3,2,2,1,6,5,4,5,6,4, &
+      5, 5, 3, 2, 2, 1, 6, 5, 4, 5, 6, 4, 11,12,10, 9,10, 9, &
-     3,5,5,4,5,6,1,2,2,4,6,5, &
+      3, 5, 5, 4, 5, 6, 1, 2, 2, 4, 6, 5, 10, 9,10, 9,11,12, &
-     5,4,6,5,3,5,2,1,2,6,4,5, &
+      5, 4, 6, 5, 3, 5, 2, 1, 2, 6, 4, 5, 10, 9,12,11, 9,10, &
-     5,6,4,6,5,4,2,2,1,5,5,3, &
+      5, 6, 4, 6, 5, 4, 2, 2, 1, 5, 5, 3, 12,11,10, 9, 9,10, &
-     4,5,6,3,5,5,4,6,5,1,2,2, &
+      4, 5, 6, 3, 5, 5, 4, 6, 5, 1, 2, 2, 10, 9, 9,10,12,11, &
-     5,3,5,5,4,6,6,4,5,2,1,2, &
+      5, 3, 5, 5, 4, 6, 6, 4, 5, 2, 1, 2, 10, 9,11,12,10, 9, &
-     6,5,4,5,6,4,5,5,3,2,2,1  &
+      6, 5, 4, 5, 6, 4, 5, 5, 3, 2, 2, 1, 12,11, 9,10,10, 9, &
-     ],pInt),shape(LATTICE_FCC_INTERACTIONSLIPSLIP),order=[2,1])                                    !< Slip--slip interaction types for fcc
+     
+      9, 9,11, 9, 9,11,10,10,12,10,10,12,  1, 7, 8, 8, 8, 8, &  
+     10,10,12,10,10,12, 9, 9,11, 9, 9,11,  7, 1, 8, 8, 8, 8, &  
+      9,11, 9,10,12,10,10,12,10, 9,11, 9,  8, 8, 1, 7, 8, 8, &   
+     10,12,10, 9,11, 9, 9,11, 9,10,12,10,  8, 8, 7, 1, 8, 8, &  
+     11, 9, 9,12,10,10,11, 9, 9,12,10,10,  8, 8, 8, 8, 1, 7, &  
+     12,10,10,11, 9, 9,12,10,10,11, 9, 9,  8, 8, 8, 8, 7, 1  &  
+     ],pInt),[LATTICE_fcc_Nslip,LATTICE_fcc_Nslip],order=[2,1])                                     !< Slip--slip interaction types for fcc
                                                                                                     !< 1: self interaction
                                                                                                     !< 2: coplanar interaction
                                                                                                     !< 3: collinear interaction
                                                                                                     !< 4: Hirth locks
                                                                                                     !< 5: glissile junctions
                                                                                                     !< 6: Lomer locks
+                                                                                                    !< 7: crossing (similar to Hirth locks in <110>{111} for two {110} planes)
+                                                                                                    !< 8: similar to Lomer locks in <110>{111} for two {110} planes
+                                                                                                    !< 9: similar to Lomer locks in <110>{111} btw one {110} and one {111} plane
+                                                                                                    !<10: similar to glissile junctions in <110>{111} btw one {110} and one {111} plane
+                                                                                                    !<11: crossing btw one {110} and one {111} plane
+                                                                                                    !<12: collinear btw one {110} and one {111} plane
  integer(pInt), dimension(LATTICE_fcc_Nslip,LATTICE_fcc_Ntwin), parameter, public :: &
    LATTICE_fcc_interactionSlipTwin = reshape(int( [&
      1,1,1,3,3,3,2,2,2,3,3,3, & ! ---> twin
@@ -198,7 +219,14 @@ module lattice
      3,3,3,3,3,3,1,1,1,2,2,2, &
      3,3,3,2,2,2,3,3,3,1,1,1, &
      2,2,2,3,3,3,3,3,3,1,1,1, &
-     3,3,3,3,3,3,2,2,2,1,1,1  &
+     3,3,3,3,3,3,2,2,2,1,1,1, &
+
+     4,4,4,4,4,4,4,4,4,4,4,4, &
+     4,4,4,4,4,4,4,4,4,4,4,4, &
+     4,4,4,4,4,4,4,4,4,4,4,4, &
+     4,4,4,4,4,4,4,4,4,4,4,4, &
+     4,4,4,4,4,4,4,4,4,4,4,4, &
+     4,4,4,4,4,4,4,4,4,4,4,4  &
      ],pInt),shape(LATTICE_FCC_INTERACTIONSLIPTWIN),order=[2,1])                                    !< Slip--twin interaction types for fcc
                                                                                                     !< 1: coplanar interaction
                                                                                                     !< 2: screw trace between slip system and twin habit plane (easy cross slip)
@@ -235,7 +263,14 @@ module lattice
      3,3,3,3,3,3,1,1,1,2,2,2, &
      3,3,3,2,2,2,3,3,3,1,1,1, &
      2,2,2,3,3,3,3,3,3,1,1,1, &
-     3,3,3,3,3,3,2,2,2,1,1,1  &
+     3,3,3,3,3,3,2,2,2,1,1,1, &
+     
+     4,4,4,4,4,4,4,4,4,4,4,4, &
+     4,4,4,4,4,4,4,4,4,4,4,4, &
+     4,4,4,4,4,4,4,4,4,4,4,4, &
+     4,4,4,4,4,4,4,4,4,4,4,4, &
+     4,4,4,4,4,4,4,4,4,4,4,4, &
+     4,4,4,4,4,4,4,4,4,4,4,4  &
      ],pInt),shape(LATTICE_FCCTOHEX_INTERACTIONSLIPTRANS),order=[2,1])                              !< Slip--trans interaction types for fcc
 
  integer(pInt), dimension(LATTICE_fcc_Ntrans,LATTICE_fcc_Nslip), parameter, public :: &

src/material.f90

@@ -169,18 +169,13 @@ module material
    homogenization_maxNgrains                                                                        !< max number of grains in any USED homogenization
 
  integer(pInt), dimension(:), allocatable, public, protected :: &
-   material_homogenizationAt, &                                                                     !< homogenization ID of each element (copy of mesh_homogenizationAt)
    phase_Nsources, &                                                                                !< number of source mechanisms active in each phase
    phase_Nkinematics, &                                                                             !< number of kinematic mechanisms active in each phase
    phase_NstiffnessDegradations, &                                                                  !< number of stiffness degradation mechanisms active in each phase
    phase_Noutput, &                                                                                 !< number of '(output)' items per phase
    phase_elasticityInstance, &                                                                      !< instance of particular elasticity of each phase
-   phase_plasticityInstance                                                                         !< instance of particular plasticity of each phase
+   phase_plasticityInstance, &                                                                      !< instance of particular plasticity of each phase
-
+   crystallite_Noutput, &                                                                           !< number of '(output)' items per crystallite setting
- integer(pInt), dimension(:), allocatable, public, protected :: &
-   crystallite_Noutput                                                                              !< number of '(output)' items per crystallite setting
-
- integer(pInt), dimension(:), allocatable, public, protected :: &
    homogenization_Ngrains, &                                                                        !< number of grains in each homogenization
    homogenization_Noutput, &                                                                        !< number of '(output)' items per homogenization
    homogenization_typeInstance, &                                                                   !< instance of particular type of each homogenization
@@ -189,7 +184,7 @@ module material
    vacancyflux_typeInstance, &                                                                      !< instance of particular type of each vacancy flux
    porosity_typeInstance, &                                                                         !< instance of particular type of each porosity model
    hydrogenflux_typeInstance, &                                                                     !< instance of particular type of each hydrogen flux
-   microstructure_crystallite                                                                       !< crystallite setting ID of each microstructure
+   microstructure_crystallite                                                                       !< crystallite setting ID of each microstructure ! DEPRECATED !!!!
 
  real(pReal), dimension(:), allocatable, public, protected :: &
    thermal_initialT, &                                                                              !< initial temperature per each homogenization
@@ -198,12 +193,27 @@ module material
    porosity_initialPhi, &                                                                           !< initial posority per each homogenization
    hydrogenflux_initialCh                                                                           !< initial hydrogen concentration per each homogenization
 
+! NEW MAPPINGS 
+ integer(pInt), dimension(:), allocatable, public, protected :: &
+   material_homogenizationAt, &                                                                     !< homogenization ID of each element (copy of mesh_homogenizationAt)
+   material_homogenizationMemberAt, &                                                               !< position of the element within its homogenization instance
+   material_aggregateAt, &                                                                          !< aggregate ID of each element FUTURE USE FOR OUTPUT
+   material_aggregatMemberAt                                                                        !< position of the element within its aggregate instance FUTURE USE FOR OUTPUT
+ integer(pInt), dimension(:,:), allocatable, public, protected :: &
+   material_phaseAt, &                                                                              !< phase ID of each element
+   material_phaseMemberAt, &                                                                        !< position of the element within its phase instance
+   material_crystalliteAt, &                                                                        !< crystallite ID of each element CURRENTLY NOT PER CONSTITUTENT
+   material_crystalliteMemberAt                                                                     !< position of the element within its crystallite instance CURRENTLY NOT PER CONSTITUTENT
+! END NEW MAPPINGS
+ 
+! DEPRECATED: use material_phaseAt
  integer(pInt), dimension(:,:,:), allocatable, public :: &
    material_phase                                                                                   !< phase (index) of each grain,IP,element
-! BEGIN DEPRECATED: use material_homogenizationAt
+! DEPRECATED: use material_homogenizationAt
  integer(pInt), dimension(:,:), allocatable, public :: &
    material_homog                                                                                   !< homogenization (index) of each IP,element
 ! END DEPRECATED
+
  type(tPlasticState), allocatable, dimension(:), public :: &
    plasticState
  type(tSourceState),  allocatable, dimension(:), public :: &
@@ -227,10 +237,6 @@ module material
    microstructure_elemhomo, &                                                                       !< flag to indicate homogeneous microstructure distribution over element's IPs
    phase_localPlasticity                                                                            !< flags phases with local constitutive law
 
-
- character(len=65536), dimension(:), allocatable, private :: &
-   texture_ODFfile                                                                                  !< name of each ODF file
-
  integer(pInt), private :: &
    microstructure_maxNconstituents, &                                                               !< max number of constituents in any phase
    texture_maxNgauss, &                                                                             !< max number of Gauss components in any texture
@@ -258,11 +264,13 @@ module material
  logical, dimension(:), allocatable, private :: &
    homogenization_active
 
+! BEGIN DEPRECATED
  integer(pInt), dimension(:,:,:), allocatable, public :: phaseAt                                    !< phase ID of every material point (ipc,ip,el)
  integer(pInt), dimension(:,:,:), allocatable, public :: phasememberAt                              !< memberID of given phase at every material point (ipc,ip,el)
- integer(pInt), dimension(:,:,:), allocatable, public, target :: mappingCrystallite
+
  integer(pInt), dimension(:,:,:), allocatable, public, target :: mappingHomogenization              !< mapping from material points to offset in heterogenous state/field
- integer(pInt), dimension(:,:),   allocatable, public, target :: mappingHomogenizationConst         !< mapping from material points to offset in constant state/field
+ integer(pInt), dimension(:,:),   allocatable, private, target :: mappingHomogenizationConst         !< mapping from material points to offset in constant state/field
+! END DEPRECATED
 
  type(tHomogMapping), allocatable, dimension(:), public :: &
    thermalMapping, &                                                                                !< mapping for thermal state/fields
@@ -283,6 +291,7 @@ module material
 
  public :: &
    material_init, &
+   material_allocatePlasticState, &
    ELASTICITY_hooke_ID ,&
    PLASTICITY_none_ID, &
    PLASTICITY_isotropic_ID, &
@@ -367,7 +376,6 @@ subroutine material_init()
  use mesh, only: &
    mesh_homogenizationAt, &
    mesh_NipsPerElem, &
-   mesh_maxNips, &
    mesh_NcpElems, &
    FE_geomtype
 
@@ -377,11 +385,10 @@ subroutine material_init()
  integer(pInt) :: &
   g, &                                                                                              !< grain number
   i, &                                                                                              !< integration point number
-  e, &                                                                                              !< element number
+  e                                                                                                 !< element number
-  phase
+ integer(pInt), dimension(:), allocatable :: &
- integer(pInt), dimension(:), allocatable :: ConstitutivePosition
+  CounterPhase, &
- integer(pInt), dimension(:), allocatable :: CrystallitePosition
+  CounterHomogenization
- integer(pInt), dimension(:), allocatable :: HomogenizationPosition
 
  myDebug = debug_level(debug_material)
 
@@ -472,30 +479,34 @@ subroutine material_init()
 
  call material_populateGrains
 
- allocate(phaseAt                   (  homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems),source=0_pInt)
+! BEGIN DEPRECATED
- allocate(phasememberAt             (  homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems),source=0_pInt)
+ allocate(phaseAt                   (  homogenization_maxNgrains,mesh_nIPsPerElem,mesh_NcpElems),source=0_pInt)
- allocate(mappingHomogenization     (2,                          mesh_maxNips,mesh_NcpElems),source=0_pInt)
+ allocate(phasememberAt             (  homogenization_maxNgrains,mesh_nIPsPerElem,mesh_NcpElems),source=0_pInt)
- allocate(mappingCrystallite        (2,homogenization_maxNgrains,             mesh_NcpElems),source=0_pInt)
+ allocate(mappingHomogenization     (2,                          mesh_nIPsPerElem,mesh_NcpElems),source=0_pInt)
- allocate(mappingHomogenizationConst(                            mesh_maxNips,mesh_NcpElems),source=1_pInt)
+ allocate(mappingHomogenizationConst(                            mesh_nIPsPerElem,mesh_NcpElems),source=1_pInt)
+! END DEPRECATED
 
- allocate(ConstitutivePosition  (size(config_phase)),         source=0_pInt)
+ allocate(material_homogenizationAt,source=mesh_homogenizationAt)
- allocate(HomogenizationPosition(size(config_homogenization)),source=0_pInt)
+ allocate(CounterPhase         (size(config_phase)),         source=0_pInt)
- allocate(CrystallitePosition   (size(config_phase)),         source=0_pInt)
+ allocate(CounterHomogenization(size(config_homogenization)),source=0_pInt)
 
- ElemLoop:do e = 1_pInt,mesh_NcpElems
+! BEGIN DEPRECATED
+ do e = 1_pInt,mesh_NcpElems
  myHomog = mesh_homogenizationAt(e)
-   IPloop:do i = 1_pInt, mesh_NipsPerElem
+   do i = 1_pInt, mesh_NipsPerElem
-     HomogenizationPosition(myHomog) = HomogenizationPosition(myHomog) + 1_pInt
+     CounterHomogenization(myHomog) = CounterHomogenization(myHomog) + 1_pInt
-     mappingHomogenization(1:2,i,e) = [HomogenizationPosition(myHomog),myHomog]
+     mappingHomogenization(1:2,i,e) = [CounterHomogenization(myHomog),myHomog]
-     GrainLoop:do g = 1_pInt,homogenization_Ngrains(myHomog)
+     do g = 1_pInt,homogenization_Ngrains(myHomog)
-       phase = material_phase(g,i,e)
+       myPhase = material_phase(g,i,e)
-       ConstitutivePosition(phase) = ConstitutivePosition(phase)+1_pInt                             ! not distinguishing between instances of same phase
+       CounterPhase(myPhase) = CounterPhase(myPhase)+1_pInt                             ! not distinguishing between instances of same phase
-       phaseAt(g,i,e)              = phase
+       phaseAt(g,i,e)              = myPhase
-       phasememberAt(g,i,e)        = ConstitutivePosition(phase)
+       phasememberAt(g,i,e)        = CounterPhase(myPhase)
-     enddo GrainLoop
+     enddo
-   enddo IPloop
+   enddo
- enddo ElemLoop
+ enddo
+! END DEPRECATED
 
+! REMOVE !!!!!
 ! hack needed to initialize field values used during constitutive and crystallite initializations
  do myHomog = 1,size(config_homogenization)
    thermalMapping     (myHomog)%p => mappingHomogenizationConst
@@ -512,7 +523,7 @@ subroutine material_init()
    allocate(vacancyConcRate (myHomog)%p(1), source=0.0_pReal)
    allocate(hydrogenConcRate(myHomog)%p(1), source=0.0_pReal)
  enddo
- 
+
 end subroutine material_init
 
 
@@ -936,9 +947,7 @@ subroutine material_parseTexture
  integer(pInt) :: section, gauss, fiber, j, t, i
  character(len=65536), dimension(:), allocatable ::  strings                                     ! Values for given key in material config 
  integer(pInt), dimension(:), allocatable :: chunkPos
- character(len=65536) :: tag
 
- allocate(texture_ODFfile(size(config_texture))); texture_ODFfile=''
  allocate(texture_symmetry(size(config_texture)), source=1_pInt)
  allocate(texture_Ngauss(size(config_texture)),   source=0_pInt)
  allocate(texture_Nfiber(size(config_texture)),   source=0_pInt)
@@ -984,9 +993,6 @@ subroutine material_parseTexture
      if(dNeq(math_det33(texture_transformation(1:3,1:3,t)),1.0_pReal)) call IO_error(157_pInt,t)
    endif
 
-   tag=''
-   texture_ODFfile(t) = config_texture(t)%getString('hybridia',defaultVal=tag)
-
    if (config_texture(t)%keyExists('symmetry')) then
      select case (config_texture(t)%getString('symmetry'))
        case('orthotropic')
@@ -1069,10 +1075,59 @@ subroutine material_parseTexture
 end subroutine material_parseTexture
 
 
+!--------------------------------------------------------------------------------------------------
+!> @brief allocates the plastic state of a phase
+!--------------------------------------------------------------------------------------------------
+subroutine material_allocatePlasticState(phase,NofMyPhase,sizeState,sizeDotState,sizeDeltaState,&
+                                         Nslip,Ntwin,Ntrans)
+ use numerics, only: &
+   numerics_integrator2 => numerics_integrator                                                      ! compatibility hack
+
+ implicit none
+ integer(pInt), intent(in) :: &
+   phase, &
+   NofMyPhase, &
+   sizeState, &
+   sizeDotState, &
+   sizeDeltaState, &
+   Nslip, &
+   Ntwin, &
+   Ntrans
+ integer(pInt) :: numerics_integrator                                                               ! compatibility hack
+ numerics_integrator = numerics_integrator2(1)                                                      ! compatibility hack
+
+ plasticState(phase)%sizeState      = sizeState
+ plasticState(phase)%sizeDotState   = sizeDotState
+ plasticState(phase)%sizeDeltaState = sizeDeltaState
+ plasticState(phase)%Nslip = Nslip
+ plasticState(phase)%Ntwin = Ntwin
+ plasticState(phase)%Ntrans= Ntrans
+
+ allocate(plasticState(phase)%aTolState           (sizeState),                source=0.0_pReal)
+ allocate(plasticState(phase)%state0              (sizeState,NofMyPhase),     source=0.0_pReal)
+ allocate(plasticState(phase)%partionedState0     (sizeState,NofMyPhase),     source=0.0_pReal)
+ allocate(plasticState(phase)%subState0           (sizeState,NofMyPhase),     source=0.0_pReal)
+ allocate(plasticState(phase)%state               (sizeState,NofMyPhase),     source=0.0_pReal)
+
+ allocate(plasticState(phase)%dotState            (sizeDotState,NofMyPhase),  source=0.0_pReal)
+ if (numerics_integrator == 1_pInt) then
+   allocate(plasticState(phase)%previousDotState  (sizeDotState,NofMyPhase),  source=0.0_pReal)
+   allocate(plasticState(phase)%previousDotState2 (sizeDotState,NofMyPhase),  source=0.0_pReal)
+ endif
+ if (numerics_integrator == 4_pInt) &
+   allocate(plasticState(phase)%RK4dotState       (sizeDotState,NofMyPhase),  source=0.0_pReal)
+ if (numerics_integrator == 5_pInt) &
+   allocate(plasticState(phase)%RKCK45dotState  (6,sizeDotState,NofMyPhase),  source=0.0_pReal)
+
+ allocate(plasticState(phase)%deltaState        (sizeDeltaState,NofMyPhase),  source=0.0_pReal)
+
+end subroutine material_allocatePlasticState
+
+
 !--------------------------------------------------------------------------------------------------
 !> @brief populates the grains
 !> @details populates the grains by identifying active microstructure/homogenization pairs,
-!! calculates the volume of the grains and deals with texture components and hybridIA
+!! calculates the volume of the grains and deals with texture components
 !--------------------------------------------------------------------------------------------------
 subroutine material_populateGrains
  use prec, only: &
@@ -1091,7 +1146,6 @@ subroutine material_populateGrains
    mesh_elemType, &
    mesh_homogenizationAt, &
    mesh_microstructureAt, &
-   mesh_maxNips, &
    mesh_NcpElems, &
    mesh_ipVolume, &
    FE_geomtype
@@ -1102,8 +1156,7 @@ subroutine material_populateGrains
    homogenization_name, &
    microstructure_name 
  use IO, only: &
-   IO_error, &
+   IO_error
-   IO_hybridIA
  use debug, only: &
    debug_level, &
    debug_material, &
@@ -1131,12 +1184,11 @@ subroutine material_populateGrains
 
  myDebug = debug_level(debug_material)
 
- allocate(material_volume(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems),       source=0.0_pReal)
+ allocate(material_volume(homogenization_maxNgrains,mesh_nIPsPerElem,mesh_NcpElems),       source=0.0_pReal)
- allocate(material_phase(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems),        source=0_pInt)
+ allocate(material_phase(homogenization_maxNgrains,mesh_nIPsPerElem,mesh_NcpElems),        source=0_pInt)
- allocate(material_homog(mesh_maxNips,mesh_NcpElems),                                  source=0_pInt)
+ allocate(material_homog(mesh_nIPsPerElem,mesh_NcpElems),                                  source=0_pInt)
- allocate(material_homogenizationAt,source=mesh_homogenizationAt)
+ allocate(material_texture(homogenization_maxNgrains,mesh_nIPsPerElem,mesh_NcpElems),      source=0_pInt)
- allocate(material_texture(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems),      source=0_pInt)
+ allocate(material_EulerAngles(3,homogenization_maxNgrains,mesh_nIPsPerElem,mesh_NcpElems),source=0.0_pReal)
- allocate(material_EulerAngles(3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems),source=0.0_pReal)
 
  allocate(Ngrains(size(config_homogenization),size(config_microstructure)),            source=0_pInt)
  allocate(Nelems (size(config_homogenization),size(config_microstructure)),            source=0_pInt)
@@ -1280,39 +1332,31 @@ subroutine material_populateGrains
                                    real(texture_symmetry(textureID),pReal),pInt)                    ! max number of unique orientations (excl. symmetry)
 
 !--------------------------------------------------------------------------------------------------
-! ...has texture components
+! has texture components
-         if (texture_ODFfile(textureID) == '') then
+         gauss: do t = 1_pInt,texture_Ngauss(textureID)                                             ! loop over Gauss components
-           gauss: do t = 1_pInt,texture_Ngauss(textureID)                                           ! loop over Gauss components
+           do g = 1_pInt,int(real(myNorientations,pReal)*texture_Gauss(5,t,textureID),pInt)         ! loop over required grain count
-             do g = 1_pInt,int(real(myNorientations,pReal)*texture_Gauss(5,t,textureID),pInt)       ! loop over required grain count
+             orientationOfGrain(:,grain+constituentGrain+g) = &
-               orientationOfGrain(:,grain+constituentGrain+g) = &
+               math_sampleGaussOri(texture_Gauss(1:3,t,textureID),&
-                 math_sampleGaussOri(texture_Gauss(1:3,t,textureID),&
+                                   texture_Gauss(  4,t,textureID))
-                                     texture_Gauss(  4,t,textureID))
+           enddo
-             enddo
+           constituentGrain = &
-             constituentGrain = &
+           constituentGrain + int(real(myNorientations,pReal)*texture_Gauss(5,t,textureID))         ! advance counter for grains of current constituent
-             constituentGrain + int(real(myNorientations,pReal)*texture_Gauss(5,t,textureID))       ! advance counter for grains of current constituent
+         enddo gauss
-           enddo gauss
 
-           fiber: do t = 1_pInt,texture_Nfiber(textureID)                                           ! loop over fiber components
+         fiber: do t = 1_pInt,texture_Nfiber(textureID)                                             ! loop over fiber components
-             do g = 1_pInt,int(real(myNorientations,pReal)*texture_Fiber(6,t,textureID),pInt)       ! loop over required grain count
+           do g = 1_pInt,int(real(myNorientations,pReal)*texture_Fiber(6,t,textureID),pInt)         ! loop over required grain count
-               orientationOfGrain(:,grain+constituentGrain+g) = &
+             orientationOfGrain(:,grain+constituentGrain+g) = &
-                 math_sampleFiberOri(texture_Fiber(1:2,t,textureID),&
+               math_sampleFiberOri(texture_Fiber(1:2,t,textureID),&
-                                     texture_Fiber(3:4,t,textureID),&
+                                   texture_Fiber(3:4,t,textureID),&
-                                     texture_Fiber(  5,t,textureID))
+                                   texture_Fiber(  5,t,textureID))
-             enddo
+           enddo
-             constituentGrain = &
+           constituentGrain = &
-             constituentGrain + int(real(myNorientations,pReal)*texture_fiber(6,t,textureID),pInt)  ! advance counter for grains of current constituent
+           constituentGrain + int(real(myNorientations,pReal)*texture_fiber(6,t,textureID),pInt)    ! advance counter for grains of current constituent
-           enddo fiber
+         enddo fiber
 
-           random: do constituentGrain = constituentGrain+1_pInt,myNorientations                    ! fill remainder with random
+         random: do constituentGrain = constituentGrain+1_pInt,myNorientations                      ! fill remainder with random
-              orientationOfGrain(:,grain+constituentGrain) = math_sampleRandomOri()
+            orientationOfGrain(:,grain+constituentGrain) = math_sampleRandomOri()
-           enddo random
+         enddo random
-!--------------------------------------------------------------------------------------------------
-! ...has hybrid IA
-         else
-           orientationOfGrain(1:3,grain+1_pInt:grain+myNorientations) = &
-                                            IO_hybridIA(myNorientations,texture_ODFfile(textureID))
-           if (all(dEq(orientationOfGrain(1:3,grain+1_pInt),-1.0_pReal))) call IO_error(156_pInt)
-         endif
 
 !--------------------------------------------------------------------------------------------------
 ! ...texture transformation
@@ -1429,12 +1473,7 @@ subroutine material_populateGrains
    enddo microstructureLoop
  enddo homogenizationLoop
 
- deallocate(volumeOfGrain)
- deallocate(phaseOfGrain)
- deallocate(textureOfGrain)
- deallocate(orientationOfGrain)
  deallocate(texture_transformation)
- deallocate(Nelems)
  deallocate(elemsOfHomogMicro)
  call config_deallocate('material.config/microstructure')
 

src/plastic_phenopowerlaw.f90

@@ -32,10 +32,7 @@ module plastic_phenopowerlaw
      totalvolfrac_twin_ID
  end enum
 
- type, private :: tParameters                                                                       !< container type for internal constitutive parameters
+ type, private :: tParameters
-   integer(pInt) :: &
-     totalNslip, &
-     totalNtwin
    real(pReal) :: &
      gdot0_slip, &                                                                                  !< reference shear strain rate for slip
      gdot0_twin, &                                                                                  !< reference shear strain rate for twin
@@ -50,45 +47,48 @@ module plastic_phenopowerlaw
      h0_TwinSlip, &                                                                                 !< reference hardening twin - slip
      h0_TwinTwin, &                                                                                 !< reference hardening twin - twin
      a_slip, &
-     aTolResistance, &                                                                              ! default absolute tolerance 1 Pa
+     aTolResistance, &                                                                              !< absolute tolerance for integration of xi
-     aTolShear, &                                                                                   ! default absolute tolerance 1e-6
+     aTolShear, &                                                                                   !< absolute tolerance for integration of gamma
-     aTolTwinfrac                                                                                   ! default absolute tolerance 1e-6
+     aTolTwinfrac                                                                                   !< absolute tolerance for integration of f
-   integer(pInt), dimension(:),   allocatable :: &
+   real(pReal),                 allocatable, dimension(:) :: &
-     Nslip, &                                                                                       !< active number of slip systems per family
-     Ntwin                                                                                          !< active number of twin systems per family
-   real(pReal),   dimension(:),   allocatable :: &
      xi_slip_0, &                                                                                   !< initial critical shear stress for slip
      xi_twin_0, &                                                                                   !< initial critical shear stress for twin
      xi_slip_sat, &                                                                                 !< maximum critical shear stress for slip
      nonSchmidCoeff, &
      H_int, &                                                                                       !< per family hardening activity (optional) !ToDo: Better name!
      gamma_twin_char                                                                                !< characteristic shear for twins
-   real(pReal),   dimension(:,:),   allocatable :: &
+   real(pReal),                 allocatable, dimension(:,:) :: &
      interaction_SlipSlip, &                                                                        !< slip resistance from slip activity
      interaction_SlipTwin, &                                                                        !< slip resistance from twin activity
      interaction_TwinSlip, &                                                                        !< twin resistance from slip activity
      interaction_TwinTwin                                                                           !< twin resistance from twin activity
-   real(pReal),   dimension(:,:,:),   allocatable :: &
+   real(pReal),                 allocatable, dimension(:,:,:) :: &
      Schmid_slip, &
      Schmid_twin, &
      nonSchmid_pos, &
      nonSchmid_neg
- integer(kind(undefined_ID)),         dimension(:),   allocatable          :: &
+   integer(pInt) :: &
-   outputID                                                                                         !< ID of each post result output
+     totalNslip, &                                                                                  !< total number of active slip system
- end type
+     totalNtwin                                                                                     !< total number of active twin systems
+   integer(pInt),               allocatable, dimension(:) :: &
+     Nslip, &                                                                                       !< number of active slip systems for each family
+     Ntwin                                                                                          !< number of active twin systems for each family
+   integer(kind(undefined_ID)), allocatable, dimension(:) :: &
+     outputID                                                                                       !< ID of each post result output
+ end type                                                                                           !< container type for internal constitutive parameters
 
  type(tParameters), dimension(:), allocatable, private :: param                                     !< containers of constitutive parameters (len Ninstance)
 
  type, private :: tPhenopowerlawState
-   real(pReal), pointer,     dimension(:,:) :: &
+   real(pReal), pointer, dimension(:) :: &
+     sumGamma, &                                                                                    ! ToDo: why not make a dependent state?
+     sumF                                                                                           ! ToDo: why not make a dependent state?
+   real(pReal), pointer, dimension(:,:) :: &
      xi_slip, &
      xi_twin, &
      gamma_slip, &
      gamma_twin, &
      whole
-   real(pReal), pointer,     dimension(:) :: &
-     sumGamma, &
-     sumF
  end type
 
  type(tPhenopowerlawState), allocatable, dimension(:), private :: &
@@ -115,6 +115,7 @@ subroutine plastic_phenopowerlaw_init
    compiler_options
 #endif
  use prec, only: &
+   pStringLen, &
    dEq0
  use debug, only: &
    debug_level, &
@@ -130,6 +131,7 @@ subroutine plastic_phenopowerlaw_init
    phase_plasticity, &
    phase_plasticityInstance, &
    phase_Noutput, &
+   material_allocatePlasticState, &
    PLASTICITY_PHENOPOWERLAW_LABEL, &
    PLASTICITY_PHENOPOWERLAW_ID, &
    material_phase, &
@@ -138,20 +140,17 @@ subroutine plastic_phenopowerlaw_init
    MATERIAL_partPhase, &
    config_phase
  use lattice
- use numerics,only: &
-   numerics_integrator
 
  implicit none
-
  integer(pInt) :: &
-   maxNinstance, &
+   Ninstance, &
-   instance,p,j,k, o, i,&
+   p, i, &
    NipcMyPhase, outputSize, &
-   sizeState,sizeDotState, &
+   sizeState, sizeDotState, &
    startIndex, endIndex
 
- integer(pInt),          dimension(0), parameter :: emptyIntArray = [integer(pInt)::]
+ integer(pInt),          dimension(0), parameter :: emptyIntArray    = [integer(pInt)::]
- real(pReal),            dimension(0), parameter :: emptyRealArray = [real(pReal)::]
+ real(pReal),            dimension(0), parameter :: emptyRealArray   = [real(pReal)::]
  character(len=65536),   dimension(0), parameter :: emptyStringArray = [character(len=65536)::]
 
  type(tParameters) :: &
@@ -163,152 +162,151 @@ subroutine plastic_phenopowerlaw_init
  integer(kind(undefined_ID)) :: &
    outputID                                                                                         !< ID of each post result output
 
- character(len=512) :: &
+ character(len=pStringLen) :: &
-   extmsg    = '', &
+   structure = '',&
-   structure = ''
+   extmsg = ''
- character(len=65536), dimension(:), allocatable :: outputs
+ character(len=65536), dimension(:), allocatable :: &
+  outputs
 
  write(6,'(/,a)')   ' <<<+-  constitutive_'//PLASTICITY_PHENOPOWERLAW_label//' init  -+>>>'
  write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
 #include "compilation_info.f90"
 
- maxNinstance = int(count(phase_plasticity == PLASTICITY_PHENOPOWERLAW_ID),pInt)
+ Ninstance = int(count(phase_plasticity == PLASTICITY_PHENOPOWERLAW_ID),pInt)
  if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
-   write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
+   write(6,'(a16,1x,i5,/)') '# instances:',Ninstance
 
- allocate(plastic_phenopowerlaw_sizePostResult(maxval(phase_Noutput),maxNinstance),source=0_pInt)
+ allocate(plastic_phenopowerlaw_sizePostResult(maxval(phase_Noutput),Ninstance),source=0_pInt)
- allocate(plastic_phenopowerlaw_output(maxval(phase_Noutput),maxNinstance))
+ allocate(plastic_phenopowerlaw_output(maxval(phase_Noutput),Ninstance))
           plastic_phenopowerlaw_output = ''
 
- allocate(param(maxNinstance))                                                                      ! one container of parameters per instance
+ allocate(param(Ninstance))
- allocate(state(maxNinstance))
+ allocate(state(Ninstance))
- allocate(dotState(maxNinstance))
+ allocate(dotState(Ninstance))
 
  do p = 1_pInt, size(phase_plasticityInstance)
    if (phase_plasticity(p) /= PLASTICITY_PHENOPOWERLAW_ID) cycle
-   instance = phase_plasticityInstance(p)
+   associate(prm => param(phase_plasticityInstance(p)), &
-   associate(prm => param(instance),stt => state(instance),dot => dotState(instance))
+             dot => dotState(phase_plasticityInstance(p)), &
-   extmsg = ''
+             stt => state(phase_plasticityInstance(p)))
 
-   structure      = config_phase(p)%getString('lattice_structure')
+   structure          = config_phase(p)%getString('lattice_structure')
 
+!--------------------------------------------------------------------------------------------------
+!  optional parameters that need to be defined
+   prm%twinB          = config_phase(p)%getFloat('twin_b',defaultVal=1.0_pReal)
+   prm%twinC          = config_phase(p)%getFloat('twin_c',defaultVal=0.0_pReal)
+   prm%twinD          = config_phase(p)%getFloat('twin_d',defaultVal=0.0_pReal)
+   prm%twinE          = config_phase(p)%getFloat('twin_e',defaultVal=0.0_pReal)
+
+   prm%aTolResistance = config_phase(p)%getFloat('atol_resistance',defaultVal=1.0_pReal)
+   prm%aTolShear      = config_phase(p)%getFloat('atol_shear',     defaultVal=1.0e-6_pReal)
+   prm%aTolTwinfrac   = config_phase(p)%getFloat('atol_twinfrac',  defaultVal=1.0e-6_pReal)
+   
+   ! sanity checks
+   if (prm%aTolResistance <= 0.0_pReal) extmsg = trim(extmsg)//'aTolresistance '
+   if (prm%aTolShear      <= 0.0_pReal) extmsg = trim(extmsg)//'aTolShear '
+   if (prm%aTolTwinfrac   <= 0.0_pReal) extmsg = trim(extmsg)//'atoltwinfrac '
+
+!--------------------------------------------------------------------------------------------------
+! slip related parameters
    prm%Nslip      = config_phase(p)%getInts('nslip',defaultVal=emptyIntArray)
    prm%totalNslip = sum(prm%Nslip)
-   if (size(prm%Nslip) > count(lattice_NslipSystem(:,p) > 0_pInt)) &
-     call IO_error(150_pInt,ext_msg='Nslip')
-   if (any(lattice_NslipSystem(1:size(prm%Nslip),p)-prm%Nslip < 0_pInt)) &
-     call IO_error(150_pInt,ext_msg='Nslip')
-
    slipActive: if (prm%totalNslip > 0_pInt) then
-
      prm%Schmid_slip          = lattice_SchmidMatrix_slip(prm%Nslip,structure(1:3),&
                                                           config_phase(p)%getFloat('c/a',defaultVal=0.0_pReal))
-     ! reading in slip related parameters
+     if(structure=='bcc') then
+       prm%nonSchmidCoeff     = config_phase(p)%getFloats('nonschmid_coefficients',&
+                                                          defaultVal = emptyRealArray)
+       prm%nonSchmid_pos      = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,+1_pInt)
+       prm%nonSchmid_neg      = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,-1_pInt)
+     else
+       prm%nonSchmid_pos      = prm%Schmid_slip
+       prm%nonSchmid_neg      = prm%Schmid_slip
+     endif
+     prm%interaction_SlipSlip = lattice_interaction_SlipSlip(prm%Nslip, &
+                                                             config_phase(p)%getFloats('interaction_slipslip'), &
+                                                             structure(1:3))
+
      prm%xi_slip_0            = config_phase(p)%getFloats('tau0_slip',   requiredShape=shape(prm%Nslip))
      prm%xi_slip_sat          = config_phase(p)%getFloats('tausat_slip', requiredShape=shape(prm%Nslip))
-     prm%interaction_SlipSlip = lattice_interaction_SlipSlip(prm%Nslip,config_phase(p)%getFloats('interaction_slipslip'), &
-                                                                      structure(1:3))
      prm%H_int                = config_phase(p)%getFloats('h_int',       requiredShape=shape(prm%Nslip), &
                                                                          defaultVal=[(0.0_pReal,i=1_pInt,size(prm%Nslip))])
-     prm%nonSchmidCoeff       = config_phase(p)%getFloats('nonschmid_coefficients',&
+ 
-                                                                         defaultVal = emptyRealArray )
-     if(structure=='bcc') then
-      prm%nonSchmid_pos        = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,+1_pInt)
-      prm%nonSchmid_neg        = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,-1_pInt)
-     else
-      prm%nonSchmid_pos        = prm%Schmid_slip
-      prm%nonSchmid_neg        = prm%Schmid_slip
-     endif
      prm%gdot0_slip           = config_phase(p)%getFloat('gdot0_slip')
      prm%n_slip               = config_phase(p)%getFloat('n_slip')
      prm%a_slip               = config_phase(p)%getFloat('a_slip')
      prm%h0_SlipSlip          = config_phase(p)%getFloat('h0_slipslip')
 
-     ! sanity checks for slip related parameters
+     ! expand: family => system
-     if (any(prm%xi_slip_0 < 0.0_pReal .and. prm%Nslip > 0_pInt)) &
+     prm%xi_slip_0   = math_expand(prm%xi_slip_0,  prm%Nslip)
-       extmsg = trim(extmsg)//"xi_slip_0 "
-     if (any(prm%xi_slip_sat < prm%xi_slip_0 .and. prm%Nslip > 0_pInt)) &
-       extmsg = trim(extmsg)//"xi_slip_sat "
-
-     if (prm%gdot0_slip <= 0.0_pReal) extmsg = trim(extmsg)//"gdot0_slip "
-     if (dEq0(prm%a_slip))            extmsg = trim(extmsg)//"a_slip "                              ! ToDo: negative values ok?
-     if (dEq0(prm%n_slip))            extmsg = trim(extmsg)//"n_slip "                              ! ToDo: negative values ok?
-
-     ! expand slip related parameters from system => family
-     prm%xi_slip_0   = math_expand(prm%xi_slip_0,prm%Nslip)
      prm%xi_slip_sat = math_expand(prm%xi_slip_sat,prm%Nslip)
-     prm%H_int       = math_expand(prm%H_int,prm%Nslip)
+     prm%H_int       = math_expand(prm%H_int,      prm%Nslip)
+
+     ! sanity checks
+     if (prm%gdot0_slip <= 0.0_pReal)           extmsg = trim(extmsg)//'gdot0_slip '
+     if (dEq0(prm%a_slip))                      extmsg = trim(extmsg)//'a_slip '                              ! ToDo: negative values ok?
+     if (dEq0(prm%n_slip))                      extmsg = trim(extmsg)//'n_slip '                              ! ToDo: negative values ok?
+     if (any(prm%xi_slip_0   <= 0.0_pReal))     extmsg = trim(extmsg)//'xi_slip_0 '
+     if (any(prm%xi_slip_sat <  prm%xi_slip_0)) extmsg = trim(extmsg)//'xi_slip_sat '
    else slipActive
      allocate(prm%interaction_SlipSlip(0,0))
      allocate(prm%xi_slip_0(0))
    endif slipActive
 
+!--------------------------------------------------------------------------------------------------
+! twin related parameters
    prm%Ntwin      = config_phase(p)%getInts('ntwin', defaultVal=emptyIntArray)
    prm%totalNtwin = sum(prm%Ntwin)
-   if (size(prm%Ntwin) > count(lattice_NtwinSystem(:,p) > 0_pInt)) &
-     call IO_error(150_pInt,ext_msg='Ntwin')
-   if (any(lattice_NtwinSystem(1:size(prm%Ntwin),p)-prm%Ntwin < 0_pInt)) &
-     call IO_error(150_pInt,ext_msg='Ntwin')
-
    twinActive: if (prm%totalNtwin > 0_pInt) then
      prm%Schmid_twin          = lattice_SchmidMatrix_twin(prm%Ntwin,structure(1:3),&
                                                           config_phase(p)%getFloat('c/a',defaultVal=0.0_pReal))
-     ! reading in twin related parameters
+     prm%interaction_TwinTwin = lattice_interaction_TwinTwin(prm%Ntwin,&
+                                                             config_phase(p)%getFloats('interaction_twintwin'), &
+                                                             structure(1:3))
+     prm%gamma_twin_char      = lattice_characteristicShear_twin(prm%Ntwin,structure(1:3),&
+                                                            config_phase(p)%getFloat('c/a'))
+
      prm%xi_twin_0            = config_phase(p)%getFloats('tau0_twin',requiredShape=shape(prm%Ntwin))
-     prm%interaction_TwinTwin = lattice_interaction_TwinTwin(prm%Ntwin,config_phase(p)%getFloats('interaction_twintwin'), &
-                                                                      structure(1:3))
 
-     prm%gdot0_twin  = config_phase(p)%getFloat('gdot0_twin')
+     prm%gdot0_twin           = config_phase(p)%getFloat('gdot0_twin')
-     prm%n_twin      = config_phase(p)%getFloat('n_twin')
+     prm%n_twin               = config_phase(p)%getFloat('n_twin')
-     prm%spr         = config_phase(p)%getFloat('s_pr')
+     prm%spr                  = config_phase(p)%getFloat('s_pr')
-     prm%h0_TwinTwin = config_phase(p)%getFloat('h0_twintwin')
+     prm%h0_TwinTwin          = config_phase(p)%getFloat('h0_twintwin')
 
-     ! sanity checks for twin related parameters
+     ! expand: family => system
-     if (any(prm%xi_twin_0 < 0.0_pReal .and. prm%Ntwin > 0_pInt)) &
+     prm%xi_twin_0   = math_expand(prm%xi_twin_0,  prm%Ntwin)
-       extmsg = trim(extmsg)//"xi_twin_0 "
-     if (prm%gdot0_twin <= 0.0_pReal) extmsg = trim(extmsg)//"gdot0_twin "
-     if (dEq0(prm%n_twin))            extmsg = trim(extmsg)//"n_twin "                              ! ToDo: negative values ok?
 
-     ! expand slip related parameters from system => family
+     ! sanity checks
-     prm%xi_twin_0   = math_expand(prm%xi_twin_0,prm%Ntwin)
+     if (prm%gdot0_twin <= 0.0_pReal)  extmsg = trim(extmsg)//'gdot0_twin '
+     if (dEq0(prm%n_twin))             extmsg = trim(extmsg)//'n_twin '                                      ! ToDo: negative values ok?
    else twinActive
      allocate(prm%interaction_TwinTwin(0,0))
      allocate(prm%xi_twin_0(0))
    endif twinActive
 
-   prm%gamma_twin_char = lattice_characteristicShear_twin(prm%Ntwin,structure(1:3),&
+!--------------------------------------------------------------------------------------------------
-   config_phase(p)%getFloat('c/a',defaultVal=0.0_pReal))
+! slip-twin related parameters
-
    slipAndTwinActive: if (prm%totalNslip > 0_pInt .and. prm%totalNtwin > 0_pInt) then
      prm%interaction_SlipTwin = lattice_interaction_SlipTwin(prm%Nslip,prm%Ntwin,&
-                                       config_phase(p)%getFloats('interaction_sliptwin'), &
+                                                             config_phase(p)%getFloats('interaction_sliptwin'), &
-                                                                      structure(1:3))
+                                                             structure(1:3))
      prm%interaction_TwinSlip = lattice_interaction_TwinSlip(prm%Ntwin,prm%Nslip,&
-                                       config_phase(p)%getFloats('interaction_twinslip'), &
+                                                             config_phase(p)%getFloats('interaction_twinslip'), &
-                                                                      structure(1:3))
+                                                             structure(1:3))
    else slipAndTwinActive
-     allocate(prm%interaction_SlipTwin(prm%totalNslip,prm%TotalNtwin))                              ! at least one dimension 0
+     allocate(prm%interaction_SlipTwin(prm%totalNslip,prm%TotalNtwin))                                        ! at least one dimension is 0
-     allocate(prm%interaction_TwinSlip(prm%totalNtwin,prm%TotalNslip))                              ! at least one dimension 0
+     allocate(prm%interaction_TwinSlip(prm%totalNtwin,prm%TotalNslip))                                        ! at least one dimension is 0
      prm%h0_TwinSlip = 0.0_pReal
    endif slipAndTwinActive
 
-   ! optional parameters that should be defined
+!--------------------------------------------------------------------------------------------------
-   prm%twinB       = config_phase(p)%getFloat('twin_b',defaultVal=1.0_pReal)
+!  exit if any parameter is out of range
-   prm%twinC       = config_phase(p)%getFloat('twin_c',defaultVal=0.0_pReal)
+   if (extmsg /= '') &
-   prm%twinD       = config_phase(p)%getFloat('twin_d',defaultVal=0.0_pReal)
+     call IO_error(211_pInt,ext_msg=trim(extmsg)//'('//PLASTICITY_PHENOPOWERLAW_label//')')
-   prm%twinE       = config_phase(p)%getFloat('twin_e',defaultVal=0.0_pReal)
-
-   prm%aTolResistance = config_phase(p)%getFloat('atol_resistance',defaultVal=1.0_pReal)
-   prm%aTolShear      = config_phase(p)%getFloat('atol_shear',     defaultVal=1.0e-6_pReal)
-   prm%aTolTwinfrac   = config_phase(p)%getFloat('atol_twinfrac',  defaultVal=1.0e-6_pReal)
-
-   if (prm%aTolResistance <= 0.0_pReal) extmsg = trim(extmsg)//"aTolresistance "
-   if (prm%aTolShear      <= 0.0_pReal) extmsg = trim(extmsg)//"aTolShear "
-   if (prm%aTolTwinfrac   <= 0.0_pReal) extmsg = trim(extmsg)//"atoltwinfrac "
-
-   if (extmsg /= '') call IO_error(211_pInt,ip=instance,&
-                                 ext_msg=trim(extmsg)//'('//PLASTICITY_PHENOPOWERLAW_label//')')
 
+!--------------------------------------------------------------------------------------------------
+!  output pararameters
    outputs = config_phase(p)%getStrings('(output)',defaultVal=emptyStringArray)
    allocate(prm%outputID(0))
    do i=1_pInt, size(outputs)
@@ -349,8 +347,8 @@ subroutine plastic_phenopowerlaw_init
       end select
 
       if (outputID /= undefined_ID) then
-        plastic_phenopowerlaw_output(i,instance) = outputs(i)
+        plastic_phenopowerlaw_output(i,phase_plasticityInstance(p)) = outputs(i)
-        plastic_phenopowerlaw_sizePostResult(i,instance) = outputSize
+        plastic_phenopowerlaw_sizePostResult(i,phase_plasticityInstance(p)) = outputSize
         prm%outputID = [prm%outputID , outputID]
       endif
 
@@ -361,32 +359,12 @@ subroutine plastic_phenopowerlaw_init
    NipcMyPhase = count(material_phase == p)                                                         ! number of IPCs containing my phase
    sizeState = size(['tau_slip  ','gamma_slip']) * prm%TotalNslip &
              + size(['tau_twin  ','gamma_twin']) * prm%TotalNtwin &
-             + size(['sum(gamma)','sum(f)    '])
+             + size(['sum(gamma)','sum(f)    '])                                                    ! ToDo: only needed if either twin or slip active!
-
-!--------------------------------------------------------------------------------------------------
-! ToDo: This could be done by a function (in constitutive?)
    sizeDotState = sizeState
-   plasticState(p)%sizeState = sizeState
-   plasticState(p)%sizeDotState = sizeDotState
-   plasticState(p)%sizePostResults = sum(plastic_phenopowerlaw_sizePostResult(:,instance))
-   plasticState(p)%nSlip = prm%totalNslip
-   plasticState(p)%nTwin = prm%totalNtwin
-   allocate(plasticState(p)%aTolState          (   sizeState),             source=0.0_pReal)
-   allocate(plasticState(p)%state0             (   sizeState,NipcMyPhase), source=0.0_pReal)
-   allocate(plasticState(p)%partionedState0    (   sizeState,NipcMyPhase), source=0.0_pReal)
-   allocate(plasticState(p)%subState0          (   sizeState,NipcMyPhase), source=0.0_pReal)
-   allocate(plasticState(p)%state              (   sizeState,NipcMyPhase), source=0.0_pReal)
 
-   allocate(plasticState(p)%dotState           (sizeDotState,NipcMyPhase), source=0.0_pReal)
+   call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0_pInt, &
-   allocate(plasticState(p)%deltaState               (0_pInt,NipcMyPhase), source=0.0_pReal)
+                                      prm%totalNslip,prm%totalNtwin,0_pInt)
-   if (any(numerics_integrator == 1_pInt)) then
+   plasticState(p)%sizePostResults = sum(plastic_phenopowerlaw_sizePostResult(:,phase_plasticityInstance(p)))
-     allocate(plasticState(p)%previousDotState (sizeDotState,NipcMyPhase),source=0.0_pReal)
-     allocate(plasticState(p)%previousDotState2(sizeDotState,NipcMyPhase),source=0.0_pReal)
-   endif
-   if (any(numerics_integrator == 4_pInt)) &
-     allocate(plasticState(p)%RK4dotState      (sizeDotState,NipcMyPhase), source=0.0_pReal)
-   if (any(numerics_integrator == 5_pInt)) &
-     allocate(plasticState(p)%RKCK45dotState (6,sizeDotState,NipcMyPhase), source=0.0_pReal)
 
 
 !--------------------------------------------------------------------------------------------------
@@ -432,7 +410,7 @@ subroutine plastic_phenopowerlaw_init
    dot%gamma_twin => plasticState(p)%dotState(startIndex:endIndex,:)
    plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolShear
 
-   plasticState(p)%state0 = plasticState(p)%state
+   plasticState(p)%state0 = plasticState(p)%state                                                 ! ToDo: this could be done centrally
    dot%whole => plasticState(p)%dotState
 
    end associate
@@ -544,25 +522,21 @@ subroutine plastic_phenopowerlaw_dotState(Mp,instance,of)
  dot%gamma_slip(:,of) = abs(gdot_slip_pos+gdot_slip_neg)
  dot%sumGamma(of) = sum(dot%gamma_slip(:,of))
  call kinetics_twin(prm,stt,of,Mp,dot%gamma_twin(:,of))
- if (stt%sumF(of) < 0.98_pReal) dot%sumF(of) = sum(dot%gamma_twin(:,of)/prm%gamma_twin_char)
+ if (prm%totalNtwin > 0_pInt) dot%sumF(of) = merge(sum(dot%gamma_twin(:,of)/prm%gamma_twin_char), &
+                                                   0.0_pReal, &
+                                                   stt%sumF(of) < 0.98_pReal)
 
 !--------------------------------------------------------------------------------------------------
 ! hardening
  hardeningSlip: do i = 1_pInt, prm%totalNslip
-   dot%xi_slip(i,of) = &
+   dot%xi_slip(i,of) = dot_product(prm%interaction_SlipSlip(i,:),right_SlipSlip*dot%gamma_slip(:,of)) &
-     c_SlipSlip * left_SlipSlip(i) &
+                           * c_SlipSlip * left_SlipSlip(i) &
-     * dot_product(prm%interaction_SlipSlip(i,:),right_SlipSlip*dot%gamma_slip(:,of)) &
+                     + dot_product(prm%interaction_SlipTwin(i,:),dot%gamma_twin(:,of))
-                    + &
-       dot_product(prm%interaction_SlipTwin(i,:),dot%gamma_twin(:,of))
  enddo hardeningSlip
 
  hardeningTwin: do i = 1_pInt, prm%totalNtwin
-   dot%xi_twin(i,of) = &
+   dot%xi_twin(i,of) = c_TwinSlip * dot_product(prm%interaction_TwinSlip(i,:),dot%gamma_slip(:,of)) &
-     c_TwinSlip &
+                     + c_TwinTwin * dot_product(prm%interaction_TwinTwin(i,:),dot%gamma_twin(:,of))
-     * dot_product(prm%interaction_TwinSlip(i,:),dot%gamma_slip(:,of)) &
-                    + &
-     c_TwinTwin &
-     * dot_product(prm%interaction_TwinTwin(i,:),dot%gamma_twin(:,of))
  enddo hardeningTwin
 
  end associate
@@ -572,7 +546,7 @@ end subroutine plastic_phenopowerlaw_dotState
 
 !--------------------------------------------------------------------------------------------------
 !> @brief calculates shear rates on slip systems and derivatives with respect to resolved stress
-!> @details: Shear rates are calculated only optionally. NOTE: Agains the common convention, the
+!> @details: Shear rates are calculated only optionally. NOTE: Against the common convention, the
 !> result (i.e. intent(out)) variables are the last to have the optional arguments at the end
 !--------------------------------------------------------------------------------------------------
 pure subroutine kinetics_slip(prm,stt,of,Mp,gdot_slip_pos,gdot_slip_neg, &
@@ -602,26 +576,39 @@ pure subroutine kinetics_slip(prm,stt,of,Mp,gdot_slip_pos,gdot_slip_neg, &
    tau_slip_pos, &
    tau_slip_neg
  integer(pInt) :: i
+ logical       :: nonSchmidActive
+
+ nonSchmidActive = size(prm%nonSchmidCoeff) > 0_pInt
 
  do i = 1_pInt, prm%totalNslip
-   tau_slip_pos(i) = math_mul33xx33(Mp,prm%nonSchmid_pos(1:3,1:3,i))
+   tau_slip_pos(i) =       math_mul33xx33(Mp,prm%nonSchmid_pos(1:3,1:3,i))
-   tau_slip_neg(i) = math_mul33xx33(Mp,prm%nonSchmid_neg(1:3,1:3,i))
+   tau_slip_neg(i) = merge(math_mul33xx33(Mp,prm%nonSchmid_neg(1:3,1:3,i)), &
+                           0.0_pReal, nonSchmidActive)
  enddo
 
- gdot_slip_pos = 0.5_pReal*prm%gdot0_slip &
+ where(dNeq0(tau_slip_pos))
-               * sign(abs(tau_slip_pos/stt%xi_slip(:,of))**prm%n_slip,  tau_slip_pos)
+   gdot_slip_pos = prm%gdot0_slip * merge(0.5_pReal,1.0_pReal, nonSchmidActive) &                   ! 1/2 if non-Schmid active
- gdot_slip_neg = 0.5_pReal*prm%gdot0_slip &
+                 * sign(abs(tau_slip_pos/stt%xi_slip(:,of))**prm%n_slip,  tau_slip_pos)
-               * sign(abs(tau_slip_neg/stt%xi_slip(:,of))**prm%n_slip,  tau_slip_neg)
+ else where
+   gdot_slip_pos = 0.0_pReal
+ end where
+
+ where(dNeq0(tau_slip_neg))
+   gdot_slip_neg = 0.5_pReal*prm%gdot0_slip &
+                 * sign(abs(tau_slip_neg/stt%xi_slip(:,of))**prm%n_slip,  tau_slip_neg)
+ else where
+   gdot_slip_neg = 0.0_pReal
+ end where
 
  if (present(dgdot_dtau_slip_pos)) then
-   where(dNeq0(tau_slip_pos))
+   where(dNeq0(gdot_slip_pos))
      dgdot_dtau_slip_pos = gdot_slip_pos*prm%n_slip/tau_slip_pos
    else where
      dgdot_dtau_slip_pos = 0.0_pReal
    end where
  endif
  if (present(dgdot_dtau_slip_neg)) then
-   where(dNeq0(tau_slip_neg))
+   where(dNeq0(gdot_slip_neg))
      dgdot_dtau_slip_neg = gdot_slip_neg*prm%n_slip/tau_slip_neg
    else where
      dgdot_dtau_slip_neg = 0.0_pReal
@@ -633,7 +620,7 @@ end subroutine kinetics_slip
 
 !--------------------------------------------------------------------------------------------------
 !> @brief calculates shear rates on twin systems and derivatives with respect to resolved stress
-!> @details: Shear rates are calculated only optionally. NOTE: Agains the common convention, the
+!> @details: Shear rates are calculated only optionally. NOTE: Against the common convention, the
 !> result (i.e. intent(out)) variables are the last to have the optional arguments at the end
 !--------------------------------------------------------------------------------------------------
 pure subroutine kinetics_twin(prm,stt,of,Mp,gdot_twin,dgdot_dtau_twin)
@@ -663,11 +650,15 @@ pure subroutine kinetics_twin(prm,stt,of,Mp,gdot_twin,dgdot_dtau_twin)
  do i = 1_pInt, prm%totalNtwin
    tau_twin(i)  = math_mul33xx33(Mp,prm%Schmid_twin(1:3,1:3,i))
  enddo
- gdot_twin = merge((1.0_pReal-stt%sumF(of))*prm%gdot0_twin*(abs(tau_twin)/stt%xi_twin(:,of))**prm%n_twin, &
+ 
-                    0.0_pReal, tau_twin>0.0_pReal)
+ where(tau_twin > 0.0_pReal)
+   gdot_twin = (1.0_pReal-stt%sumF(of))*prm%gdot0_twin*(abs(tau_twin)/stt%xi_twin(:,of))**prm%n_twin
+ else where            
+   gdot_twin = 0.0_pReal
+ end where
 
  if (present(dgdot_dtau_twin)) then
-   where(dNeq0(tau_twin))
+   where(dNeq0(gdot_twin))
      dgdot_dtau_twin = gdot_twin*prm%n_twin/tau_twin
    else where
      dgdot_dtau_twin = 0.0_pReal
@@ -681,14 +672,8 @@ end subroutine kinetics_twin
 !> @brief return array of constitutive results
 !--------------------------------------------------------------------------------------------------
 function plastic_phenopowerlaw_postResults(Mp,instance,of) result(postResults)
- use material, only: &
-   material_phase, &
-   plasticState, &
-   phasememberAt, &
-   phase_plasticityInstance
  use math, only: &
-   math_mul33xx33, &
+   math_mul33xx33
-   math_Mandel6to33
 
  implicit none
  real(pReal), dimension(3,3), intent(in) :: &
@@ -701,7 +686,7 @@ function plastic_phenopowerlaw_postResults(Mp,instance,of) result(postResults)
    postResults
 
  integer(pInt) :: &
-   o,c,i,j
+   o,c,i
  real(pReal), dimension(param(instance)%totalNslip) :: &
    gdot_slip_pos,gdot_slip_neg