Merge branch 'development' into Dislotwin-climb2

PRIVATE

@@ -1 +1 @@
-Subproject commit 18a976753be06aca6e15f580998e713daa08bb41
+Subproject commit 5c5adbd8ccc0210fd6507431db8ec82ecec75352

VERSION

@@ -1 +1 @@
-v2.0.3-614-g4d6a047b
+v2.0.3-625-gbace77db

installation/symlink_Processing.py

@@ -1,8 +1,9 @@
 #!/usr/bin/env python3
-# -*- coding: UTF-8 no BOM -*-
 
-# Makes postprocessing routines acessible from everywhere.
+# Makes postprocessing routines accessible from everywhere.
-import os,sys
+import os
+import sys
+
 import damask
 
 damaskEnv = damask.Environment()

processing/post/addOrientations.py

@@ -27,6 +27,7 @@ Additional (globally fixed) rotations of the lab frame and/or crystal frame can
 representations = {
                   'quaternion': ['qu',4],
                   'rodrigues':  ['ro',4],
+                  'Rodrigues':  ['Ro',3],
                   'eulers':     ['eu',3],
                   'matrix':     ['om',9],
                   'angleaxis':  ['ax',4],
@@ -80,15 +81,20 @@ parser.add_option('-z',
                   dest = 'z',
                   metavar = 'string',
                   help = 'label of lab z vector (expressed in crystal coords)')
+parser.add_option('--lattice',
+                  dest = 'lattice',
+                  metavar = 'string',
+                  help = 'lattice structure to reduce rotation into fundamental zone')
 
 parser.set_defaults(output = [],
                     labrotation     = (1.,1.,1.,0.),                                                # no rotation about (1,1,1)
                     crystalrotation = (1.,1.,1.,0.),                                                # no rotation about (1,1,1)
+                    lattice = None,
                    )
 
 (options, filenames) = parser.parse_args()
 
-options.output = list(map(lambda x: x.lower(), options.output))
+#options.output = list(map(lambda x: x.lower(), options.output))
 if options.output == [] or (not set(options.output).issubset(set(representations))):
   parser.error('output must be chosen from {}.'.format(', '.join(representations)))
 
@@ -121,7 +127,7 @@ if filenames == []: filenames = [None]
 for name in filenames:
   try:    table = damask.ASCIItable(name = name,
                                     buffered = False)
-  except: continue
+  except Exception: continue
   damask.util.report(scriptName,name)
 
 # ------------------------------------------ read header ------------------------------------------
@@ -156,16 +162,16 @@ for name in filenames:
   outputAlive = True
   while outputAlive and table.data_read():                                                          # read next data line of ASCII table
     if   inputtype == 'eulers':
-      l = representations['eulers'][1]
+      d = representations['eulers'][1]
-      o = damask.Rotation.fromEulers(list(map(float,table.data[column:column+l])),options.degrees)
+      o = damask.Rotation.fromEulers(list(map(float,table.data[column:column+d])),options.degrees)
       
     elif inputtype == 'rodrigues':
-      l = representations['rodrigues'][1]
+      d = representations['rodrigues'][1]
-      o = damask.Rotation.fromRodrigues(list(map(float,table.data[column:column+l])))
+      o = damask.Rotation.fromRodrigues(list(map(float,table.data[column:column+d])))
       
     elif inputtype == 'matrix':
-      l = representations['matrix'][1]
+      d = representations['matrix'][1]
-      o = damask.Rotation.fromMatrix(list(map(float,table.data[column:column+l])))
+      o = damask.Rotation.fromMatrix(list(map(float,table.data[column:column+d])))
 
     elif inputtype == 'frame':
       M = np.array(list(map(float,table.data[column[0]:column[0]+3] + \
@@ -174,14 +180,18 @@ for name in filenames:
       o = damask.Rotation.fromMatrix(M/np.linalg.norm(M,axis=0))
       
     elif inputtype == 'quaternion':
-      l = representations['quaternion'][1]
+      d = representations['quaternion'][1]
-      o = damask.Rotation.fromQuaternion(list(map(float,table.data[column:column+l])))
+      o = damask.Rotation.fromQuaternion(list(map(float,table.data[column:column+d])))
 
-    o= r*o*R                                                                                        # apply additional lab and crystal frame rotations
+    o = r*o*R                                                                                       # apply additional lab and crystal frame rotations
+
+    if options.lattice is not None:
+      o = damask.Orientation(rotation = o,lattice = options.lattice).reduced().rotation
 
     for output in options.output:
       if   output == 'quaternion': table.data_append(o.asQuaternion())
       elif output == 'rodrigues':  table.data_append(o.asRodrigues())
+      elif output == 'Rodrigues':  table.data_append(o.asRodrigues(vector=True))
       elif output == 'eulers':     table.data_append(o.asEulers(degrees=options.degrees))
       elif output == 'matrix':     table.data_append(o.asMatrix())
       elif output == 'angleaxis':  table.data_append(o.asAxisAngle(degrees=options.degrees))

python/damask/Lambert.py

@@ -1,5 +1,6 @@
 ####################################################################################################
-# Code below available according to the followin conditions on https://github.com/MarDiehl/3Drotations
+# Code below available according to the following conditions on
+# https://github.com/MarDiehl/3Drotations
 ####################################################################################################
 # Copyright (c) 2017-2019, Martin Diehl/Max-Planck-Institut für Eisenforschung GmbH
 # Copyright (c) 2013-2014, Marc De Graef/Carnegie Mellon University
@@ -36,7 +37,20 @@ beta = np.pi**(5./6.)/6.**(1./6.)/2.
 R1   = (3.*np.pi/4.)**(1./3.)
 
 def CubeToBall(cube):
-  
+    """
+    Map a point in a uniform refinable cubical grid to a point on a uniform refinable grid on a ball.
+
+    Parameters
+    ----------
+    cube : numpy.ndarray
+        coordinates of a point in a uniform refinable cubical grid.
+        
+    References
+    ----------
+    D. Roşca et al., Modelling and Simulation in Materials Science and Engineering 22:075013, 2014
+    https://doi.org/10.1088/0965-0393/22/7/075013
+
+    """
     if np.abs(np.max(cube))>np.pi**(2./3.) * 0.5:
       raise ValueError
   
@@ -45,7 +59,7 @@ def CubeToBall(cube):
       ball = np.zeros(3)
     else:
       # get pyramide and scale by grid parameter ratio
-      p = GetPyramidOrder(cube)
+      p = get_order(cube)
       XYZ = cube[p] * sc
 
       # intercept all the points along the z-axis
@@ -74,7 +88,20 @@ def CubeToBall(cube):
 
 
 def BallToCube(ball):
-  
+    """
+    Map a point on a uniform refinable grid on a ball to a point in a uniform refinable cubical grid.
+
+    Parameters
+    ----------
+    ball : numpy.ndarray
+       coordinates of a point on a uniform refinable grid on a ball.
+
+    References
+    ----------
+    D. Roşca et al., Modelling and Simulation in Materials Science and Engineering 22:075013, 2014
+    https://doi.org/10.1088/0965-0393/22/7/075013
+
+    """
     rs = np.linalg.norm(ball)
     if rs > R1: 
       raise ValueError
@@ -82,7 +109,7 @@ def BallToCube(ball):
     if np.allclose(ball,0.0,rtol=0.0,atol=1.0e-300):
       cube = np.zeros(3)
     else:
-      p = GetPyramidOrder(ball)
+      p = get_order(ball)
       xyz3 = ball[p] 
 
       # inverse M_3
@@ -110,8 +137,24 @@ def BallToCube(ball):
       
     return cube
 
-def GetPyramidOrder(xyz):
+def get_order(xyz):
- 
+    """
+    Get order of the coordinates.
+
+    Depending on the pyramid in which the point is located, the order need to be adjusted.
+    
+    Parameters
+    ----------
+    xyz : numpy.ndarray
+       coordinates of a point on a uniform refinable grid on a ball or
+       in a uniform refinable cubical grid.
+
+    References
+    ----------
+    D. Roşca et al., Modelling and Simulation in Materials Science and Engineering 22:075013, 2014
+    https://doi.org/10.1088/0965-0393/22/7/075013
+
+    """
     if   (abs(xyz[0])<= xyz[2]) and (abs(xyz[1])<= xyz[2]) or \
          (abs(xyz[0])<=-xyz[2]) and (abs(xyz[1])<=-xyz[2]):
       return [0,1,2]

python/damask/colormaps.py

@@ -3,13 +3,7 @@ import math
 import numpy as np
 
 class Color():
-  """
+  """Color representation in and conversion between different color-spaces."""
-  Conversion of colors between different color-spaces.
-
-  Colors should be given in the form Color('model',[vector]).
-  To convert or copy color from one space to other, use the methods
-  convertTo('model') or expressAs('model'), respectively.
-  """
 
   __slots__ = [
                'model',
@@ -22,7 +16,17 @@ class Color():
   def __init__(self,
                model = 'RGB',
                color = np.zeros(3,'d')):
+    """
+    Create a Color object.
+    
+    Parameters
+    ----------
+    model : string
+        color model
+    color : numpy.ndarray
+        vector representing the color according to the selected model
 
+    """
     self.__transforms__ = \
                    {'HSV':    {'index': 0, 'next': self._HSV2HSL},
                     'HSL':    {'index': 1, 'next': self._HSL2RGB,     'prev': self._HSL2HSV},
@@ -49,18 +53,27 @@ class Color():
 
 # ------------------------------------------------------------------
   def __repr__(self):
-    """Color model and values"""
+    """Color model and values."""
     return 'Model: %s Color: %s'%(self.model,str(self.color))
 
 
 # ------------------------------------------------------------------
   def __str__(self):
-    """Color model and values"""
+    """Color model and values."""
     return self.__repr__()
 
 
 # ------------------------------------------------------------------
   def convertTo(self,toModel = 'RGB'):
+    """
+    Change the color model permanently.
+    
+    Parameters
+    ----------
+    toModel : string
+        color model
+
+    """
     toModel = toModel.upper()
     if toModel not in list(self.__transforms__.keys()): return
 
@@ -79,16 +92,25 @@ class Color():
 
 # ------------------------------------------------------------------
   def expressAs(self,asModel = 'RGB'):
+    """
+    Return the color in a different model.
+    
+    Parameters
+    ----------
+    asModel : string
+        color model
+
+    """
     return self.__class__(self.model,self.color).convertTo(asModel)
 
 
 
   def _HSV2HSL(self):
     """
-    Convert H(ue) S(aturation) V(alue or brightness) to H(ue) S(aturation) L(uminance)
+    Convert H(ue) S(aturation) V(alue or brightness) to H(ue) S(aturation) L(uminance).
 
-    with all values in the range of 0 to 1
+    All values are in the range [0,1]
-    http://codeitdown.com/hsl-hsb-hsv-color/
+    http://codeitdown.com/hsl-hsb-hsv-color
     """
     if self.model != 'HSV': return
 
@@ -105,10 +127,10 @@ class Color():
 
   def _HSL2HSV(self):
     """
-    Convert H(ue) S(aturation) L(uminance) to H(ue) S(aturation) V(alue or brightness)
+    Convert H(ue) S(aturation) L(uminance) to H(ue) S(aturation) V(alue or brightness).
 
-    with all values in the range of 0 to 1
+    All values are in the range [0,1]
-    http://codeitdown.com/hsl-hsb-hsv-color/
+    http://codeitdown.com/hsl-hsb-hsv-color
     """
     if self.model != 'HSL': return
 
@@ -124,9 +146,9 @@ class Color():
 
   def _HSL2RGB(self):
     """
-    Convert H(ue) S(aturation) L(uminance) to R(red) G(reen) B(lue)
+    Convert H(ue) S(aturation) L(uminance) to R(red) G(reen) B(lue).
 
-    with all values in the range of 0 to 1
+    All values are in the range [0,1]
     from http://en.wikipedia.org/wiki/HSL_and_HSV
     """
     if self.model != 'HSL': return
@@ -150,9 +172,9 @@ class Color():
 
   def _RGB2HSL(self):
     """
-    Convert R(ed) G(reen) B(lue) to H(ue) S(aturation) L(uminance)
+    Convert R(ed) G(reen) B(lue) to H(ue) S(aturation) L(uminance).
 
-    with all values in the range of 0 to 1
+    All values are in the range [0,1]
     from http://130.113.54.154/~monger/hsl-rgb.html
     """
     if self.model != 'RGB': return
@@ -190,9 +212,9 @@ class Color():
 
   def _RGB2XYZ(self):
     """
-    Convert R(ed) G(reen) B(lue) to CIE XYZ
+    Convert R(ed) G(reen) B(lue) to CIE XYZ.
 
-    with all values in the range of 0 to 1
+    All values are in the range [0,1]
     from http://www.cs.rit.edu/~ncs/color/t_convert.html
     """
     if self.model != 'RGB': return
@@ -219,9 +241,9 @@ class Color():
 
   def _XYZ2RGB(self):
     """
-    Convert  CIE XYZ to R(ed) G(reen) B(lue)
+    Convert  CIE XYZ to R(ed) G(reen) B(lue).
 
-    with all values in the range of 0 to 1
+    All values are in the range [0,1]
     from http://www.cs.rit.edu/~ncs/color/t_convert.html
     """
     if self.model != 'XYZ':
@@ -251,9 +273,9 @@ class Color():
 
   def _CIELAB2XYZ(self):
     """
-    Convert  CIE Lab to CIE XYZ
+    Convert  CIE Lab to CIE XYZ.
 
-    with XYZ in the range of 0 to 1
+    All values are in the range [0,1]
     from http://www.easyrgb.com/index.php?X=MATH&H=07#text7
     """
     if self.model != 'CIELAB': return
@@ -275,11 +297,11 @@ class Color():
 
   def _XYZ2CIELAB(self):
     """
-    Convert CIE XYZ to CIE Lab
+    Convert CIE XYZ to CIE Lab.
 
-    with XYZ in the range of 0 to 1
+    All values are in the range [0,1]
     from http://en.wikipedia.org/wiki/Lab_color_space,
-    http://www.cs.rit.edu/~ncs/color/t_convert.html
+         http://www.cs.rit.edu/~ncs/color/t_convert.html
     """
     if self.model != 'XYZ': return
 
@@ -299,7 +321,7 @@ class Color():
 
   def _CIELAB2MSH(self):
     """
-    Convert CIE Lab to Msh colorspace
+    Convert CIE Lab to Msh colorspace.
 
     from http://www.cs.unm.edu/~kmorel/documents/ColorMaps/DivergingColorMapWorkshop.xls
     """
@@ -319,7 +341,7 @@ class Color():
 
   def _MSH2CIELAB(self):
     """
-    Convert Msh colorspace to CIE Lab
+    Convert Msh colorspace to CIE Lab.
 
     with s,h in radians
     from http://www.cs.unm.edu/~kmorel/documents/ColorMaps/DivergingColorMapWorkshop.xls
@@ -337,7 +359,7 @@ class Color():
 
 
 class Colormap():
-  """Perceptually uniform diverging or sequential colormaps."""
+  """Perceptually uniform diverging or sequential colormap."""
 
   __slots__ = [
                'left',
@@ -394,7 +416,21 @@ class Colormap():
                interpolate = 'perceptualuniform',
                predefined = None
                ):
+    """
+    Create a Colormap object.
+    
+    Parameters
+    ----------
+    left : Color
+        left color (minimum value)
+    right : Color
+        right color (maximum value)
+    interpolate : str
+        interpolation scheme (either 'perceptualuniform' or 'linear')
+    predefined : bool
+         ignore other arguments and use predefined definition
 
+    """
     if predefined is not None:
       left = self.__predefined__[predefined.lower()]['left']
       right= self.__predefined__[predefined.lower()]['right']
@@ -412,16 +448,22 @@ class Colormap():
 
 # ------------------------------------------------------------------
   def __repr__(self):
-    """Left and right value of colormap"""
+    """Left and right value of colormap."""
     return 'Left: %s Right: %s'%(self.left,self.right)
 
 
 # ------------------------------------------------------------------
   def invert(self):
+    """Switch left/minimum with right/maximum."""
     (self.left, self.right) = (self.right, self.left)
     return self
 
 
+# ------------------------------------------------------------------
+  def show_predefined(self):
+    """Show the labels of the predefined colormaps."""
+    print('\n'.join(self.__predefined__.keys()))
+
 # ------------------------------------------------------------------
   def color(self,fraction = 0.5):
 
@@ -490,21 +532,22 @@ class Colormap():
     frac = 0.5*(np.array(crop) + 1.0)                                                               # rescale crop range to fractions
     colors = [self.color(float(i)/(steps-1)*(frac[1]-frac[0])+frac[0]).expressAs(model).color for i in range(steps)]
     if   format == 'paraview':
-      colormap = ['[\n {{\n  "ColorSpace" : "RGB", "Name" : "{}",\n  "RGBPoints" : ['.format(name)] \
+      colormap = ['[\n {{\n  "ColorSpace": "RGB", "Name": "{}", "DefaultMap": true,\n  "RGBPoints" : ['.format(name)] \
-               + ['    {:4d},{:8.6f},{:8.6f},{:8.6f},'.format(i,color[0],color[1],color[2],)
+               + ['    {:4d},{:8.6f},{:8.6f},{:8.6f},'.format(i,color[0],color[1],color[2],) \
-                                                                        for i,color in enumerate(colors[:-1])]\
+                                                                        for i,color in enumerate(colors[:-1])] \
-               + ['    {:4d},{:8.6f},{:8.6f},{:8.6f} '.format(len(colors),colors[-1][0],colors[-1][1],colors[-1][2],)]\
+               + ['    {:4d},{:8.6f},{:8.6f},{:8.6f} '.format(len(colors),colors[-1][0],colors[-1][1],colors[-1][2],)] \
                + ['   ]\n }\n]']
+
     elif format == 'gmsh':
       colormap = ['View.ColorTable = {'] \
                + [',\n'.join(['{%s}'%(','.join([str(x*255.0) for x in color])) for color in colors])] \
                + ['}']
 
     elif format == 'gom':
-      colormap = ['1 1 ' + str(name) \
+      colormap = ['1 1 ' + str(name)
-                 + ' 9 ' + str(name) \
+                 + ' 9 ' + str(name)
-                 + ' 0 1 0 3 0 0 -1 9 \ 0 0 0 255 255 255 0 0 255 ' \
+                 + ' 0 1 0 3 0 0 -1 9 \\ 0 0 0 255 255 255 0 0 255 '
-                 + '30 NO_UNIT 1 1 64 64 64 255 1 0 0 0 0 0 0 3 0 ' + str(len(colors)) \
+                 + '30 NO_UNIT 1 1 64 64 64 255 1 0 0 0 0 0 0 3 0 ' + str(len(colors))
                  + ' '.join([' 0 %s 255 1'%(' '.join([str(int(x*255.0)) for x in color])) for color in reversed(colors)])]
 
     elif format == 'raw':

python/damask/geom.py

@@ -10,10 +10,27 @@ from . import version
 
 
 class Geom():
-  """Geometry definition for grid solvers"""
+  """Geometry definition for grid solvers."""
 
   def __init__(self,microstructure,size,origin=[0.0,0.0,0.0],homogenization=1,comments=[]):
-    """New geometry definition from array of microstructures and size"""
+    """
+    New geometry definition from array of microstructures and size.
+
+    Parameters
+    ----------
+    microstructure : numpy.ndarray
+        microstructure array (3D)
+    size : list or numpy.ndarray
+        physical size of the microstructure in meter.
+    origin : list or numpy.ndarray, optional
+        physical origin of the microstructure in meter.
+    homogenization : integer, optional
+        homogenization index.
+    comments : list of str, optional
+        comments lines.
+
+    """
+    self.__transforms__ = \
     self.set_microstructure(microstructure)
     self.set_size(size)
     self.set_origin(origin)
@@ -22,7 +39,7 @@ class Geom():
     
 
   def __repr__(self):
-    """Basic information on geometry definition"""
+    """Basic information on geometry definition."""
     return util.srepr([
            'grid     a b c:      {}'.format(' x '.join(map(str,self.get_grid  ()))),
            'size     x y z:      {}'.format(' x '.join(map(str,self.get_size  ()))),
@@ -33,7 +50,21 @@ class Geom():
           ])
 
   def update(self,microstructure=None,size=None,origin=None,rescale=False):
-    """Updates microstructure and size"""
+    """
+    Updates microstructure and size.
+
+    Parameters
+    ----------
+    microstructure : numpy.ndarray, optional
+        microstructure array (3D).
+    size : list or numpy.ndarray, optional
+        physical size of the microstructure in meter.
+    origin : list or numpy.ndarray, optional
+        physical origin of the microstructure in meter.
+    rescale : bool, optional
+        ignore size parameter and rescale according to change of grid points.
+
+    """
     grid_old    = self.get_grid()
     size_old    = self.get_size()
     origin_old  = self.get_origin()
@@ -44,46 +75,77 @@ class Geom():
       raise ValueError('Either set size explicitly or rescale automatically')
 
     self.set_microstructure(microstructure)
-    self.set_size(self.get_grid()/grid_old*self.size if rescale else size)
     self.set_origin(origin)
-   
+
+    if size is not None:
+      self.set_size(size)
+    elif rescale:
+      self.set_size(self.get_grid()/grid_old*self.size)
+
     message = ['grid     a b c:      {}'.format(' x '.join(map(str,grid_old)))]
     if np.any(grid_old != self.get_grid()):
       message[-1] = util.delete(message[-1])
-      message.append('grid     a b c:      {}'.format(' x '.join(map(str,self.get_grid()))))
+      message.append(util.emph('grid     a b c:      {}'.format(' x '.join(map(str,self.get_grid())))))
 
     message.append('size     x y z:      {}'.format(' x '.join(map(str,size_old))))
     if np.any(size_old != self.get_size()):
       message[-1] = util.delete(message[-1])
-      message.append('size     x y z:      {}'.format(' x '.join(map(str,self.get_size()))))
+      message.append(util.emph('size     x y z:      {}'.format(' x '.join(map(str,self.get_size())))))
 
     message.append('origin   x y z:      {}'.format('   '.join(map(str,origin_old))))
     if np.any(origin_old != self.get_origin()):
       message[-1] = util.delete(message[-1])
-      message.append('origin   x y z:      {}'.format('   '.join(map(str,self.get_origin()))))
+      message.append(util.emph('origin   x y z:      {}'.format('   '.join(map(str,self.get_origin())))))
 
     message.append('homogenization:      {}'.format(self.get_homogenization()))
 
     message.append('# microstructures:   {}'.format(unique_old))
     if unique_old != len(np.unique(self.microstructure)):
       message[-1] = util.delete(message[-1])
-      message.append('# microstructures:   {}'.format(len(np.unique(self.microstructure))))
+      message.append(util.emph('# microstructures:   {}'.format(len(np.unique(self.microstructure)))))
 
     message.append('max microstructure:  {}'.format(max_old))
     if max_old != np.nanmax(self.microstructure):
       message[-1] = util.delete(message[-1])
-      message.append('max microstructure:  {}'.format(np.nanmax(self.microstructure)))
+      message.append(util.emph('max microstructure:  {}'.format(np.nanmax(self.microstructure))))
     
     return util.return_message(message)
 
   def set_comments(self,comments):
+    """
+    Replaces all existing comments.
+
+    Parameters
+    ----------
+    comments : list of str
+        new comments.
+
+    """
     self.comments = []
     self.add_comments(comments)
     
   def add_comments(self,comments):
+    """
+    Appends comments to existing comments.
+
+    Parameters
+    ----------
+    comments : list of str
+        new comments.
+
+    """
     self.comments += [str(c) for c in comments] if isinstance(comments,list) else [str(comments)]
 
   def set_microstructure(self,microstructure):
+    """
+    Replaces the existing microstructure representation.
+
+    Parameters
+    ----------
+    microstructure : numpy.ndarray
+        microstructure array (3D).
+
+    """
     if microstructure is not None:
       if len(microstructure.shape) != 3:
         raise ValueError('Invalid microstructure shape {}'.format(*microstructure.shape))
@@ -93,6 +155,15 @@ class Geom():
         self.microstructure = np.copy(microstructure)
 
   def set_size(self,size):
+    """
+    Replaces the existing size information.
+
+    Parameters
+    ----------
+    size : list or numpy.ndarray
+        physical size of the microstructure in meter.
+
+    """
     if size is None:
       grid = np.asarray(self.microstructure.shape)
       self.size = grid/np.max(grid)
@@ -103,6 +174,15 @@ class Geom():
         self.size = np.array(size)
 
   def set_origin(self,origin):
+    """
+    Replaces the existing origin information.
+
+    Parameters
+    ----------
+    origin : list or numpy.ndarray
+        physical origin of the microstructure in meter
+
+    """
     if origin is not None:
       if len(origin) != 3:
         raise ValueError('Invalid origin {}'.format(*origin))
@@ -110,6 +190,15 @@ class Geom():
         self.origin = np.array(origin)
 
   def set_homogenization(self,homogenization):
+    """
+    Replaces the existing homogenization index.
+
+    Parameters
+    ----------
+    homogenization : integer
+        homogenization index
+
+    """
     if homogenization is not None:
       if not isinstance(homogenization,int) or homogenization < 1:
         raise TypeError('Invalid homogenization {}'.format(homogenization))
@@ -118,24 +207,31 @@ class Geom():
 
 
   def get_microstructure(self):
+    """Return the microstructure representation."""
     return np.copy(self.microstructure)
 
   def get_size(self):
+    """Return the physical size in meter."""
     return np.copy(self.size)
 
   def get_origin(self):
+    """Return the origin in meter."""
     return np.copy(self.origin)
 
   def get_grid(self):
+    """Return the grid discretization."""
     return np.array(self.microstructure.shape)
 
   def get_homogenization(self):
+    """Return the homogenization index."""
     return self.homogenization
 
   def get_comments(self):
+    """Return the comments."""
     return self.comments[:]
 
   def get_header(self):
+    """Return the full header (grid, size, origin, homogenization, comments)."""
     header =  ['{} header'.format(len(self.comments)+4)] + self.comments
     header.append('grid   a {} b {} c {}'.format(*self.get_grid()))
     header.append('size   x {} y {} z {}'.format(*self.get_size()))
@@ -145,7 +241,15 @@ class Geom():
       
   @classmethod
   def from_file(cls,fname):
-    """Reads a geom file"""
+    """
+    Reads a geom file.
+
+    Parameters
+    ----------
+    fname : str or file handle
+        geometry file to read.
+
+    """
     with (open(fname) if isinstance(fname,str) else fname) as f:
       f.seek(0)
       header_length,keyword = f.readline().split()[:2]
@@ -190,13 +294,21 @@ class Geom():
       raise TypeError('Invalid file: expected {} entries,found {}'.format(grid.prod(),i))
     
     microstructure = microstructure.reshape(grid,order='F')
-    if not np.any(np.mod(microstructure.flatten(),1) != 0.0):                                     # no float present
+    if not np.any(np.mod(microstructure.flatten(),1) != 0.0):                                       # no float present
       microstructure = microstructure.astype('int')
     
     return cls(microstructure.reshape(grid),size,origin,homogenization,comments)
 
   def to_file(self,fname):
-    """Writes to file"""
+    """
+    Writes a geom file.
+
+    Parameters
+    ----------
+    fname : str or file handle
+        geometry file to write.
+
+    """
     header = self.get_header()
     grid   = self.get_grid()
     format_string = '%{}i'.format(1+int(np.floor(np.log10(np.nanmax(self.microstructure))))) if self.microstructure.dtype == int \
@@ -208,7 +320,15 @@ class Geom():
 
                
   def to_vtk(self,fname=None):
-    """Generates vtk file. If file name is given, stored in file otherwise returned as string"""
+    """
+    Generates vtk file.
+
+    Parameters
+    ----------
+    fname : str, optional
+        vtk file to write. If no file is given, a string is returned.
+
+    """
     grid = self.get_grid() + np.ones(3,dtype=int)
     size = self.get_size()
     origin = self.get_origin()
@@ -262,7 +382,7 @@ class Geom():
 
                
   def show(self):
-    """Show raw content (as in file)"""
+    """Show raw content (as in file)."""
     f=StringIO()
     self.to_file(f)
     f.seek(0)

python/damask/orientation.py

@@ -10,35 +10,44 @@ from .quaternion import P
 ####################################################################################################
 class Rotation:
     u"""
-    Orientation stored as unit quaternion.
+    Orientation stored with functionality for conversion to different representations.
-    
-    Following: D Rowenhorst et al. Consistent representations of and conversions between 3D rotations
-               10.1088/0965-0393/23/8/083501
-    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 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, π]
-    Convention 6: P = -1 (as default)
-    
-    q is the real part, p = (x, y, z) are the imaginary parts.
 
+    References
+    ----------
+    D. Rowenhorst et al., Modelling and Simulation in Materials Science and Engineering 23:083501, 2015
+    https://doi.org/10.1088/0965-0393/23/8/083501
+
+    Conventions
+    -----------
+    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 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, π].
+    Convention 6: P = -1 (as default).
+    
+    Usage
+    -----
     Vector "a" (defined in coordinate system "A") is passively rotated
                resulting in new coordinates "b" when expressed in system "B".
     b = Q * a
     b = np.dot(Q.asMatrix(),a)
+
     """
 
     __slots__ = ['quaternion']
     
     def __init__(self,quaternion = np.array([1.0,0.0,0.0,0.0])):
       """
-      Initializes to identity unless specified
+      Initializes to identity unless specified.
       
-      If a quaternion is given, it needs to comply with the convection. Use .fromQuaternion
+      Parameters
-      to check the input.
+      ----------
+      quaternion : numpy.ndarray, optional
+          Unit quaternion that follows the conventions. Use .fromQuaternion to perform a sanity check.
+
       """
       if isinstance(quaternion,Quaternion):
         self.quaternion = quaternion.copy()
@@ -46,14 +55,14 @@ class Rotation:
         self.quaternion = Quaternion(q=quaternion[0],p=quaternion[1:4])
         
     def __copy__(self):
-      """Copy"""
+      """Copy."""
       return self.__class__(self.quaternion)
       
     copy = __copy__
 
 
     def __repr__(self):
-      """Value in selected representation"""
+      """Orientation displayed as unit quaternion, rotation matrix, and Bunge-Euler angles."""
       return '\n'.join([
             '{}'.format(self.quaternion),
             'Matrix:\n{}'.format( '\n'.join(['\t'.join(list(map(str,self.asMatrix()[i,:]))) for i in range(3)]) ),
@@ -62,9 +71,18 @@ class Rotation:
               
     def __mul__(self, other):
       """
-      Multiplication
+      Multiplication.
       
-      Rotation: Details needed (active/passive), rotation of (3,3,3,3)-matrix should be considered
+      Parameters
+      ----------
+      other : numpy.ndarray or Rotation
+          Vector, second or fourth order tensor, or rotation object that is rotated.
+
+      Todo
+      ----
+      Document details active/passive)
+      considere rotation of (3,3,3,3)-matrix
+
       """
       if isinstance(other, Rotation):                                                               # rotate a rotation
         return self.__class__(self.quaternion * other.quaternion).standardize()
@@ -104,32 +122,48 @@ class Rotation:
         
     
     def inverse(self):
-      """In-place inverse rotation/backward rotation"""
+      """In-place inverse rotation/backward rotation."""
       self.quaternion.conjugate()
       return self
       
     def inversed(self):
-      """Inverse rotation/backward rotation"""
+      """Inverse rotation/backward rotation."""
       return self.copy().inverse()
 
 
     def standardize(self):
-      """In-place quaternion representation with positive q"""
+      """In-place quaternion representation with positive q."""
       if self.quaternion.q < 0.0: self.quaternion.homomorph()
       return self
 
     def standardized(self):
-      """Quaternion representation with positive q"""
+      """Quaternion representation with positive q."""
       return self.copy().standardize()
 
 
     def misorientation(self,other):
-      """Misorientation"""
+      """
+      Get Misorientation.
+
+      Parameters
+      ----------
+      other : Rotation
+          Rotation to which the misorientation is computed.
+
+      """
       return other*self.inversed()
 
 
     def average(self,other):
-      """Calculate the average rotation"""
+      """
+      Calculate the average rotation.
+
+      Parameters
+      ----------
+      other : Rotation
+          Rotation from which the average is rotated.
+
+      """
       return Rotation.fromAverage([self,other])
       
 
@@ -137,41 +171,75 @@ class Rotation:
     # convert to different orientation representations (numpy arrays)
 
     def asQuaternion(self):
-      """Unit quaternion: (q, [p_1, p_2, p_3])"""
+      """Unit quaternion: (q, p_1, p_2, p_3)."""
       return self.quaternion.asArray()
       
     def asEulers(self,
                  degrees = False):
-      """Bunge-Euler angles: (φ_1, ϕ, φ_2)"""
+      """
+      Bunge-Euler angles: (φ_1, ϕ, φ_2).
+
+      Parameters
+      ----------
+      degrees : bool, optional
+          return angles in degrees.
+
+      """
       eu = qu2eu(self.quaternion.asArray())
       if degrees: eu = np.degrees(eu)
       return eu
     
     def asAxisAngle(self,
                     degrees = False):
-      """Axis-angle pair: ([n_1, n_2, n_3], ω)"""
+      """
+      Axis angle pair: ([n_1, n_2, n_3], ω).
+
+      Parameters
+      ----------
+      degrees : bool, optional
+          return rotation angle in degrees.
+      
+      """
       ax = qu2ax(self.quaternion.asArray())
       if degrees: ax[3] = np.degrees(ax[3])
       return ax
       
     def asMatrix(self):
-      """Rotation matrix"""
+      """Rotation matrix."""
       return qu2om(self.quaternion.asArray())
 
-    def asRodrigues(self,vector=False):
+    def asRodrigues(self,
-      """Rodrigues-Frank vector: ([n_1, n_2, n_3], tan(ω/2))"""
+                    vector=False):
+      """
+      Rodrigues-Frank vector: ([n_1, n_2, n_3], tan(ω/2)).
+
+      Parameters
+      ----------
+      vector : bool, optional
+          return as array of length 3, i.e. scale the unit vector giving the rotation axis.
+
+      """
       ro = qu2ro(self.quaternion.asArray())
       return ro[:3]*ro[3] if vector else ro
   
     def asHomochoric(self):
-      """Homochoric vector: (h_1, h_2, h_3)"""
+      """Homochoric vector: (h_1, h_2, h_3)."""
       return qu2ho(self.quaternion.asArray())
          
     def asCubochoric(self):
+      """Cubochoric vector: (c_1, c_2, c_3)."""
       return qu2cu(self.quaternion.asArray())
       
     def asM(self):
-      """Intermediate representation supporting quaternion averaging (see F. Landis Markley et al.)"""
+      """
+      Intermediate representation supporting quaternion averaging.
+
+      References
+      ----------
+      F. Landis Markley et al., Journal of Guidance, Control, and Dynamics 30(4):1193-1197, 2007
+      https://doi.org/10.2514/1.28949
+
+      """
       return self.quaternion.asM()
       
 
@@ -299,14 +367,20 @@ class Rotation:
                     rotations,
                     weights = []):
       """
-      Average rotation
+      Average rotation.
-
+      
-      ref: F. Landis Markley, Yang Cheng, John Lucas Crassidis, and Yaakov Oshman.
+      References
-           Averaging Quaternions,
+      ----------
-           Journal of Guidance, Control, and Dynamics, Vol. 30, No. 4 (2007), pp. 1193-1197.
+      F. Landis Markley et al., Journal of Guidance, Control, and Dynamics 30(4):1193-1197, 2007
-           doi: 10.2514/1.28949
+      https://doi.org/10.2514/1.28949
            
-      usage: input a list of rotations and optional weights
+      Parameters
+      ----------
+      rotations : list of Rotations
+          Rotations to average from
+      weights : list of floats, optional
+          Weights for each rotation used for averaging
+
       """
       if not all(isinstance(item, Rotation) for item in rotations):
         raise TypeError("Only instances of Rotation can be averaged.")
@@ -339,42 +413,78 @@ class Rotation:
 # ******************************************************************************************
 class Symmetry:
   """
-  Symmetry operations for lattice systems
+  Symmetry operations for lattice systems.
   
+  References
+  ----------
   https://en.wikipedia.org/wiki/Crystal_system
+  
   """
 
   lattices = [None,'orthorhombic','tetragonal','hexagonal','cubic',]
 
   def __init__(self, symmetry = None):
-    if isinstance(symmetry, str) and symmetry.lower() in Symmetry.lattices:
+    """
-      self.lattice = symmetry.lower()
+    Symmetry Definition.
-    else:
+
-      self.lattice = None
+    Parameters
+    ----------
+    symmetry : str, optional
+        label of the crystal system
+
+    """
+    if symmetry is not None and symmetry.lower() not in Symmetry.lattices:
+      raise KeyError('Symmetry/crystal system "{}" is unknown'.format(symmetry))
+
+    self.lattice = symmetry.lower() if isinstance(symmetry,str) else symmetry
 
 
   def __copy__(self):
-    """Copy"""
+    """Copy."""
     return self.__class__(self.lattice)
 
   copy = __copy__
 
 
   def __repr__(self):
-    """Readable string"""
+    """Readable string."""
     return '{}'.format(self.lattice)
 
 
   def __eq__(self, other):
-    """Equal to other"""
+    """
+    Equal to other.
+    
+    Parameters
+    ----------
+    other : Symmetry
+        Symmetry to check for equality.
+        
+    """
     return self.lattice == other.lattice
 
   def __neq__(self, other):
-    """Not equal to other"""
+    """
+    Not Equal to other.
+    
+    Parameters
+    ----------
+    other : Symmetry
+        Symmetry to check for inequality.
+        
+    """
     return not self.__eq__(other)
 
   def __cmp__(self,other):
-    """Linear ordering"""
+    """
+    Linear ordering.
+
+    Parameters
+    ----------
+    other : Symmetry
+        Symmetry to check for for order.
+        
+    """
     myOrder    = Symmetry.lattices.index(self.lattice)
     otherOrder = Symmetry.lattices.index(other.lattice)
     return (myOrder > otherOrder) - (myOrder < otherOrder)
@@ -458,7 +568,7 @@ class Symmetry:
 
   def inFZ(self,rodrigues):
     """
-    Check whether given Rodrigues vector falls into fundamental zone of own symmetry.
+    Check whether given Rodriques-Frank vector falls into fundamental zone of own symmetry.
     
     Fundamental zone in Rodrigues space is point symmetric around origin.
     """
@@ -491,11 +601,13 @@ class Symmetry:
 
   def inDisorientationSST(self,rodrigues):
     """
-    Check whether given Rodrigues vector (of misorientation) falls into standard stereographic triangle of own symmetry.
+    Check whether given Rodriques-Frank vector (of misorientation) falls into standard stereographic triangle of own symmetry.
+
+    References
+    ----------
+    A. Heinz and P. Neumann, Acta Crystallographica Section A 47:780-789, 1991
+    https://doi.org/10.1107/S0108767391006864
 
-    Determination of disorientations follow the work of A. Heinz and P. Neumann:
-    Representation of Orientation and Disorientation Data for Cubic, Hexagonal, Tetragonal and Orthorhombic Crystals
-    Acta Cryst. (1991). A47, 780-789
     """
     if (len(rodrigues) != 3):
       raise ValueError('Input is not a Rodriques-Frank vector.\n')
@@ -611,11 +723,15 @@ class Symmetry:
 # ******************************************************************************************
 class Lattice:
   """
-  Lattice system
+  Lattice system.
   
   Currently, this contains only a mapping from Bravais lattice to symmetry
   and orientation relationships. It could include twin and slip systems.
+
+  References
+  ----------
   https://en.wikipedia.org/wiki/Bravais_lattice
+
   """
 
   lattices = {
@@ -628,18 +744,27 @@ class Lattice:
 
 
   def __init__(self, lattice):
+    """
+    New lattice of given type.
+      
+    Parameters
+    ----------
+    lattice : str
+        Bravais lattice.
+
+    """
     self.lattice  = lattice
     self.symmetry = Symmetry(self.lattices[lattice]['symmetry'])
     
     
   def __repr__(self):
-    """Report basic lattice information"""
+    """Report basic lattice information."""
     return 'Bravais lattice {} ({} symmetry)'.format(self.lattice,self.symmetry)
     
     
   # Kurdjomov--Sachs orientation relationship for fcc <-> bcc transformation
-  # from S. Morito et al. Journal of Alloys and Compounds 577 (2013) 587-S592
+  # from S. Morito et al., Journal of Alloys and Compounds 577:s587-s592, 2013 
-  # also see K. Kitahara et al. Acta Materialia 54 (2006) 1279-1288
+  # also see K. Kitahara et al., Acta Materialia 54:1279-1288, 2006 
   KS = {'mapping':{'fcc':0,'bcc':1},
       'planes': np.array([
       [[  1,  1,  1],[  0,  1,  1]],
@@ -693,7 +818,7 @@ class Lattice:
       [[  1,  0,  1],[ -1,  1, -1]]],dtype='float')}
       
   # Greninger--Troiano orientation relationship for fcc <-> bcc transformation
-  # from Y. He et al. Journal of Applied Crystallography 39 (2006) 72-81
+  # from Y. He et al., Journal of Applied Crystallography 39:72-81, 2006 
   GT = {'mapping':{'fcc':0,'bcc':1},
       'planes': np.array([
       [[  1,  1,  1],[  1,  0,  1]],
@@ -747,7 +872,7 @@ class Lattice:
       [[-17,-12,  5],[-17,  7, 17]]],dtype='float')}
  
   # Greninger--Troiano' orientation relationship for fcc <-> bcc transformation
-  # from Y. He et al. Journal of Applied Crystallography 39 (2006) 72-81
+  # from Y. He et al., Journal of Applied Crystallography 39:72-81, 2006 
   GTprime = {'mapping':{'fcc':0,'bcc':1},
       'planes': np.array([
       [[  7, 17, 17],[ 12,  5, 17]],
@@ -801,7 +926,7 @@ class Lattice:
       [[  1,  1,  0],[  1,  1, -1]]],dtype='float')}
     
   # Nishiyama--Wassermann orientation relationship for fcc <-> bcc transformation
-  # from H. Kitahara et al. Materials Characterization 54 (2005) 378-386
+  # from H. Kitahara et al., Materials Characterization 54:378-386, 2005
   NW = {'mapping':{'fcc':0,'bcc':1},
       'planes': np.array([
       [[  1,  1,  1],[  0,  1,  1]],
@@ -831,7 +956,7 @@ class Lattice:
       [[ -1, -1, -2],[  0, -1,  1]]],dtype='float')}
                
   # Pitsch orientation relationship for fcc <-> bcc transformation              
-  # from Y. He et al. Acta Materialia 53 (2005) 1179-1190               
+  # from Y. He et al., Acta Materialia 53:1179-1190, 2005 
   Pitsch = {'mapping':{'fcc':0,'bcc':1},
       'planes': np.array([
       [[  0,  1,  0],[ -1,  0,  1]],
@@ -861,7 +986,7 @@ class Lattice:
       [[  1,  1,  0],[  1,  1, -1]]],dtype='float')}
                               
   # Bain orientation relationship for fcc <-> bcc transformation                        
-  # from Y. He et al./Journal of Applied Crystallography (2006). 39, 72-81
+  # from Y. He et al., Journal of Applied Crystallography 39:72-81, 2006
   Bain = {'mapping':{'fcc':0,'bcc':1},
       'planes': np.array([
       [[  1,  0,  0],[  1,  0,  0]],
@@ -873,12 +998,32 @@ class Lattice:
       [[  1,  0,  0],[  1,  1,  0]]],dtype='float')}
   
   def relationOperations(self,model):
+    """
+    Crystallographic orientation relationships for phase transformations.
     
+    References
+    ----------
+    S. Morito et al., Journal of Alloys and Compounds 577:s587-s592, 2013
+    https://doi.org/10.1016/j.jallcom.2012.02.004
+    
+    K. Kitahara et al., Acta Materialia 54(5):1279-1288, 2006
+    https://doi.org/10.1016/j.actamat.2005.11.001
+    
+    Y. He et al., Journal of Applied Crystallography 39:72-81, 2006
+    https://doi.org/10.1107/S0021889805038276
+    
+    H. Kitahara et al., Materials Characterization 54(4-5):378-386, 2005
+    https://doi.org/10.1016/j.matchar.2004.12.015
+    
+    Y. He et al., Acta Materialia 53(4):1179-1190, 2005 
+    https://doi.org/10.1016/j.actamat.2004.11.021
+
+    """
     models={'KS':self.KS, 'GT':self.GT, "GT'":self.GTprime, 
             'NW':self.NW, 'Pitsch': self.Pitsch, 'Bain':self.Bain}
     try:
       relationship = models[model]
-    except:
+    except KeyError :
       raise KeyError('Orientation relationship "{}" is unknown'.format(model))
 
     if self.lattice not in relationship['mapping']:
@@ -909,19 +1054,29 @@ class Lattice:
 
 class Orientation:
   """
-  Crystallographic orientation
+  Crystallographic orientation.
   
-  A crystallographic orientation contains a rotation and a lattice
+  A crystallographic orientation contains a rotation and a lattice.
   """
 
   __slots__ = ['rotation','lattice']
   
   def __repr__(self):
-    """Report lattice type and orientation"""
+    """Report lattice type and orientation."""
     return self.lattice.__repr__()+'\n'+self.rotation.__repr__()
 
   def __init__(self, rotation, lattice):
-  
+    """
+    New orientation from rotation and lattice.
+
+    Parameters
+    ----------
+    rotation : Rotation
+        Rotation specifying the lattice orientation.
+    lattice : Lattice
+        Lattice type of the crystal.
+
+    """
     if isinstance(lattice, Lattice):
       self.lattice = lattice
     else:
@@ -971,7 +1126,7 @@ class Orientation:
     return self.lattice.symmetry.inFZ(self.rotation.asRodrigues(vector=True))
   
   def equivalentOrientations(self,members=[]):
-    """List of orientations which are symmetrically equivalent"""
+    """List of orientations which are symmetrically equivalent."""
     try:
       iter(members)                                                                                 # asking for (even empty) list of members?
     except TypeError:
@@ -981,12 +1136,12 @@ class Orientation:
                                     for q in self.lattice.symmetry.symmetryOperations(members)]           # yes, return list of rotations
     
   def relatedOrientations(self,model):
-    """List of orientations related by the given orientation relationship"""
+    """List of orientations related by the given orientation relationship."""
     r = self.lattice.relationOperations(model)
     return [self.__class__(self.rotation*o,r['lattice']) for o in r['rotations']]
     
   def reduced(self):
-    """Transform orientation to fall into fundamental zone according to symmetry"""
+    """Transform orientation to fall into fundamental zone according to symmetry."""
     for me in self.equivalentOrientations():
       if self.lattice.symmetry.inFZ(me.rotation.asRodrigues(vector=True)): break
 
@@ -996,7 +1151,7 @@ class Orientation:
                   axis,
                   proper = False,
                   SST = True):
-    """Axis rotated according to orientation (using crystal symmetry to ensure location falls into SST)"""
+    """Axis rotated according to orientation (using crystal symmetry to ensure location falls into SST)."""
     if SST:                                                                                         # pole requested to be within SST
       for i,o in enumerate(self.equivalentOrientations()):                                          # test all symmetric equivalent quaternions
         pole = o.rotation*axis                                                                      # align crystal direction to axis
@@ -1008,7 +1163,7 @@ class Orientation:
     
     
   def IPFcolor(self,axis):
-    """TSL color of inverse pole figure for given axis"""
+    """TSL color of inverse pole figure for given axis."""
     color = np.zeros(3,'d')
 
     for o in self.equivalentOrientations():
@@ -1023,7 +1178,7 @@ class Orientation:
   def fromAverage(cls,
                   orientations,
                   weights = []):
-    """Create orientation from average of list of orientations"""
+    """Create orientation from average of list of orientations."""
     if not all(isinstance(item, Orientation) for item in orientations):
       raise TypeError("Only instances of Orientation can be averaged.")
 
@@ -1039,7 +1194,7 @@ class Orientation:
 
 
   def average(self,other):
-    """Calculate the average rotation"""
+    """Calculate the average rotation."""
     return Orientation.fromAverage([self,other])
 
 
@@ -1080,10 +1235,10 @@ def isone(a):
 def iszero(a):
   return np.isclose(a,0.0,atol=1.0e-12,rtol=0.0)
   
-#---------- quaternion ----------
+#---------- Quaternion ----------
 
 def qu2om(qu):
-  """Quaternion to orientation matrix"""
+  """Quaternion to rotation matrix."""
   qq = qu[0]**2-(qu[1]**2 + qu[2]**2 + qu[3]**2)
   om = np.diag(qq + 2.0*np.array([qu[1],qu[2],qu[3]])**2)
 
@@ -1097,7 +1252,7 @@ def qu2om(qu):
   
   
 def qu2eu(qu):
-  """Quaternion to Euler angles""" 
+  """Quaternion to Bunge-Euler angles.""" 
   q03 = qu[0]**2+qu[3]**2
   q12 = qu[1]**2+qu[2]**2
   chi = np.sqrt(q03*q12)
@@ -1117,7 +1272,7 @@ def qu2eu(qu):
 
 def qu2ax(qu):
   """
-  Quaternion to axis angle
+  Quaternion to axis angle pair.
   
   Modified version of the original formulation, should be numerically more stable
   """
@@ -1134,7 +1289,7 @@ def qu2ax(qu):
 
 
 def qu2ro(qu):
-  """Quaternion to Rodrigues vector"""
+  """Quaternion to Rodriques-Frank vector."""
   if iszero(qu[0]):
     ro = [qu[1], qu[2], qu[3], np.inf]
   else:
@@ -1146,7 +1301,7 @@ def qu2ro(qu):
 
 
 def qu2ho(qu):
-  """Quaternion to homochoric"""
+  """Quaternion to homochoric vector."""
   omega = 2.0 * np.arccos(np.clip(qu[0],-1.0,1.0))                                                  # avoid numerical difficulties
   
   if iszero(omega):
@@ -1160,15 +1315,15 @@ def qu2ho(qu):
 
 
 def qu2cu(qu):
-  """Quaternion to cubochoric"""
+  """Quaternion to cubochoric vector."""
   return ho2cu(qu2ho(qu))
 
 
-#---------- orientation matrix ----------
+#---------- Rotation matrix ----------
 
 def om2qu(om):
   """
-  Orientation matrix to quaternion
+  Rotation matrix to quaternion.
   
   The original formulation (direct conversion) had (numerical?) issues
   """
@@ -1176,7 +1331,7 @@ def om2qu(om):
 
 
 def om2eu(om):
-  """Orientation matrix to Euler angles"""
+  """Rotation matrix to Bunge-Euler angles."""
   if abs(om[2,2]) < 1.0:
     zeta = 1.0/np.sqrt(1.0-om[2,2]**2)
     eu = np.array([np.arctan2(om[2,0]*zeta,-om[2,1]*zeta),
@@ -1191,7 +1346,7 @@ def om2eu(om):
 
 
 def om2ax(om):
-  """Orientation matrix to axis angle"""
+  """Rotation matrix to axis angle pair."""
   ax=np.empty(4)
 
   # first get the rotation angle
@@ -1212,24 +1367,24 @@ def om2ax(om):
   
 
 def om2ro(om):
-  """Orientation matrix to Rodriques vector"""
+  """Rotation matrix to Rodriques-Frank vector."""
   return eu2ro(om2eu(om))
   
 
 def om2ho(om):
-  """Orientation matrix to homochoric"""
+  """Rotation matrix to homochoric vector."""
   return ax2ho(om2ax(om))
 
 
 def om2cu(om):
-  """Orientation matrix to cubochoric"""
+  """Rotation matrix to cubochoric vector."""
   return ho2cu(om2ho(om))
   
   
-#---------- Euler angles ----------
+#---------- Bunge-Euler angles ----------
 
 def eu2qu(eu):
-  """Euler angles to quaternion"""
+  """Bunge-Euler angles to quaternion."""
   ee = 0.5*eu
   cPhi = np.cos(ee[1])
   sPhi = np.sin(ee[1])
@@ -1242,7 +1397,7 @@ def eu2qu(eu):
 
 
 def eu2om(eu):
-  """Euler angles to orientation matrix"""
+  """Bunge-Euler angles to rotation matrix."""
   c = np.cos(eu)
   s = np.sin(eu)
 
@@ -1255,7 +1410,7 @@ def eu2om(eu):
 
 
 def eu2ax(eu):
-  """Euler angles to axis angle""" 
+  """Bunge-Euler angles to axis angle pair.""" 
   t = np.tan(eu[1]*0.5)
   sigma = 0.5*(eu[0]+eu[2])
   delta = 0.5*(eu[0]-eu[2])
@@ -1266,7 +1421,7 @@ def eu2ax(eu):
   if iszero(alpha):
     ax = np.array([ 0.0, 0.0, 1.0, 0.0 ])
   else:
-    ax = -P/tau * np.array([ t*np.cos(delta), t*np.sin(delta), np.sin(sigma) ])                     # passive axis-angle pair so a minus sign in front
+    ax = -P/tau * np.array([ t*np.cos(delta), t*np.sin(delta), np.sin(sigma) ])                     # passive axis angle pair so a minus sign in front
     ax = np.append(ax,alpha) 
     if alpha < 0.0: ax *= -1.0                                                                      # ensure alpha is positive
   
@@ -1274,8 +1429,8 @@ def eu2ax(eu):
 
 
 def eu2ro(eu):
-  """Euler angles to Rodrigues vector"""
+  """Bunge-Euler angles to Rodriques-Frank vector."""
-  ro = eu2ax(eu)                                                                                    # convert to axis angle representation
+  ro = eu2ax(eu)                                                                                    # convert to axis angle pair representation
   if ro[3] >= np.pi:                                                                                # Differs from original implementation. check convention 5
     ro[3] = np.inf
   elif iszero(ro[3]):
@@ -1287,19 +1442,19 @@ def eu2ro(eu):
 
 
 def eu2ho(eu):
-  """Euler angles to homochoric"""
+  """Bunge-Euler angles to homochoric vector."""
   return ax2ho(eu2ax(eu))
 
 
 def eu2cu(eu):
-  """Euler angles to cubochoric"""
+  """Bunge-Euler angles to cubochoric vector."""
   return ho2cu(eu2ho(eu))
 
 
-#---------- axis angle ----------
+#---------- Axis angle pair ----------
 
 def ax2qu(ax):
-  """Axis angle to quaternion"""
+  """Axis angle pair to quaternion."""
   if iszero(ax[3]):
     qu = np.array([ 1.0, 0.0, 0.0, 0.0 ])
   else:
@@ -1311,7 +1466,7 @@ def ax2qu(ax):
 
 
 def ax2om(ax):
-  """Axis angle to orientation matrix"""
+  """Axis angle pair to rotation matrix."""
   c = np.cos(ax[3])
   s = np.sin(ax[3])
   omc = 1.0-c
@@ -1326,12 +1481,12 @@ def ax2om(ax):
 
   
 def ax2eu(ax):
-  """Orientation matrix to Euler angles"""
+  """Rotation matrix to Bunge Euler angles."""
   return om2eu(ax2om(ax))
 
 
 def ax2ro(ax):
-  """Axis angle to Rodrigues vector""" 
+  """Axis angle pair to Rodriques-Frank vector.""" 
   if iszero(ax[3]):
     ro = [ 0.0, 0.0, P, 0.0 ]
   else:
@@ -1344,36 +1499,36 @@ def ax2ro(ax):
 
 
 def ax2ho(ax):
-  """Axis angle to homochoric""" 
+  """Axis angle pair to homochoric vector.""" 
   f = (0.75 * ( ax[3] - np.sin(ax[3]) ))**(1.0/3.0)
   ho = ax[0:3] * f
   return ho
 
 
 def ax2cu(ax):
-  """Axis angle to cubochoric"""
+  """Axis angle pair to cubochoric vector."""
   return ho2cu(ax2ho(ax))
 
 
-#---------- Rodrigues--Frank ----------
+#---------- Rodrigues-Frank vector ----------
 
 def ro2qu(ro):
-  """Rodrigues vector to quaternion"""
+  """Rodriques-Frank vector to quaternion."""
   return ax2qu(ro2ax(ro))
 
 
 def ro2om(ro):
- """Rodgrigues vector to orientation matrix"""
+ """Rodgrigues-Frank vector to rotation matrix."""
  return ax2om(ro2ax(ro))
 
 
 def ro2eu(ro):
-  """Rodrigues vector to orientation matrix"""
+  """Rodriques-Frank vector to Bunge-Euler angles."""
   return om2eu(ro2om(ro))
  
  
 def ro2ax(ro):
-  """Rodrigues vector to axis angle"""
+  """Rodriques-Frank vector to axis angle pair."""
   ta = ro[3]
   
   if iszero(ta):
@@ -1389,7 +1544,7 @@ def ro2ax(ro):
 
 
 def ro2ho(ro):
-  """Rodrigues vector to homochoric""" 
+  """Rodriques-Frank vector to homochoric vector.""" 
   if iszero(np.sum(ro[0:3]**2.0)):
     ho = [ 0.0, 0.0, 0.0 ]
   else:
@@ -1400,29 +1555,29 @@ def ro2ho(ro):
 
 
 def ro2cu(ro):
-  """Rodrigues vector to cubochoric"""
+  """Rodriques-Frank vector to cubochoric vector."""
   return ho2cu(ro2ho(ro))
 
 
-#---------- homochoric ----------
+#---------- Homochoric vector----------
 
 def ho2qu(ho):
-  """Homochoric to quaternion"""
+  """Homochoric vector to quaternion."""
   return ax2qu(ho2ax(ho))
 
 
 def ho2om(ho):
-  """Homochoric to orientation matrix"""
+  """Homochoric vector to rotation matrix."""
   return ax2om(ho2ax(ho))
 
 
 def ho2eu(ho):
-  """Homochoric to Euler angles"""
+  """Homochoric vector to Bunge-Euler angles."""
   return ax2eu(ho2ax(ho))
   
 
 def ho2ax(ho):
-  """Homochoric to axis angle"""
+  """Homochoric vector to axis angle pair."""
   tfit = np.array([+1.0000000000018852,      -0.5000000002194847,
                    -0.024999992127593126,    -0.003928701544781374,
                    -0.0008152701535450438,   -0.0002009500426119712,
@@ -1448,42 +1603,42 @@ def ho2ax(ho):
 
 
 def ho2ro(ho):
-  """Axis angle to Rodriques vector"""
+  """Axis angle pair to Rodriques-Frank vector."""
   return ax2ro(ho2ax(ho))
 
 
 def ho2cu(ho):
-  """Homochoric to cubochoric"""
+  """Homochoric vector to cubochoric vector."""
   return  Lambert.BallToCube(ho)
 
 
-#---------- cubochoric ----------
+#---------- Cubochoric ----------
 
 def cu2qu(cu):
-  """Cubochoric to quaternion"""
+  """Cubochoric vector to quaternion."""
   return ho2qu(cu2ho(cu))
 
 
 def cu2om(cu):
-  """Cubochoric to orientation matrix"""
+  """Cubochoric vector to rotation matrix."""
   return ho2om(cu2ho(cu))
 
 
 def cu2eu(cu):
-  """Cubochoric to Euler angles"""
+  """Cubochoric vector to Bunge-Euler angles."""
   return ho2eu(cu2ho(cu))
 
 
 def cu2ax(cu):
-  """Cubochoric to axis angle"""
+  """Cubochoric vector to axis angle pair."""
   return ho2ax(cu2ho(cu))
 
 
 def cu2ro(cu):
-  """Cubochoric to Rodrigues vector"""
+  """Cubochoric vector to Rodriques-Frank vector."""
   return ho2ro(cu2ho(cu))
 
 
 def cu2ho(cu):
-  """Cubochoric to homochoric"""
+  """Cubochoric vector to homochoric vector."""
   return  Lambert.CubeToBall(cu)

python/damask/solver/__init__.py

@@ -1,7 +1,5 @@
-# -*- coding: UTF-8 no BOM -*-
+"""Tools to control the various solvers."""
-"""Tools to control the various BVP solvers"""
 
 from .solver   import Solver     # noqa
-from .spectral import Spectral   # noqa
 from .marc     import Marc       # noqa
 from .abaqus   import Abaqus     # noqa

python/damask/solver/abaqus.py

@@ -1,17 +1,24 @@
-# -*- coding: UTF-8 no BOM -*-
+import subprocess
 
 from .solver import Solver
 import damask
-import subprocess
+
 
 class Abaqus(Solver):
+  """Wrapper to run DAMASK with Abaqus."""
 
-  def __init__(self,version=''):                                                                    # example version string: 2017
+  def __init__(self,version=int(damask.Environment().options['ABAQUS_VERSION'])):
-    self.solver='Abaqus'
+    """
-    if version =='':
+    Create a Abaqus solver object.
-      version = damask.Environment().options['ABAQUS_VERSION']
+    
-    else:
+    Parameters
-      self.version = version
+    ----------
+    version : integer
+        Abaqus version
+
+    """
+    self.solver  ='Abaqus'
+    self.version = int(version)
 
   def return_run_command(self,model):
     env=damask.Environment()
@@ -21,7 +28,7 @@ class Abaqus(Solver):
     except OSError:                                                                                 # link to abqXXX not existing
       cmd='abaqus'
       process = subprocess.Popen(['abaqus','information=release'],stdout = subprocess.PIPE,stderr = subprocess.PIPE)
-      detectedVersion = process.stdout.readlines()[1].split()[1].decode('utf-8')
+      detectedVersion = int(process.stdout.readlines()[1].split()[1].decode('utf-8'))
       if self.version != detectedVersion:
         raise Exception('found Abaqus version {}, but requested {}'.format(detectedVersion,self.version))
     return '{} -job {} -user {}/src/DAMASK_abaqus interactive'.format(cmd,model,env.rootDir())

python/damask/solver/marc.py

@@ -1,49 +1,47 @@
-# -*- coding: UTF-8 no BOM -*-
+import os
-
+import subprocess
+import shlex
 
 from .solver import Solver
-
+import damask
 
 class Marc(Solver):
+  """Wrapper to run DAMASK with MSCMarc."""
 
-  def __init__(self):
+  def __init__(self,version=float(damask.Environment().options['MARC_VERSION'])):
-    self.solver = 'Marc'
+    """
-
+    Create a Marc solver object.
-
-#--------------------------
-  def version(self):
-    import damask.environment
-
-    return damask.environment.Environment().options['MARC_VERSION']
-      
     
-#--------------------------
+    Parameters
-  def libraryPath(self,release = ''):
+    ----------
-    import os,damask.environment
+    version : float
+        Marc version
 
-    MSCpath     = damask.environment.Environment().options['MSC_ROOT']
+    """
-    if len(release) == 0: release = self.version()
+    self.solver  ='Marc'
-    
+    self.version = float(damask.environment.Environment().options['MARC_VERSION'])
-    path = '{}/mentat{}/shlib/linux64'.format(MSCpath,release)
-
-    return path if os.path.exists(path) else ''
 
 
 #--------------------------
-  def toolsPath(self,release = ''):
+  def libraryPath(self):
-    import os,damask.environment
 
-    MSCpath = damask.environment.Environment().options['MSC_ROOT']
+    path_MSC = damask.environment.Environment().options['MSC_ROOT']   
-    if len(release) == 0: release = self.version()
+    path_lib = '{}/mentat{}/shlib/linux64'.format(path_MSC,self.version)
 
-    path = '%s/marc%s/tools'%(MSCpath,release)
+    return path_lib if os.path.exists(path_lib) else ''
+
+
+#--------------------------
+  def toolsPath(self):
+
+    path_MSC   = damask.environment.Environment().options['MSC_ROOT']
+    path_tools = '{}/marc{}/tools'.format(path_MSC,self.version)
+
+    return path_tools if os.path.exists(path_tools) else ''
 
-    return path if os.path.exists(path) else ''
-  
 
 #--------------------------
   def submit_job(self,
-                 release      = '',
                  model        = 'model',
                  job          = 'job1',
                  logfile      = None,
@@ -51,25 +49,22 @@ class Marc(Solver):
                  optimization ='',
                 ):
 
-    import os,damask.environment
+
-    import subprocess,shlex
-    
-    if len(release) == 0: release = self.version()
     damaskEnv = damask.environment.Environment()
     
     user = 'not found'
 
     if compile:
-      if os.path.isfile(os.path.join(damaskEnv.relPath('src/'),'DAMASK_marc{}.f90'.format(release))):
+      if os.path.isfile(os.path.join(damaskEnv.relPath('src/'),'DAMASK_marc{}.f90'.format(self.version))):
-        user = os.path.join(damaskEnv.relPath('src/'),'DAMASK_marc{}'.format(release))
+        user = os.path.join(damaskEnv.relPath('src/'),'DAMASK_marc{}'.format(self.version))
     else:
-      if os.path.isfile(os.path.join(damaskEnv.relPath('src/'),'DAMASK_marc{}.marc'.format(release))):
+      if os.path.isfile(os.path.join(damaskEnv.relPath('src/'),'DAMASK_marc{}.marc'.format(self.version))):
-        user = os.path.join(damaskEnv.relPath('src/'),'DAMASK_marc{}'.format(release))
+        user = os.path.join(damaskEnv.relPath('src/'),'DAMASK_marc{}'.format(self.version))
 
     # Define options [see Marc Installation and Operation Guide, pp 23]
     script = 'run_damask_{}mp'.format(optimization)
     
-    cmd = os.path.join(self.toolsPath(release),script) + \
+    cmd = os.path.join(self.toolsPath(),script) + \
           ' -jid ' + model + '_' + job + \
           ' -nprocd 1  -autorst 0 -ci n  -cr n  -dcoup 0 -b no -v no'
 

python/damask/solver/solver.py

@@ -1,21 +1,6 @@
-# -*- coding: UTF-8 no BOM -*-
-
-
-import damask.solver
-
 class Solver():
   """
   General class for solver specific functionality.
 
   Sub-classed by the individual solvers.
   """
-  
-  def __init__(self,solver=''):
-    solverClass = {
-                      'spectral': damask.solver.Spectral,
-                      'marc':     damask.solver.Marc,
-                    }
-    if solver.lower() in list(solverClass.keys()):
-      self.__class__=solverClass[solver.lower()]
-      self.__init__()
-

python/damask/solver/spectral.py

@@ -1,9 +0,0 @@
-# -*- coding: UTF-8 no BOM -*-
-
-
-from .solver import Solver
-
-class Spectral(Solver):
-
-  def __init__(self):
-    self.solver='Spectral'

src/crystallite.f90

@@ -1274,7 +1274,7 @@ logical function integrateStress(ipc,ip,el,timeFraction)
       if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0 &
           .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) &
                  .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0)) then
-        write(6,'(a,i3,/)')                   '<< CRYST integrateStress >> iteration ', NiterationStressLp
+        write(6,'(a,i3,/)')                   '<< CRYST integrateStress >> Lp iteration ', NiterationStressLp
         write(6,'(a,/,3(12x,3(e20.10,1x)/))') '<< CRYST integrateStress >> Lpguess', transpose(Lpguess)
         write(6,'(a,/,3(12x,3(e20.10,1x)/))') '<< CRYST integrateStress >> Lp_constitutive', transpose(Lp_constitutive)
         write(6,'(a,/,3(12x,3(e20.10,1x)/))') '<< CRYST integrateStress >> Fi', transpose(Fi_new)
@@ -1377,6 +1377,7 @@ logical function integrateStress(ipc,ip,el,timeFraction)
       if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0 &
           .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) &
                  .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0)) then
+        write(6,'(a,i3,/)')                  '<< CRYST integrateStress >> Li iteration ', NiterationStressLi
         write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST integrateStress >> Li_constitutive', transpose(Li_constitutive)
         write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST integrateStress >> Liguess', transpose(Liguess)
       endif
@@ -1405,6 +1406,13 @@ logical function integrateStress(ipc,ip,el,timeFraction)
     else                                                                                            ! not converged and residuum not improved...
       steplengthLi   = num%subStepSizeLi * steplengthLi                                             ! ...try with smaller step length in same direction
       Liguess        = Liguess_old + steplengthLi * deltaLi
+#ifdef DEBUG
+        if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0 &
+            .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) &
+                   .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0)) then
+          write(6,'(a,1x,f7.4)') '<< CRYST integrateStress >> linear search for Liguess with step', steplengthLi
+        endif
+#endif
       cycle LiLoop
     endif
  
@@ -1450,6 +1458,13 @@ logical function integrateStress(ipc,ip,el,timeFraction)
     jacoCounterLi = jacoCounterLi + 1
  
     Liguess = Liguess + steplengthLi * deltaLi
+#ifdef DEBUG
+        if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0 &
+            .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) &
+                   .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0)) then
+          write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST integrateStress >> corrected Liguess by', transpose(deltaLi)
+        endif
+#endif
   enddo LiLoop
  
   !* calculate new plastic and elastic deformation gradient

src/grid/DAMASK_grid.f90

@@ -118,6 +118,7 @@ program DAMASK_spectral
 ! assign mechanics solver depending on selected type
  select case (trim(config_numerics%getString('spectral_solver',defaultVal='basic')))
    case ('basic')
+     write(6,'(a)')   ' Basic solution algorithm'
      mech_init     => grid_mech_spectral_basic_init
      mech_forward  => grid_mech_spectral_basic_forward
      mech_solution => grid_mech_spectral_basic_solution
@@ -125,6 +126,7 @@ program DAMASK_spectral
    case ('polarisation')
      if(iand(debug_level(debug_spectral),debug_levelBasic)/= 0) &
        call IO_warning(42, ext_msg='debug Divergence')
+     write(6,'(a)')   ' Polarisation solution algorithm'
      mech_init     => grid_mech_spectral_polarisation_init
      mech_forward  => grid_mech_spectral_polarisation_forward
      mech_solution => grid_mech_spectral_polarisation_solution
@@ -132,6 +134,7 @@ program DAMASK_spectral
    case ('fem')
      if(iand(debug_level(debug_spectral),debug_levelBasic)/= 0) &
        call IO_warning(42, ext_msg='debug Divergence')
+     write(6,'(a)')   ' FEM solution algorithm'
      mech_init     => grid_mech_FEM_init
      mech_forward  => grid_mech_FEM_forward
      mech_solution => grid_mech_FEM_solution

src/plastic_isotropic.f90

@@ -283,49 +283,52 @@ end subroutine plastic_isotropic_LpAndItsTangent
 
 !--------------------------------------------------------------------------------------------------
 !> @brief calculates plastic velocity gradient and its tangent
-! ToDo: Rename Tstar to Mi?
 !--------------------------------------------------------------------------------------------------
-subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dTstar,Tstar,instance,of)
+subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,instance,of)
  
   real(pReal), dimension(3,3), intent(out) :: &
     Li                                                                                              !< inleastic velocity gradient
   real(pReal), dimension(3,3,3,3), intent(out)  :: &
-    dLi_dTstar                                                                                      !< derivative of Li with respect to the Mandel stress
+    dLi_dMi                                                                                         !< derivative of Li with respect to the Mandel stress
  
   real(pReal), dimension(3,3),   intent(in) :: &
-    Tstar                                                                                           !< Mandel stress ToDo: Mi?
+    Mi                                                                                              !< Mandel stress 
   integer,                       intent(in) :: &
     instance, &
     of
  
   real(pReal), dimension(3,3) :: &
-    Tstar_sph                                                                                       !< sphiatoric part of the Mandel stress
+    Mi_sph                                                                                          !< spherical part of the Mandel stress
   real(pReal) :: &
-    dot_gamma, &                                                                                    !< strainrate
+    dot_gamma, &                                                                                    !< shear rate
-    norm_Tstar_sph, &                                                                               !< euclidean norm of Tstar_sph
+    tr                                                                                               !< pressure
-    squarenorm_Tstar_sph                                                                            !< square of the euclidean norm of Tstar_sph
   integer :: &
     k, l, m, n
  
   associate(prm => param(instance), stt => state(instance))
  
-  Tstar_sph = math_spherical33(Tstar)
+  tr=math_trace33(math_spherical33(Mi))
-  squarenorm_Tstar_sph = math_mul33xx33(Tstar_sph,Tstar_sph)
+
-  norm_Tstar_sph = sqrt(squarenorm_Tstar_sph)
+  if (prm%dilatation .and. abs(tr) > 0.0_pReal) then                                                 ! no stress or J2 plasticity --> Li and its derivative are zero
+    Li = math_I3 &
+       * prm%dot_gamma_0/prm%M * (3.0_pReal*prm%M*stt%xi(of))**(-prm%n) &
+       * tr * abs(tr)**(prm%n-1.0_pReal)
  
-  if (prm%dilatation .and. norm_Tstar_sph > 0.0_pReal) then                                         ! no stress or J2 plastitiy --> Li and its derivative are zero
+#ifdef DEBUG
-    dot_gamma = prm%dot_gamma_0 * (sqrt(1.5_pReal) * norm_Tstar_sph /(prm%M*stt%xi(of))) **prm%n
+    if (iand(debug_level(debug_constitutive), debug_levelExtensive) /= 0 &
- 
+        .and. (of == prm%of_debug .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0)) then
-    Li = math_I3/sqrt(3.0_pReal) * dot_gamma/prm%M
+      write(6,'(/,a,/,f12.5)') '<< CONST isotropic >> pressure / MPa', tr/3.0_pReal*1.0e-6_pReal
+      write(6,'(/,a,/,f12.5)') '<< CONST isotropic >> gdot', prm%dot_gamma_0 * (3.0_pReal*prm%M*stt%xi(of))**(-prm%n) &
+                                                           * tr * abs(tr)**(prm%n-1.0_pReal)
+    end if
+#endif
+
     forall (k=1:3,l=1:3,m=1:3,n=1:3) &
-      dLi_dTstar(k,l,m,n) = (prm%n-1.0_pReal) * Tstar_sph(k,l)*Tstar_sph(m,n) / squarenorm_Tstar_sph
+      dLi_dMi(k,l,m,n) = n / tr * Li(k,l) * math_I3(m,n)
-    forall (k=1:3,l=1:3) &
+
-      dLi_dTstar(k,l,k,l) = dLi_dTstar(k,l,k,l) + 1.0_pReal
- 
-    dLi_dTstar = dot_gamma / prm%M * dLi_dTstar / norm_Tstar_sph
   else
-    Li = 0.0_pReal
+    Li      = 0.0_pReal
-    dLi_dTstar = 0.0_pReal
+    dLi_dMi = 0.0_pReal
   endif
  
   end associate

src/plastic_phenopowerlaw.f90

@@ -250,6 +250,7 @@ subroutine plastic_phenopowerlaw_init
 !--------------------------------------------------------------------------------------------------
 ! slip-twin related parameters
    slipAndTwinActive: if (prm%totalNslip > 0 .and. prm%totalNtwin > 0) then
+     prm%h0_TwinSlip          = config%getFloat('h0_twinslip')
      prm%interaction_SlipTwin = lattice_interaction_SlipByTwin(prm%Nslip,prm%Ntwin,&
                                                                config%getFloats('interaction_sliptwin'), &
                                                                config%getString('lattice_structure'))