Merge branch 'development' into 19-NewStylePhenopowerlaw

.gitlab-ci.yml

@@ -158,12 +158,12 @@ Post_AverageDown:
     - master
     - release
 
-Post_General:
+#Post_General:
-  stage: postprocessing
+#  stage: postprocessing
-  script: PostProcessing/test.py
+#  script: PostProcessing/test.py
-  except:
+#  except:
-    - master
+#    - master
-    - release
+#    - release
 
 Post_GeometryReconstruction:
   stage: postprocessing
@@ -364,12 +364,12 @@ Phenopowerlaw_singleSlip:
     - master
     - release
 
-TextureComponents:
+#TextureComponents:
-  stage: spectral
+#  stage: spectral
-  script: TextureComponents/test.py
+#  script: TextureComponents/test.py
-  except:
+#  except:
-    - master
+#    - master
-    - release
+#    - release
 
 
 ###################################################################################################

DAMASK_prerequisites.sh

@@ -119,6 +119,9 @@ for executable in mpirun mpiexec; do
   getDetails $executable '--version'
 done
 
+firstLevel "CMake"
+getDetails cmake --version
+
 firstLevel "Abaqus"
 cd installation/mods_Abaqus                                                                         # to have the right environment file
 for executable in abaqus abq2017 abq2018; do

VERSION

@@ -1 +1 @@
-v2.0.2-983-g0a648459
+v2.0.2-1124-g0ff9e9fa

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

lib/damask/orientation.py

@@ -33,7 +33,7 @@ class Quaternion:
     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
+                  when viewing from the end point of the rotation axis 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π]
@@ -48,206 +48,138 @@ class Quaternion:
     """
 
     def __init__(self,
-                 quatArray = [1.0,0.0,0.0,0.0]):
+                 quat = None,
-      """Initializes to identity if not given"""
+                 q    =  1.0,
-      self.w, \
+                 p    = np.zeros(3,dtype=float)):
-      self.x, \
+      """Initializes to identity unless specified"""
-      self.y, \
+      self.q = quat[0] if quat is not None else q
-      self.z = quatArray
+      self.p = np.array(quat[1:4]) if quat is not None else p
       self.homomorph()
 
     def __iter__(self):
       """Components"""
-      return iter([self.w,self.x,self.y,self.z])
+      return iter(self.asList())
 
     def __copy__(self):
-      """Create copy"""
+      """Copy"""
-      Q = Quaternion([self.w,self.x,self.y,self.z])
+      return self.__class__(q=self.q,p=self.p.copy())
-      return Q
 
     copy = __copy__
 
     def __repr__(self):
-      """Readbable string"""
+      """Readable string"""
-      return 'Quaternion(real=%+.6f, imag=<%+.6f, %+.6f, %+.6f>)' % \
+      return 'Quaternion(real={q:+.6f}, imag=<{p[0]:+.6f}, {p[1]:+.6f}, {p[2]:+.6f}>)'.format(q=self.q,p=self.p)
-          (self.w, self.x, self.y, self.z)
 
     def __pow__(self, exponent):
       """Power"""
-      omega = math.acos(self.w)
+      omega = math.acos(self.q)
-      vRescale = math.sin(exponent*omega)/math.sin(omega)
+      return self.__class__(q=         math.cos(exponent*omega),
-      Q = Quaternion()
+                            p=self.p * math.sin(exponent*omega)/math.sin(omega))
-      Q.w = math.cos(exponent*omega)
-      Q.x = self.x * vRescale
-      Q.y = self.y * vRescale
-      Q.z = self.z * vRescale
-      return Q
 
     def __ipow__(self, exponent):
       """In-place power"""
-      omega    = math.acos(self.w)
+      omega = math.acos(self.q)
-      vRescale = math.sin(exponent*omega)/math.sin(omega)
+      self.q  = math.cos(exponent*omega)
-      self.w  = np.cos(exponent*omega)
+      self.p *= math.sin(exponent*omega)/math.sin(omega)
-      self.x *= vRescale
-      self.y *= vRescale
-      self.z *= vRescale
       return self
 
     def __mul__(self, other):
       """Multiplication"""
+      # Rowenhorst_etal2015 MSMSE: value of P is selected as -1
+      P = -1.0
       try:                                                          # quaternion
-          Aw = self.w
+          return self.__class__(q=self.q*other.q - np.dot(self.p,other.p),
-          Ax = self.x
+                                p=self.q*other.p + other.q*self.p + P * np.cross(self.p,other.p))
-          Ay = self.y
-          Az = self.z
-          Bw = other.w
-          Bx = other.x
-          By = other.y
-          Bz = other.z
-          Q = Quaternion()
-          Q.w = - Ax * Bx - Ay * By - Az * Bz + Aw * Bw
-          Q.x = + Ax * Bw + Ay * Bz - Az * By + Aw * Bx
-          Q.y = - Ax * Bz + Ay * Bw + Az * Bx + Aw * By
-          Q.z = + Ax * By - Ay * Bx + Az * Bw + Aw * Bz
-          return Q
       except: pass
-      try:                                                         # vector (perform active rotation, i.e. q*v*q.conjugated)
+      try:                                                         # vector (perform passive rotation)
-          w = self.w
+          ( x, y, z) = self.p
-          x = self.x
+          (Vx,Vy,Vz) = other[0:3]
-          y = self.y
+          A = self.q*self.q - np.dot(self.p,self.p)
-          z = self.z
+          B = 2.0 * (x*Vx + y*Vy + z*Vz)
-          Vx = other[0]
+          C = 2.0 * P*self.q
-          Vy = other[1]
-          Vz = other[2]
 
-          return np.array([\
+          return np.array([
-             w * w * Vx + 2 * y * w * Vz - 2 * z * w * Vy + \
+            A*Vx + B*x + C*(y*Vz - z*Vy),
-             x * x * Vx + 2 * y * x * Vy + 2 * z * x * Vz - \
+            A*Vy + B*y + C*(z*Vx - x*Vz),
-             z * z * Vx - y * y * Vx,
+            A*Vz + B*z + C*(x*Vy - y*Vx),
-             2 * x * y * Vx + y * y * Vy + 2 * z * y * Vz + \
+            ])
-             2 * w * z * Vx - z * z * Vy + w * w * Vy - \
-             2 * x * w * Vz - x * x * Vy,
-             2 * x * z * Vx + 2 * y * z * Vy + \
-             z * z * Vz - 2 * w * y * Vx - y * y * Vz + \
-             2 * w * x * Vy - x * x * Vz + w * w * Vz ])
       except: pass
       try:                                                        # scalar
-          Q = self.copy()
+          return self.__class__(q=self.q*other,
-          Q.w *= other
+                                p=self.p*other)
-          Q.x *= other
-          Q.y *= other
-          Q.z *= other
-          return Q
       except:
           return self.copy()
 
     def __imul__(self, other):
       """In-place multiplication"""
+      # Rowenhorst_etal2015 MSMSE: value of P is selected as -1
+      P = -1.0
       try:                                                        # Quaternion
-          Aw = self.w
+          self.q = self.q*other.q - np.dot(self.p,other.p)
-          Ax = self.x
+          self.p = self.q*other.p + other.q*self.p + P * np.cross(self.p,other.p)
-          Ay = self.y
-          Az = self.z
-          Bw = other.w
-          Bx = other.x
-          By = other.y
-          Bz = other.z
-          self.w = - Ax * Bx - Ay * By - Az * Bz + Aw * Bw
-          self.x = + Ax * Bw + Ay * Bz - Az * By + Aw * Bx
-          self.y = - Ax * Bz + Ay * Bw + Az * Bx + Aw * By
-          self.z = + Ax * By - Ay * Bx + Az * Bw + Aw * Bz
       except: pass
       return self
 
     def __div__(self, other):
       """Division"""
       if isinstance(other, (int,float)):
-        w = self.w / other
+        return self.__class__(q=self.q / other,
-        x = self.x / other
+                              p=self.p / other)
-        y = self.y / other
-        z = self.z / other
-        return self.__class__([w,x,y,z])
       else:
           return NotImplemented
 
     def __idiv__(self, other):
       """In-place division"""
       if isinstance(other, (int,float)):
-          self.w /= other
+          self.q /= other
-          self.x /= other
+          self.p /= other
-          self.y /= other
-          self.z /= other
       return self
 
     def __add__(self, other):
       """Addition"""
       if isinstance(other, Quaternion):
-        w = self.w + other.w
+        return self.__class__(q=self.q + other.q,
-        x = self.x + other.x
+                              p=self.p + other.p)
-        y = self.y + other.y
-        z = self.z + other.z
-        return self.__class__([w,x,y,z])
       else:
           return NotImplemented
 
     def __iadd__(self, other):
       """In-place addition"""
       if isinstance(other, Quaternion):
-          self.w += other.w
+          self.q += other.q
-          self.x += other.x
+          self.p += other.p
-          self.y += other.y
-          self.z += other.z
       return self
 
     def __sub__(self, other):
       """Subtraction"""
       if isinstance(other, Quaternion):
-          Q = self.copy()
+          return self.__class__(q=self.q - other.q,
-          Q.w -= other.w
+                                p=self.p - other.p)
-          Q.x -= other.x
-          Q.y -= other.y
-          Q.z -= other.z
-          return Q
       else:
-          return self.copy()
+          return NotImplemented
 
     def __isub__(self, other):
       """In-place subtraction"""
       if isinstance(other, Quaternion):
-          self.w -= other.w
+          self.q -= other.q
-          self.x -= other.x
+          self.p -= other.p
-          self.y -= other.y
-          self.z -= other.z
       return self
 
     def __neg__(self):
       """Additive inverse"""
-      self.w = -self.w
+      self.q = -self.q
-      self.x = -self.x
+      self.p = -self.p
-      self.y = -self.y
-      self.z = -self.z
       return self
 
     def __abs__(self):
       """Norm"""
-      return math.sqrt(self.w ** 2 + \
+      return math.sqrt(self.q ** 2 + np.dot(self.p,self.p))
-                       self.x ** 2 + \
-                       self.y ** 2 + \
-                       self.z ** 2)
 
     magnitude = __abs__
 
     def __eq__(self,other):
       """Equal at e-8 precision"""
-      return (abs(self.w-other.w) < 1e-8 and \
+      return (self-other).magnitude() < 1e-8 or (-self-other).magnitude() < 1e-8
-              abs(self.x-other.x) < 1e-8 and \
-              abs(self.y-other.y) < 1e-8 and \
-              abs(self.z-other.z) < 1e-8) \
-              or \
-             (abs(-self.w-other.w) < 1e-8 and \
-              abs(-self.x-other.x) < 1e-8 and \
-              abs(-self.y-other.y) < 1e-8 and \
-              abs(-self.z-other.z) < 1e-8)
 
     def __ne__(self,other):
       """Not equal at e-8 precision"""
@@ -255,31 +187,26 @@ class Quaternion:
 
     def __cmp__(self,other):
       """Linear ordering"""
-      return (self.Rodrigues()>other.Rodrigues()) - (self.Rodrigues()<other.Rodrigues())
+      return (1 if np.linalg.norm(self.asRodrigues()) > np.linalg.norm(other.asRodrigues()) else 0) \
+           - (1 if np.linalg.norm(self.asRodrigues()) < np.linalg.norm(other.asRodrigues()) else 0)
 
     def magnitude_squared(self):
-      return self.w ** 2 + \
+      return self.q ** 2 + np.dot(self.p,self.p)
-             self.x ** 2 + \
-             self.y ** 2 + \
-             self.z ** 2
 
     def identity(self):
-      self.w = 1.
+      self.q = 1.
-      self.x = 0.
+      self.p = np.zeros(3,dtype=float)
-      self.y = 0.
-      self.z = 0.
       return self
 
     def normalize(self):
       d = self.magnitude()
       if d > 0.0:
-          self /= d
+          self.q /= d
+          self.p /= d
       return self
 
     def conjugate(self):
-      self.x = -self.x
+      self.p = -self.p
-      self.y = -self.y
-      self.z = -self.z
       return self
 
     def inverse(self):
@@ -290,11 +217,9 @@ class Quaternion:
       return self
 
     def homomorph(self):
-      if self.w < 0.0:
+      if self.q < 0.0:
-        self.w = -self.w
+        self.q = -self.q
-        self.x = -self.x
+        self.p = -self.p
-        self.y = -self.y
-        self.z = -self.z
       return self
 
     def normalized(self):
@@ -310,27 +235,35 @@ class Quaternion:
       return self.copy().homomorph()
 
     def asList(self):
-      return [i for i in self]
+      return [self.q]+list(self.p)
 
     def asM(self):                                                # to find Averaging Quaternions (see F. Landis Markley et al.)
-      return np.outer([i for i in self],[i for i in self])
+      return np.outer(self.asList(),self.asList())
       
     def asMatrix(self):
-      qbarhalf = 0.5*(self.w**2 - self.x**2 - self.y**2 - self.z**2)
+      # Rowenhorst_etal2015 MSMSE: value of P is selected as -1
+      P = -1.0
+      qbarhalf = 0.5*(self.q**2 - np.dot(self.p,self.p))
       return 2.0*np.array(
-        [[      qbarhalf + self.x**2    , self.x*self.y - self.w*self.z, self.x*self.z + self.w*self.y],
+        [[        qbarhalf + self.p[0]**2          ,
-         [ self.x*self.y + self.w*self.z,      qbarhalf + self.y**2    , self.y*self.z - self.w*self.x],
+           self.p[0]*self.p[1] -P* self.q*self.p[2],
-         [ self.x*self.z - self.w*self.y, self.y*self.z + self.w*self.x,      qbarhalf + self.z**2    ],
+           self.p[0]*self.p[2] +P* self.q*self.p[1] ],
+         [ self.p[0]*self.p[1] +P* self.q*self.p[2],
+                  qbarhalf + self.p[1]**2          ,
+           self.p[1]*self.p[2] -P* self.q*self.p[0] ],
+         [ self.p[0]*self.p[2] -P* self.q*self.p[1],
+           self.p[1]*self.p[2] +P* self.q*self.p[0],
+                  qbarhalf + self.p[2]**2           ],
         ])
 
     def asAngleAxis(self,
                     degrees = False):
-      if self.w > 1:
+      if self.q > 1.:
           self.normalize()
 
-      s = math.sqrt(1. - self.w**2)
+      s = math.sqrt(1. - self.q**2)
-      x = 2*self.w**2 - 1.
+      x = 2*self.q**2 - 1.
-      y = 2*self.w * s
+      y = 2*self.q * s
 
       angle = math.atan2(y,x)
       if angle < 0.0:
@@ -338,26 +271,28 @@ class Quaternion:
         s     *= -1.
 
       return (np.degrees(angle) if degrees else angle,
-              np.array([1.0, 0.0, 0.0] if np.abs(angle) < 1e-6 else [self.x / s, self.y / s, self.z / s]))
+              np.array([1.0, 0.0, 0.0] if np.abs(angle) < 1e-6 else self.p / s))
 
     def asRodrigues(self):
-      return np.inf*np.ones(3) if self.w == 0.0 else np.array([self.x, self.y, self.z])/self.w
+      return np.inf*np.ones(3) if self.q == 0.0 else self.p/self.q
 
     def asEulers(self,
                  degrees = False):
       """Orientation as Bunge-Euler angles."""
-      q03 = self.w**2+self.z**2
+      # Rowenhorst_etal2015 MSMSE: value of P is selected as -1
-      q12 = self.x**2+self.y**2
+      P   = -1.0
+      q03 = self.q**2 + self.p[2]**2
+      q12 = self.p[0]**2 + self.p[1]**2
       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])
+        eulers = np.array([math.atan2(-2*P*self.q*self.p[2],self.q**2-self.p[2]**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])
+        eulers = np.array([math.atan2( 2  *self.p[0]*self.p[1],self.p[0]**2-self.p[1]**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),
+        eulers = np.array([math.atan2((self.p[0]*self.p[2]-P*self.q*self.p[1])/chi,(-P*self.q*self.p[0]-self.p[1]*self.p[2])/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),
+                           math.atan2((P*self.q*self.p[1]+self.p[0]*self.p[2])/chi,( self.p[1]*self.p[2]-P*self.q*self.p[0])/chi),
                           ])
 
       return np.degrees(eulers) if degrees else eulers
@@ -371,25 +306,26 @@ class Quaternion:
 
     @classmethod
     def fromRandom(cls,randomSeed = None):
+      import binascii
       if randomSeed is None:
-        randomSeed = int(os.urandom(4).encode('hex'), 16)
+        randomSeed = int(binascii.hexlify(os.urandom(4)),16)
       np.random.seed(randomSeed)
       r = np.random.random(3)
       w = math.cos(2.0*math.pi*r[0])*math.sqrt(r[2])
       x = math.sin(2.0*math.pi*r[1])*math.sqrt(1.0-r[2])
       y = math.cos(2.0*math.pi*r[1])*math.sqrt(1.0-r[2])
       z = math.sin(2.0*math.pi*r[0])*math.sqrt(r[2])
-      return cls([w,x,y,z])
+      return cls(quat=[w,x,y,z])
 
 
     @classmethod
     def fromRodrigues(cls, rodrigues):
       if not isinstance(rodrigues, np.ndarray): rodrigues = np.array(rodrigues)
-      halfangle = math.atan(np.linalg.norm(rodrigues))
+      norm = np.linalg.norm(rodrigues)
+      halfangle = math.atan(norm)
+      s = math.sin(halfangle)
       c = math.cos(halfangle)
-      w = c
+      return cls(q=c,p=s*rodrigues/norm)
-      x,y,z = rodrigues/c
-      return cls([w,x,y,z])
 
 
     @classmethod
@@ -397,22 +333,19 @@ class Quaternion:
                       angle,
                       axis,
                       degrees = False):
-      if not isinstance(axis, np.ndarray): axis = np.array(axis,dtype='d')
+      if not isinstance(axis, np.ndarray): axis = np.array(axis,dtype=float)
       axis = axis.astype(float)/np.linalg.norm(axis)
       angle = np.radians(angle) if degrees else angle
       s = math.sin(0.5 * angle)
-      w = math.cos(0.5 * angle)
+      c = math.cos(0.5 * angle)
-      x = axis[0] * s
+      return cls(q=c,p=axis*s)
-      y = axis[1] * s
-      z = axis[2] * s
-      return cls([w,x,y,z])
 
 
     @classmethod
     def fromEulers(cls,
                    eulers,
                    degrees = False):
-      if not isinstance(eulers, np.ndarray): eulers = np.array(eulers,dtype='d')
+      if not isinstance(eulers, np.ndarray): eulers = np.array(eulers,dtype=float)
       eulers = np.radians(eulers) if degrees else eulers
 
       sigma = 0.5*(eulers[0]+eulers[2])
@@ -420,11 +353,13 @@ class Quaternion:
       c = np.cos(0.5*eulers[1])
       s = np.sin(0.5*eulers[1])
 
+      # Rowenhorst_etal2015 MSMSE: value of P is selected as -1
+      P = -1.0
       w =       c * np.cos(sigma) 
-      x = -s * np.cos(delta) 
+      x =  -P * s * np.cos(delta) 
-      y = -s * np.sin(delta)
+      y =  -P * s * np.sin(delta)
-      z = -c * np.sin(sigma) 
+      z =  -P * c * np.sin(sigma) 
-      return cls([w,x,y,z])
+      return cls(quat=[w,x,y,z])
 
 
 # Modified Method to calculate Quaternion from Orientation Matrix,
@@ -435,16 +370,18 @@ class Quaternion:
       if m.shape != (3,3) and np.prod(m.shape) == 9:
         m = m.reshape(3,3)
 
+      # Rowenhorst_etal2015 MSMSE: value of P is selected as -1
+      P = -1.0
       w =  0.5*math.sqrt(1.+m[0,0]+m[1,1]+m[2,2])
-      x = 0.5*math.sqrt(1.+m[0,0]-m[1,1]-m[2,2])
+      x = P*0.5*math.sqrt(1.+m[0,0]-m[1,1]-m[2,2])
-      y = 0.5*math.sqrt(1.-m[0,0]+m[1,1]-m[2,2])
+      y = P*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])
+      z = P*0.5*math.sqrt(1.-m[0,0]-m[1,1]+m[2,2])
 
       x *= -1 if m[2,1] < m[1,2] else 1
       y *= -1 if m[0,2] < m[2,0] else 1
       z *= -1 if m[1,0] < m[0,1] else 1
 
-      return cls( np.array([w,x,y,z])/math.sqrt(w**2 + x**2 + y**2 + z**2))
+      return cls(quat=np.array([w,x,y,z])/math.sqrt(w**2 + x**2 + y**2 + z**2))
 
 
     @classmethod
@@ -458,36 +395,30 @@ class Quaternion:
         assert isinstance(q1, Quaternion) and isinstance(q2, Quaternion)
         Q = cls()
 
-        costheta = q1.w * q2.w + q1.x * q2.x + q1.y * q2.y + q1.z * q2.z
+        costheta = q1.q*q2.q + np.dot(q1.p,q2.p)
         if costheta < 0.:
             costheta = -costheta
             q1 = q1.conjugated()
-        elif costheta > 1:
+        elif costheta > 1.:
-            costheta = 1
+            costheta = 1.
 
         theta = math.acos(costheta)
         if abs(theta) < 0.01:
-            Q.w = q2.w
+            Q.q = q2.q
-            Q.x = q2.x
+            Q.p = q2.p
-            Q.y = q2.y
-            Q.z = q2.z
             return Q
 
         sintheta = math.sqrt(1.0 - costheta * costheta)
         if abs(sintheta) < 0.01:
-            Q.w = (q1.w + q2.w) * 0.5
+            Q.q = (q1.q + q2.q) * 0.5
-            Q.x = (q1.x + q2.x) * 0.5
+            Q.p = (q1.p + q2.p) * 0.5
-            Q.y = (q1.y + q2.y) * 0.5
-            Q.z = (q1.z + q2.z) * 0.5
             return Q
 
-        ratio1 = math.sin((1 - t) * theta) / sintheta
+        ratio1 = math.sin((1.0 - t) * theta) / sintheta
-        ratio2 = math.sin(t * theta) / sintheta
+        ratio2 = math.sin(       t  * theta) / sintheta
 
-        Q.w = q1.w * ratio1 + q2.w * ratio2
+        Q.q = q1.q * ratio1 + q2.q * ratio2
-        Q.x = q1.x * ratio1 + q2.x * ratio2
+        Q.p = q1.p * ratio1 + q2.p * ratio2
-        Q.y = q1.y * ratio1 + q2.y * ratio2
-        Q.z = q1.z * ratio1 + q2.z * ratio2
         return Q
 
 
@@ -513,7 +444,7 @@ class Symmetry:
 
   def __repr__(self):
     """Readbable string"""
-    return '%s' % (self.lattice)
+    return '{}'.format(self.lattice)
 
 
   def __eq__(self, other):
@@ -722,7 +653,7 @@ class Symmetry:
       else:
         return True
 
-    v = np.array(vector,dtype = float)
+    v = np.array(vector,dtype=float)
     if proper:                                                                                      # check both improper ...
       theComponents = np.dot(basis['improper'],v)
       inSST = np.all(theComponents >= 0.0)
@@ -737,10 +668,10 @@ class Symmetry:
     if color:                                                                                       # have to return color array
       if inSST:
         rgb = np.power(theComponents/np.linalg.norm(theComponents),0.5)                             # smoothen color ramps
-        rgb = np.minimum(np.ones(3,'d'),rgb)                                                        # limit to maximum intensity
+        rgb = np.minimum(np.ones(3,dtype=float),rgb)                                                        # limit to maximum intensity
         rgb /= max(rgb)                                                                             # normalize to (HS)V = 1
       else:
-        rgb = np.zeros(3,'d')
+        rgb = np.zeros(3,dtype=float)
       return (inSST,rgb)
     else:
       return inSST
@@ -780,8 +711,9 @@ class Orientation:
       self.quaternion = Quaternion.fromRodrigues(Rodrigues)
     elif isinstance(quaternion, Quaternion):                                                        # based on given quaternion
       self.quaternion = quaternion.homomorphed()
-    elif isinstance(quaternion, np.ndarray) and quaternion.shape == (4,):                           # based on given quaternion-like array
+    elif (isinstance(quaternion, np.ndarray) and quaternion.shape == (4,)) or \
-      self.quaternion = Quaternion(quaternion).homomorphed()
+         (isinstance(quaternion, list)       and len(quaternion)  ==  4  ):                         # based on given quaternion-like array
+      self.quaternion = Quaternion(quat=quaternion).homomorphed()
 
     self.symmetry = Symmetry(symmetry)
 
@@ -794,10 +726,12 @@ class Orientation:
 
   def __repr__(self):
     """Value as all implemented representations"""
-    return 'Symmetry: %s\n' % (self.symmetry) + \
+    return '\n'.join([
-           'Quaternion: %s\n' % (self.quaternion) + \
+              'Symmetry: {}'.format(self.symmetry),
-           'Matrix:\n%s\n' % ( '\n'.join(['\t'.join(map(str,self.asMatrix()[i,:])) for i in range(3)]) ) + \
+              'Quaternion: {}'.format(self.quaternion),
-           'Bunge Eulers / deg: %s' % ('\t'.join(map(str,self.asEulers(degrees=True))) )
+              'Matrix:\n{}'.format( '\n'.join(['\t'.join(list(map(str,self.asMatrix()[i,:]))) for i in range(3)]) ),
+              'Bunge Eulers / deg: {}'.format('\t'.join(list(map(str,self.asEulers(degrees=True)))) ),
+              ])
 
   def asQuaternion(self):
     return self.quaternion.asList()
@@ -927,7 +861,7 @@ class Orientation:
       M = closest.quaternion.asM() * n if i == 0 else M + closest.quaternion.asM() * n              # noqa add (multiples) of this orientation to average noqa
     eig, vec = np.linalg.eig(M/N)
 
-    return Orientation(quaternion = Quaternion(quatArray = np.real(vec.T[eig.argmax()])),
+    return Orientation(quaternion = Quaternion(quat = np.real(vec.T[eig.argmax()])),
                        symmetry = reference.symmetry.lattice)
 
 

processing/post/addGrainID.py

@@ -200,9 +200,9 @@ for name in filenames:
       if gID != -1 and gID not in alreadyChecked:                                                   # indexed point belonging to a grain not yet tested?
         alreadyChecked[gID] = True                                                                  # remember not to check again
         disorientation = o.disorientation(orientations[gID],SST = False)[0]                         # compare against other orientation
-        if disorientation.quaternion.w >  cos_disorientation:                                       # within threshold ...
+        if disorientation.quaternion.q >  cos_disorientation:                                       # within threshold ...
           candidates.append(gID)                                                                    # remember as potential candidate
-          if disorientation.quaternion.w >= bestDisorientation.w:                                   # ... and better than current best? 
+          if disorientation.quaternion.q >= bestDisorientation.q:                                   # ... and better than current best? 
             matched = True
             matchedID = gID                                                                         # remember that grain
             bestDisorientation = disorientation.quaternion

processing/pre/geom_addPrimitive.py

@@ -64,11 +64,11 @@ if options.dimension is None:
   parser.error('no dimension specified.')   
 if options.angleaxis is not None:
   options.angleaxis = list(map(float,options.angleaxis))
-  rotation = damask.Quaternion().fromAngleAxis(np.radians(options.angleaxis[0]) if options.degrees else options.angleaxis[0],
+  rotation = damask.Quaternion.fromAngleAxis(np.radians(options.angleaxis[0]) if options.degrees else options.angleaxis[0],
                                              options.angleaxis[1:4])
 elif options.quaternion is not None:
-  options.quaternion = map(float,options.quaternion)
+  options.quaternion = list(map(float,options.quaternion))
-  rotation = damask.Quaternion(options.quaternion)
+  rotation = damask.Quaternion(quat=options.quaternion)
 else:
   rotation = damask.Quaternion()
 

processing/pre/geom_fromTable.py

@@ -43,7 +43,7 @@ parser.add_option('-e', '--eulers',
 parser.add_option('-d', '--degrees',
                   dest = 'degrees',
                   action = 'store_true',
-                  help = 'angles are given in degrees [%default]')
+                  help = 'all angles are in degrees')
 parser.add_option('-m', '--matrix',
                   dest = 'matrix',
                   type = 'string', metavar = 'string',
@@ -71,7 +71,7 @@ parser.add_option('--axes',
 parser.add_option('-s', '--symmetry',
                   dest = 'symmetry',
                   action = 'extend', metavar = '<string LIST>',
-                  help = 'crystal symmetry %default {{{}}} '.format(', '.join(damask.Symmetry.lattices[1:])))
+                  help = 'crystal symmetry of each phase %default {{{}}} '.format(', '.join(damask.Symmetry.lattices[1:])))
 parser.add_option('--homogenization',
                   dest = 'homogenization',
                   type = 'int', metavar = 'int',
@@ -234,7 +234,7 @@ for name in filenames:
             o = damask.Orientation(Eulers = myData[colOri:colOri+3]*toRadians,
                                    symmetry = mySym)
           elif inputtype == 'matrix':
-            o = damask.Orientation(matrix = myData[colOri:colOri+9].reshape(3,3).transpose(),
+            o = damask.Orientation(matrix = myData[colOri:colOri+9].reshape(3,3),
                                    symmetry = mySym)
           elif inputtype == 'frame':
             o = damask.Orientation(matrix = np.hstack((myData[colOri[0]:colOri[0]+3],
@@ -246,7 +246,7 @@ for name in filenames:
             o = damask.Orientation(quaternion = myData[colOri:colOri+4],
                                    symmetry = mySym)
           
-          cos_disorientations = -np.ones(1,dtype='f')                                               # largest possible disorientation
+          cos_disorientations = -np.ones(1,dtype=float)                                             # largest possible disorientation
           closest_grain = -1                                                                        # invalid neighbor
 
           if options.tolerance > 0.0:                                                               # only try to compress orientations if asked to
@@ -258,7 +258,7 @@ for name in filenames:
 
             if len(grains) > 0:                                                                     # check immediate neighborhood first
               cos_disorientations = np.array([o.disorientation(orientations[grainID],
-                                                               SST = False)[0].quaternion.w \
+                                                               SST = False)[0].quaternion.q \
                                               for grainID in grains])                               # store disorientation per grainID
               closest_grain = np.argmax(cos_disorientations)                                        # grain among grains with closest orientation to myself
               match = 'local'
@@ -269,7 +269,7 @@ for name in filenames:
 
               if len(grains) > 0:
                 cos_disorientations = np.array([o.disorientation(orientations[grainID],
-                                                                 SST = False)[0].quaternion.w \
+                                                                 SST = False)[0].quaternion.q \
                                                 for grainID in grains])                             # store disorientation per grainID
                 closest_grain = np.argmax(cos_disorientations)                                      # grain among grains with closest orientation to myself
                 match = 'global'

processing/pre/hybridIA_linODFsampling.py

@@ -244,7 +244,7 @@ for name in filenames:
     continue
   damask.util.report(scriptName,name)
 
-  randomSeed = int(os.urandom(4).encode('hex'), 16)  if options.randomSeed is None else options.randomSeed         # random seed per file for second phase
+  randomSeed = int(os.urandom(4).hex(), 16)  if options.randomSeed is None else options.randomSeed         # random seed per file for second phase
   random.seed(randomSeed)
 
 # ------------------------------------------ read header and data ---------------------------------

src/IO.f90

@@ -191,7 +191,9 @@ recursive function IO_recursiveRead(fileName,cnt) result(fileContent)
     l,i, &
     myStat
 
-  if (merge(cnt,0_pInt,present(cnt))>10_pInt) call IO_error(106_pInt,ext_msg=trim(fileName))
+  if (present(cnt)) then
+    if (cnt>10_pInt) call IO_error(106_pInt,ext_msg=trim(fileName))
+  endif
 
 !--------------------------------------------------------------------------------------------------
 ! read data as stream
@@ -684,7 +686,11 @@ function IO_stringValue(string,chunkPos,myChunk,silent)
 
  logical                                                  :: warn
 
- warn = merge(silent,.false.,present(silent))
+ if (present(silent)) then
+   warn = silent
+ else
+   warn = .false.
+ endif
 
  IO_stringValue = ''
  valuePresent: if (myChunk > chunkPos(1) .or. myChunk < 1_pInt) then

src/config.f90

@@ -513,8 +513,12 @@ character(len=65536) function getString(this,key,defaultVal,raw)
  type(tPartitionedStringList),  pointer                      :: item
  logical                                                     :: found, &
                                                                 whole
+ if (present(raw)) then
+   whole = raw
+ else
+   whole = .false.
+ endif
 
- whole = merge(raw,.false.,present(raw))                                                            ! whole string or white space splitting
  found = present(defaultVal)
  if (found) then
    getString = trim(defaultVal)
@@ -661,7 +665,11 @@ function getStrings(this,key,defaultVal,requiredShape,raw)
                                                              cumulative
 
  cumulative = (key(1:1) == '(' .and. key(len_trim(key):len_trim(key)) == ')')
- whole = merge(raw,.false.,present(raw))
+ if (present(raw)) then
+   whole = raw
+ else
+   whole = .false.
+ endif
  found = .false.
 
  item => this

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/kinematics_thermal_expansion.f90

@@ -225,7 +225,7 @@ subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar3333, ipc,
              + lattice_thermalExpansion33(1:3,1:3,2,phase)*(T - TRef)**1 &                           ! linear    coefficient
              + lattice_thermalExpansion33(1:3,1:3,3,phase)*(T - TRef)**2 &                           ! quadratic coefficient
              ) / &
-      (1.0_pReal \
+      (1.0_pReal &
             + lattice_thermalExpansion33(1:3,1:3,1,phase)*(T - TRef)**1 / 1. &
             + lattice_thermalExpansion33(1:3,1:3,2,phase)*(T - TRef)**2 / 2. &
             + lattice_thermalExpansion33(1:3,1:3,3,phase)*(T - TRef)**3 / 3. &

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/math.f90

@@ -2625,12 +2625,9 @@ real(pReal) pure function math_clip(a, left, right)
  real(pReal), intent(in) :: a
  real(pReal), intent(in), optional :: left, right
 
-   
+ math_clip = a
- math_clip = min ( &
+ if (present(left))   math_clip = max(left,math_clip)
-                   max (merge(left, -huge(a), present(left)), a), &
+ if (present(right))  math_clip = min(right,math_clip)
-                        merge(right, huge(a), present(right)) &
-                  )
-
  if (present(left) .and. present(right)) &
    math_clip = merge (IEEE_value(1.0_pReal,IEEE_quiet_NaN),math_clip, left>right)
 

src/spectral_damage.f90

@@ -120,8 +120,8 @@ subroutine spectral_damage_init()
      trim(snes_type) == 'vinewtonssls') then
    call DMGetGlobalVector(damage_grid,lBound,ierr); CHKERRQ(ierr)
    call DMGetGlobalVector(damage_grid,uBound,ierr); CHKERRQ(ierr)
-   call VecSet(lBound,0.0,ierr); CHKERRQ(ierr)
+   call VecSet(lBound,0.0_pReal,ierr); CHKERRQ(ierr)
-   call VecSet(uBound,1.0,ierr); CHKERRQ(ierr)
+   call VecSet(uBound,1.0_pReal,ierr); CHKERRQ(ierr)
    call SNESVISetVariableBounds(damage_snes,lBound,uBound,ierr)                                     !< variable bounds for variational inequalities like contact mechanics, damage etc.
    call DMRestoreGlobalVector(damage_grid,lBound,ierr); CHKERRQ(ierr)
    call DMRestoreGlobalVector(damage_grid,uBound,ierr); CHKERRQ(ierr)
@@ -134,7 +134,7 @@ subroutine spectral_damage_init()
  xend = xstart + xend - 1
  yend = ystart + yend - 1
  zend = zstart + zend - 1 
- call VecSet(solution,1.0,ierr); CHKERRQ(ierr)
+ call VecSet(solution,1.0_pReal,ierr); CHKERRQ(ierr)
  allocate(damage_current(grid(1),grid(2),grid3), source=1.0_pReal)
  allocate(damage_lastInc(grid(1),grid(2),grid3), source=1.0_pReal)
  allocate(damage_stagInc(grid(1),grid(2),grid3), source=1.0_pReal)