Merge branch 'DADF5_point_calculations-2' into development

PRIVATE

@@ -1 +1 @@
-Subproject commit ec615d249d39e5d01446b01ab9a5b7e7601340ad
+Subproject commit 6db5f4666fc651b4de3b44ceaed3f2b848170ac9

processing/post/addCauchy.py

@@ -42,8 +42,8 @@ for name in filenames:
 
     table = damask.Table.from_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name)
     table.add('Cauchy',
-              damask.mechanics.Cauchy(table.get(options.defgrad).reshape(-1,3,3),
+              damask.mechanics.Cauchy(table.get(options.stress ).reshape(-1,3,3),
-                                      table.get(options.stress ).reshape(-1,3,3)).reshape(-1,9),
+                                      table.get(options.defgrad).reshape(-1,3,3)).reshape(-1,9),
               scriptID+' '+' '.join(sys.argv[1:]))
 
     table.to_ASCII(sys.stdout if name is None else name)

processing/post/addPK2.py

@@ -43,8 +43,8 @@ for name in filenames:
     table = damask.Table.from_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name)
 
     table.add('S',
-              damask.mechanics.PK2(table.get(options.defgrad).reshape(-1,3,3),
+              damask.mechanics.PK2(table.get(options.stress ).reshape(-1,3,3),
-                                   table.get(options.stress ).reshape(-1,3,3)).reshape(-1,9),
+                                   table.get(options.defgrad).reshape(-1,3,3)).reshape(-1,9),
               scriptID+' '+' '.join(sys.argv[1:]))
 
     table.to_ASCII(sys.stdout if name is None else name)

processing/post/addSpectralDecomposition.py

@@ -2,6 +2,7 @@
 
 import os
 import sys
+from io import StringIO
 from optparse import OptionParser
 
 import numpy as np
@@ -33,69 +34,27 @@ parser.add_option('--no-check',
 
 parser.set_defaults(rh = True,
                     )
+
 (options,filenames) = parser.parse_args()
-
-if options.tensor is None:
-  parser.error('no data column specified.')
-
-# --- loop over input files -------------------------------------------------------------------------
-
 if filenames == []: filenames = [None]
 
 for name in filenames:
-  try:
+    damask.util.report(scriptName,name)
-    table = damask.ASCIItable(name = name,
-                              buffered = False)
-  except: continue
-  damask.util.report(scriptName,name)
 
-# ------------------------------------------ read header ------------------------------------------
+    table = damask.Table.from_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name)
 
-  table.head_read()
+    for tensor in options.tensor:
 
+        t = table.get(tensor).reshape(-1,3,3)
+        (u,v) = np.linalg.eigh(damask.mechanics.symmetric(t))
+        if options.rh: v[np.linalg.det(v) < 0.0,:,2] *= -1.0
 
-# ------------------------------------------ assemble header 1 ------------------------------------
+        for i,o in enumerate(['Min','Mid','Max']):
+            table.add('eigval{}({})'.format(o,tensor),u[:,i],
+                      scriptID+' '+' '.join(sys.argv[1:]))
 
-  items = {
+        for i,o in enumerate(['Min','Mid','Max']):
-            'tensor': {'dim': 9, 'shape': [3,3], 'labels':options.tensor, 'column': []},
+            table.add('eigvec{}({})'.format(o,tensor),v[:,:,i],
-          }
+                      scriptID+' '+' '.join(sys.argv[1:]))
-  errors  = []
+        
-  remarks = []
+    table.to_ASCII(sys.stdout if name is None else name)
-  
-  for type, data in items.items():
-    for what in data['labels']:
-      dim = table.label_dimension(what)
-      if dim != data['dim']: remarks.append('column {} is not a {}...'.format(what,type))
-      else:
-        items[type]['column'].append(table.label_index(what))
-        for order in ['Min','Mid','Max']:
-          table.labels_append(['eigval{}({})'.format(order,what)])                                         # extend ASCII header with new labels
-        for order in ['Min','Mid','Max']:
-          table.labels_append(['{}_eigvec{}({})'.format(i+1,order,what) for i in range(3)])                # extend ASCII header with new labels
-
-  if remarks != []: damask.util.croak(remarks)
-  if errors  != []:
-    damask.util.croak(errors)
-    table.close(dismiss = True)
-    continue
-
-# ------------------------------------------ assemble header 2 ------------------------------------
-
-  table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:]))
-  table.head_write()
-
-# ------------------------------------------ process data -----------------------------------------
-
-  outputAlive = True
-  while outputAlive and table.data_read():                                                          # read next data line of ASCII table
-    for type, data in items.items():
-      for column in data['column']:
-        (u,v) = np.linalg.eigh(np.array(list(map(float,table.data[column:column+data['dim']]))).reshape(data['shape']))
-        if options.rh and np.dot(np.cross(v[:,0], v[:,1]), v[:,2]) < 0.0 : v[:, 2] *= -1.0          # ensure right-handed eigenvector basis
-        table.data_append(list(u))                                                                  # vector of max,mid,min eigval
-        table.data_append(list(v.transpose().reshape(data['dim'])))                                 # 3x3=9 combo vector of max,mid,min eigvec coordinates
-    outputAlive = table.data_write()                                                                # output processed line in accordance with column labeling
-
-# ------------------------------------------ output finalization -----------------------------------
-
-  table.close()                                                                                     # close input ASCII table (works for stdin)

python/damask/__init__.py

@@ -15,6 +15,7 @@ from .config      import Material         # noqa
 from .colormaps   import Colormap, Color  # noqa
 from .orientation import Symmetry, Lattice, Rotation, Orientation # noqa
 from .dadf5       import DADF5            # noqa
+from .dadf5       import DADF5 as Result  # noqa
 
 from .geom        import Geom             # noqa
 from .solver      import Solver           # noqa

python/damask/environment.py

@@ -8,7 +8,7 @@ class Environment():
     def __init__(self):
         """Read and provide values of DAMASK configuration."""
         self.options = {}
-        self.get_options()
+        self.__get_options()
 
     def relPath(self,relative = '.'):
         return os.path.join(self.rootDir(),relative)
@@ -16,7 +16,7 @@ class Environment():
     def rootDir(self):
         return os.path.normpath(os.path.join(os.path.realpath(__file__),'../../../'))
 
-    def get_options(self):
+    def __get_options(self):
         for item in ['DAMASK_NUM_THREADS',
                      'MSC_ROOT',
                      'MARC_VERSION',

python/damask/mechanics.py

@@ -1,11 +1,11 @@
 import numpy as np
 
-def Cauchy(F,P):
+def Cauchy(P,F):
     """
-    Return Cauchy stress calculated from 1. Piola-Kirchhoff stress and deformation gradient.
+    Return Cauchy stress calculated from first Piola-Kirchhoff stress and deformation gradient.
-      
+
     Resulting tensor is symmetrized as the Cauchy stress needs to be symmetric.
-    
+
     Parameters
     ----------
     F : numpy.array of shape (:,3,3) or (3,3)
@@ -21,67 +21,10 @@ def Cauchy(F,P):
     return symmetric(sigma)
 
 
-def PK2(F,P):
-    """
-    Return 2. Piola-Kirchhoff stress calculated from 1. Piola-Kirchhoff stress and deformation gradient.
-      
-    Parameters
-    ----------
-    F : numpy.array of shape (:,3,3) or (3,3)
-      Deformation gradient.
-    P : numpy.array of shape (:,3,3) or (3,3)
-      1. Piola-Kirchhoff stress.
-
-    """
-    if np.shape(F) == np.shape(P) == (3,3):
-        S = np.dot(np.linalg.inv(F),P)
-    else:
-        S = np.einsum('ijk,ikl->ijl',np.linalg.inv(F),P)
-    return S
-
-
-def strain_tensor(F,t,m):
-    """
-    Return strain tensor calculated from deformation gradient.
-      
-    For details refer to https://en.wikipedia.org/wiki/Finite_strain_theory and
-    https://de.wikipedia.org/wiki/Verzerrungstensor
-    
-    Parameters
-    ----------
-    F : numpy.array of shape (:,3,3) or (3,3)
-      Deformation gradient.
-    t : {‘V’, ‘U’}
-      Type of the polar decomposition, ‘V’ for left stretch tensor and ‘U’ for right stretch tensor.
-    m : float
-      Order of the strain.
-
-    """
-    F_ = F.reshape((1,3,3)) if F.shape == (3,3) else F
-    if   t == 'V':
-        B   = np.matmul(F_,transpose(F_))
-        w,n = np.linalg.eigh(B)
-    elif t == 'U':
-        C   = np.matmul(transpose(F_),F_)
-        w,n = np.linalg.eigh(C)
-    
-    if   m > 0.0:
-        eps = 1.0/(2.0*abs(m)) * (+ np.matmul(n,np.einsum('ij,ikj->ijk',w**m,n)) 
-                                  - np.broadcast_to(np.eye(3),[F_.shape[0],3,3]))
-    elif m < 0.0:
-        eps = 1.0/(2.0*abs(m)) * (- np.matmul(n,np.einsum('ij,ikj->ijk',w**m,n))
-                                  + np.broadcast_to(np.eye(3),[F_.shape[0],3,3]))
-    else:
-        eps = np.matmul(n,np.einsum('ij,ikj->ijk',0.5*np.log(w),n))
-      
-    return eps.reshape((3,3)) if np.shape(F) == (3,3) else \
-           eps
-
-
 def deviatoric_part(x):
     """
     Return deviatoric part of a tensor.
-      
+
     Parameters
     ----------
     x : numpy.array of shape (:,3,3) or (3,3)
@@ -89,13 +32,151 @@ def deviatoric_part(x):
 
     """
     return x - np.eye(3)*spherical_part(x) if np.shape(x) == (3,3) else \
-           x - np.einsum('ijk,i->ijk',np.broadcast_to(np.eye(3),[x.shape[0],3,3]),spherical_part(x)) 
+           x - np.einsum('ijk,i->ijk',np.broadcast_to(np.eye(3),[x.shape[0],3,3]),spherical_part(x))
+
+
+def eigenvalues(x):
+    """
+    Return the eigenvalues, i.e. principal components, of a symmetric tensor.
+
+    The eigenvalues are sorted in ascending order, each repeated according to
+    its multiplicity.
+
+    Parameters
+    ----------
+    x : numpy.array of shape (:,3,3) or (3,3)
+      Symmetric tensor of which the eigenvalues are computed.
+
+    """
+    return np.linalg.eigvalsh(symmetric(x))
+
+
+def eigenvectors(x,RHS=False):
+    """
+    Return eigenvectors of a symmetric tensor.
+
+    The eigenvalues are sorted in ascending order of their associated eigenvalues.
+
+    Parameters
+    ----------
+    x : numpy.array of shape (:,3,3) or (3,3)
+      Symmetric tensor of which the eigenvectors are computed.
+    RHS: bool, optional
+      Enforce right-handed coordinate system. Default is False.
+
+    """
+    (u,v) = np.linalg.eigh(symmetric(x))
+
+    if RHS:
+        if np.shape(x) == (3,3):
+            if np.linalg.det(v) < 0.0: v[:,2] *= -1.0
+        else:
+            v[np.linalg.det(v) < 0.0,:,2] *= -1.0
+    return v
+
+
+def left_stretch(x):
+    """
+    Return the left stretch of a tensor.
+
+    Parameters
+    ----------
+    x : numpy.array of shape (:,3,3) or (3,3)
+      Tensor of which the left stretch is computed.
+
+    """
+    return __polar_decomposition(x,'V')[0]
+
+
+def maximum_shear(x):
+    """
+    Return the maximum shear component of a symmetric tensor.
+
+    Parameters
+    ----------
+    x : numpy.array of shape (:,3,3) or (3,3)
+      Symmetric tensor of which the maximum shear is computed.
+
+    """
+    w = eigenvalues(x)
+    return (w[0] - w[2])*0.5 if np.shape(x) == (3,3) else \
+           (w[:,0] - w[:,2])*0.5
+
+
+def Mises_strain(epsilon):
+    """
+    Return the Mises equivalent of a strain tensor.
+
+    Parameters
+    ----------
+    epsilon : numpy.array of shape (:,3,3) or (3,3)
+      Symmetric strain tensor of which the von Mises equivalent is computed.
+
+    """
+    return __Mises(epsilon,2.0/3.0)
+
+
+def Mises_stress(sigma):
+    """
+    Return the Mises equivalent of a stress tensor.
+
+    Parameters
+    ----------
+    sigma : numpy.array of shape (:,3,3) or (3,3)
+      Symmetric stress tensor of which the von Mises equivalent is computed.
+
+    """
+    return __Mises(sigma,3.0/2.0)
+
+
+def PK2(P,F):
+    """
+    Calculate second Piola-Kirchhoff stress from first Piola-Kirchhoff stress and deformation gradient.
+
+    Parameters
+    ----------
+    P : numpy.array of shape (:,3,3) or (3,3)
+      1. Piola-Kirchhoff stress.
+    F : numpy.array of shape (:,3,3) or (3,3)
+      Deformation gradient.
+
+    """
+    if np.shape(F) == np.shape(P) == (3,3):
+        S = np.dot(np.linalg.inv(F),P)
+    else:
+        S = np.einsum('ijk,ikl->ijl',np.linalg.inv(F),P)
+    return symmetric(S)
+
+def right_stretch(x):
+    """
+    Return the right stretch of a tensor.
+
+    Parameters
+    ----------
+    x : numpy.array of shape (:,3,3) or (3,3)
+      Tensor of which the right stretch is computed.
+
+    """
+    return __polar_decomposition(x,'U')[0]
+
+
+def rotational_part(x):
+    """
+    Return the rotational part of a tensor.
+
+    Parameters
+    ----------
+    x : numpy.array of shape (:,3,3) or (3,3)
+      Tensor of which the rotational part is computed.
+
+    """
+    return __polar_decomposition(x,'R')[0]
 
 
 def spherical_part(x,tensor=False):
     """
     Return spherical (hydrostatic) part of a tensor.
-    
+
     Parameters
     ----------
     x : numpy.array of shape (:,3,3) or (3,3)
@@ -113,42 +194,50 @@ def spherical_part(x,tensor=False):
             return sph
         else:
             return np.einsum('ijk,i->ijk',np.broadcast_to(np.eye(3),(x.shape[0],3,3)),sph)
-    
+
-    
+
-def Mises_stress(sigma):
+def strain_tensor(F,t,m):
     """
-    Return the Mises equivalent of a stress tensor.
+    Return strain tensor calculated from deformation gradient.
-    
+
+    For details refer to https://en.wikipedia.org/wiki/Finite_strain_theory and
+    https://de.wikipedia.org/wiki/Verzerrungstensor
+
     Parameters
     ----------
-    sigma : numpy.array of shape (:,3,3) or (3,3)
+    F : numpy.array of shape (:,3,3) or (3,3)
-      Symmetric stress tensor of which the von Mises equivalent is computed.
+      Deformation gradient.
+    t : {‘V’, ‘U’}
+      Type of the polar decomposition, ‘V’ for left stretch tensor and ‘U’ for right stretch tensor.
+    m : float
+      Order of the strain.
 
     """
-    s = deviatoric_part(sigma)
+    F_ = F.reshape((1,3,3)) if F.shape == (3,3) else F
-    return np.sqrt(3.0/2.0*(np.sum(s**2.0))) if np.shape(sigma) == (3,3) else \
+    if   t == 'V':
-           np.sqrt(3.0/2.0*np.einsum('ijk->i',s**2.0))
+        B   = np.matmul(F_,transpose(F_))
-    
+        w,n = np.linalg.eigh(B)
-    
+    elif t == 'U':
-def Mises_strain(epsilon):
+        C   = np.matmul(transpose(F_),F_)
-    """
+        w,n = np.linalg.eigh(C)
-    Return the Mises equivalent of a strain tensor.
-    
-    Parameters
-    ----------
-    epsilon : numpy.array of shape (:,3,3) or (3,3)
-      Symmetric strain tensor of which the von Mises equivalent is computed.
 
-    """
+    if   m > 0.0:
-    s = deviatoric_part(epsilon)
+        eps = 1.0/(2.0*abs(m)) * (+ np.matmul(n,np.einsum('ij,ikj->ijk',w**m,n))
-    return np.sqrt(2.0/3.0*(np.sum(s**2.0))) if np.shape(epsilon) == (3,3) else \
+                                  - np.broadcast_to(np.eye(3),[F_.shape[0],3,3]))
-           np.sqrt(2.0/3.0*np.einsum('ijk->i',s**2.0))
+    elif m < 0.0:
+        eps = 1.0/(2.0*abs(m)) * (- np.matmul(n,np.einsum('ij,ikj->ijk',w**m,n))
+                                  + np.broadcast_to(np.eye(3),[F_.shape[0],3,3]))
+    else:
+        eps = np.matmul(n,np.einsum('ij,ikj->ijk',0.5*np.log(w),n))
+
+    return eps.reshape((3,3)) if np.shape(F) == (3,3) else \
+           eps
 
 
 def symmetric(x):
     """
     Return the symmetrized tensor.
-    
+
     Parameters
     ----------
     x : numpy.array of shape (:,3,3) or (3,3)
@@ -158,43 +247,10 @@ def symmetric(x):
     return (x+transpose(x))*0.5
 
 
-def maximum_shear(x):
-    """
-    Return the maximum shear component of a symmetric tensor.
-    
-    Parameters
-    ----------
-    x : numpy.array of shape (:,3,3) or (3,3)
-      Symmetric tensor of which the maximum shear is computed.
-
-    """
-    w = np.linalg.eigvalsh(symmetric(x))                                                            # eigenvalues in ascending order
-    return (w[2] - w[0])*0.5 if np.shape(x) == (3,3) else \
-           (w[:,2] - w[:,0])*0.5
-    
-    
-def principal_components(x):
-    """
-    Return the principal components of a symmetric tensor.
-    
-    The principal components (eigenvalues) are sorted in descending order, each repeated according to
-    its multiplicity.
-    
-    Parameters
-    ----------
-    x : numpy.array of shape (:,3,3) or (3,3)
-      Symmetric tensor of which the principal compontents are computed.
-
-    """
-    w = np.linalg.eigvalsh(symmetric(x))                                                            # eigenvalues in ascending order
-    return w[::-1] if np.shape(x) == (3,3) else \
-           w[:,::-1]
-    
-    
 def transpose(x):
     """
     Return the transpose of a tensor.
-      
+
     Parameters
     ----------
     x : numpy.array of shape (:,3,3) or (3,3)
@@ -205,62 +261,23 @@ def transpose(x):
            np.transpose(x,(0,2,1))
 
 
-def rotational_part(x):
-    """
-    Return the rotational part of a tensor.
-      
-    Parameters
-    ----------
-    x : numpy.array of shape (:,3,3) or (3,3)
-      Tensor of which the rotational part is computed.
-
-    """
-    return __polar_decomposition(x,'R')[0]
-
-
-def left_stretch(x):
-    """
-    Return the left stretch of a tensor.
-      
-    Parameters
-    ----------
-    x : numpy.array of shape (:,3,3) or (3,3)
-      Tensor of which the left stretch is computed.
-
-    """
-    return __polar_decomposition(x,'V')[0]
-  
-  
-def right_stretch(x):
-    """
-    Return the right stretch of a tensor.
-      
-    Parameters
-    ----------
-    x : numpy.array of shape (:,3,3) or (3,3)
-      Tensor of which the right stretch is computed.
-
-    """
-    return __polar_decomposition(x,'U')[0]
-
-
 def __polar_decomposition(x,requested):
     """
     Singular value decomposition.
-      
+
     Parameters
     ----------
     x : numpy.array of shape (:,3,3) or (3,3)
       Tensor of which the singular values are computed.
     requested : iterable of str
-      Requested outputs: ‘R’ for the rotation tensor, 
+      Requested outputs: ‘R’ for the rotation tensor,
       ‘V’ for left stretch tensor and ‘U’ for right stretch tensor.
 
     """
     u, s, vh = np.linalg.svd(x)
     R = np.dot(u,vh) if np.shape(x) == (3,3) else \
         np.einsum('ijk,ikl->ijl',u,vh)
-    
+
     output = []
     if 'R' in requested:
         output.append(R)
@@ -268,5 +285,22 @@ def __polar_decomposition(x,requested):
         output.append(np.dot(x,R.T) if np.shape(x) == (3,3) else np.einsum('ijk,ilk->ijl',x,R))
     if 'U' in requested:
         output.append(np.dot(R.T,x) if np.shape(x) == (3,3) else np.einsum('ikj,ikl->ijl',R,x))
-    
+
     return tuple(output)
+
+
+def __Mises(x,s):
+    """
+    Base equation for Mises equivalent of a stres or strain tensor.
+
+    Parameters
+    ----------
+    x : numpy.array of shape (:,3,3) or (3,3)
+      Symmetric tensor of which the von Mises equivalent is computed.
+    s : float
+      Scaling factor (2/3 for strain, 3/2 for stress).
+    
+    """
+    d = deviatoric_part(x)
+    return np.sqrt(s*(np.sum(d**2.0))) if np.shape(x) == (3,3) else \
+           np.sqrt(s*np.einsum('ijk->i',d**2.0))

python/damask/util.py

@@ -3,14 +3,18 @@ import time
 import os
 import subprocess
 import shlex
+from fractions import Fraction
+from functools import reduce
 from optparse import Option
 from queue import Queue
 from threading import Thread
 
+import numpy as np
+
 class bcolors:
     """
     ASCII Colors (Blender code).
-    
+
     https://svn.blender.org/svnroot/bf-blender/trunk/blender/build_files/scons/tools/bcolors.py
     http://stackoverflow.com/questions/287871/print-in-terminal-with-colors-using-python
     """
@@ -36,7 +40,7 @@ class bcolors:
         self.BOLD = ''
         self.UNDERLINE = ''
         self.CROSSOUT = ''
-    
+
 
 # -----------------------------
 def srepr(arg,glue = '\n'):
@@ -159,11 +163,29 @@ def progressBar(iteration, total, prefix='', bar_length=50):
 
   if iteration == total: sys.stderr.write('\n')
   sys.stderr.flush()
- 
+
+
+def scale_to_coprime(v):
+  """Scale vector to co-prime (relatively prime) integers."""
+  MAX_DENOMINATOR = 1000
+
+  def get_square_denominator(x):
+    """Denominator of the square of a number."""
+    return Fraction(x ** 2).limit_denominator(MAX_DENOMINATOR).denominator
+
+  def lcm(a, b):
+    """Least common multiple."""
+    return a * b // np.gcd(a, b)
+
+  denominators = [int(get_square_denominator(i)) for i in v]
+  s = reduce(lcm, denominators) ** 0.5
+  m = (np.array(v)*s).astype(np.int)
+  return m//reduce(np.gcd,m)
+
 
 class return_message():
   """Object with formatted return message."""
-  
+
   def __init__(self,message):
     """
     Sets return message.
@@ -175,7 +197,7 @@ class return_message():
 
     """
     self.message = message
-  
+
   def __repr__(self):
     """Return message suitable for interactive shells."""
     return srepr(self.message)

python/tests/test_DADF5.py

@@ -40,7 +40,7 @@ class TestDADF5:
         assert np.allclose(in_memory,in_file)
 
     def test_add_calculation(self,default):
-        default.add_calculation('2.0*np.abs(#F#)-1.0','x','-','test')
+        default.add_calculation('x','2.0*np.abs(#F#)-1.0','-','my notes')
         loc = {'F':    default.get_dataset_location('F'),
                'x':    default.get_dataset_location('x')}
         in_memory = 2.0*np.abs(default.read_dataset(loc['F'],0))-1.0
@@ -52,8 +52,8 @@ class TestDADF5:
         loc = {'F':    default.get_dataset_location('F'),
                'P':    default.get_dataset_location('P'),
                'sigma':default.get_dataset_location('sigma')}
-        in_memory = mechanics.Cauchy(default.read_dataset(loc['F'],0),
+        in_memory = mechanics.Cauchy(default.read_dataset(loc['P'],0),
-                                     default.read_dataset(loc['P'],0))
+                                     default.read_dataset(loc['F'],0))
         in_file   = default.read_dataset(loc['sigma'],0)
         assert np.allclose(in_memory,in_file)
 
@@ -73,6 +73,54 @@ class TestDADF5:
         in_file   = default.read_dataset(loc['s_P'],0)
         assert np.allclose(in_memory,in_file)
 
+    def test_add_eigenvalues(self,default):
+        default.add_Cauchy('P','F')
+        default.add_eigenvalues('sigma')
+        loc = {'sigma'        :default.get_dataset_location('sigma'),
+               'lambda(sigma)':default.get_dataset_location('lambda(sigma)')}
+        in_memory = mechanics.eigenvalues(default.read_dataset(loc['sigma'],0))
+        in_file   = default.read_dataset(loc['lambda(sigma)'],0)
+        assert np.allclose(in_memory,in_file)
+
+    def test_add_eigenvectors(self,default):
+        default.add_Cauchy('P','F')
+        default.add_eigenvectors('sigma')
+        loc = {'sigma'   :default.get_dataset_location('sigma'),
+               'v(sigma)':default.get_dataset_location('v(sigma)')}
+        in_memory = mechanics.eigenvectors(default.read_dataset(loc['sigma'],0))
+        in_file   = default.read_dataset(loc['v(sigma)'],0)
+        assert np.allclose(in_memory,in_file)
+
+    def test_add_maximum_shear(self,default):
+        default.add_Cauchy('P','F')
+        default.add_maximum_shear('sigma')
+        loc = {'sigma'           :default.get_dataset_location('sigma'),
+               'max_shear(sigma)':default.get_dataset_location('max_shear(sigma)')}
+        in_memory = mechanics.maximum_shear(default.read_dataset(loc['sigma'],0)).reshape(-1,1)
+        in_file   = default.read_dataset(loc['max_shear(sigma)'],0)
+        assert np.allclose(in_memory,in_file)
+
+    def test_add_Mises_strain(self,default):
+        t = ['V','U'][np.random.randint(0,2)]
+        m = np.random.random()*2.0 - 1.0
+        default.add_strain_tensor('F',t,m)
+        label = 'epsilon_{}^{}(F)'.format(t,m)
+        default.add_Mises(label)
+        loc = {label      :default.get_dataset_location(label),
+               label+'_vM':default.get_dataset_location(label+'_vM')}
+        in_memory = mechanics.Mises_strain(default.read_dataset(loc[label],0)).reshape(-1,1)
+        in_file   = default.read_dataset(loc[label+'_vM'],0)
+        assert np.allclose(in_memory,in_file)
+
+    def test_add_Mises_stress(self,default):
+        default.add_Cauchy('P','F')
+        default.add_Mises('sigma')
+        loc = {'sigma'   :default.get_dataset_location('sigma'),
+               'sigma_vM':default.get_dataset_location('sigma_vM')}
+        in_memory = mechanics.Mises_stress(default.read_dataset(loc['sigma'],0)).reshape(-1,1)
+        in_file   = default.read_dataset(loc['sigma_vM'],0)
+        assert np.allclose(in_memory,in_file)
+
     def test_add_norm(self,default):
         default.add_norm('F',1)
         loc = {'F':    default.get_dataset_location('F'),
@@ -81,6 +129,24 @@ class TestDADF5:
         in_file   = default.read_dataset(loc['|F|_1'],0)
         assert np.allclose(in_memory,in_file)
 
+    def test_add_PK2(self,default):
+        default.add_PK2('P','F')
+        loc = {'F':default.get_dataset_location('F'),
+               'P':default.get_dataset_location('P'),
+               'S':default.get_dataset_location('S')}
+        in_memory = mechanics.PK2(default.read_dataset(loc['P'],0),
+                                  default.read_dataset(loc['F'],0))
+        in_file   = default.read_dataset(loc['S'],0)
+        assert np.allclose(in_memory,in_file)
+
+    def test_add_rotational_part(self,default):
+        default.add_rotational_part('F')
+        loc = {'F':    default.get_dataset_location('F'),
+               'R(F)': default.get_dataset_location('R(F)')}
+        in_memory = mechanics.rotational_part(default.read_dataset(loc['F'],0))
+        in_file   = default.read_dataset(loc['R(F)'],0)
+        assert np.allclose(in_memory,in_file)
+
     def test_add_spherical(self,default):
         default.add_spherical('P')
         loc = {'P':   default.get_dataset_location('P'),
@@ -88,3 +154,30 @@ class TestDADF5:
         in_memory = mechanics.spherical_part(default.read_dataset(loc['P'],0)).reshape(-1,1)
         in_file   = default.read_dataset(loc['p_P'],0)
         assert np.allclose(in_memory,in_file)
+
+    def test_add_strain(self,default):
+        t = ['V','U'][np.random.randint(0,2)]
+        m = np.random.random()*2.0 - 1.0
+        default.add_strain_tensor('F',t,m)
+        label = 'epsilon_{}^{}(F)'.format(t,m)
+        loc = {'F':   default.get_dataset_location('F'),
+               label: default.get_dataset_location(label)}
+        in_memory = mechanics.strain_tensor(default.read_dataset(loc['F'],0),t,m)
+        in_file   = default.read_dataset(loc[label],0)
+        assert np.allclose(in_memory,in_file)
+
+    def test_add_stretch_right(self,default):
+        default.add_stretch_tensor('F','U')
+        loc = {'F':    default.get_dataset_location('F'),
+               'U(F)': default.get_dataset_location('U(F)')}
+        in_memory = mechanics.right_stretch(default.read_dataset(loc['F'],0))
+        in_file   = default.read_dataset(loc['U(F)'],0)
+        assert np.allclose(in_memory,in_file)
+
+    def test_add_stretch_left(self,default):
+        default.add_stretch_tensor('F','V')
+        loc = {'F':    default.get_dataset_location('F'),
+               'V(F)': default.get_dataset_location('V(F)')}
+        in_memory = mechanics.left_stretch(default.read_dataset(loc['F'],0))
+        in_file   = default.read_dataset(loc['V(F)'],0)
+        assert np.allclose(in_memory,in_file)

python/tests/test_mechanics.py

@@ -2,187 +2,224 @@ import numpy as np
 from damask import mechanics
 
 class TestMechanics:
-   
+
     n = 1000
     c = np.random.randint(n)
-   
+
 
     def test_vectorize_Cauchy(self):
-         P = np.random.random((self.n,3,3))
+        P = np.random.random((self.n,3,3))
-         F = np.random.random((self.n,3,3))
+        F = np.random.random((self.n,3,3))
-         assert np.allclose(mechanics.Cauchy(F,P)[self.c],
+        assert np.allclose(mechanics.Cauchy(P,F)[self.c],
-                            mechanics.Cauchy(F[self.c],P[self.c]))
+                           mechanics.Cauchy(P[self.c],F[self.c]))
-
-
-    def test_vectorize_strain_tensor(self):
-         F = np.random.random((self.n,3,3))
-         t = ['V','U'][np.random.randint(0,2)]
-         m = np.random.random()*10. -5.0
-         assert np.allclose(mechanics.strain_tensor(F,t,m)[self.c],
-                            mechanics.strain_tensor(F[self.c],t,m))
-
 
     def test_vectorize_deviatoric_part(self):
-         x = np.random.random((self.n,3,3))
+        x = np.random.random((self.n,3,3))
-         assert np.allclose(mechanics.deviatoric_part(x)[self.c],
+        assert np.allclose(mechanics.deviatoric_part(x)[self.c],
-                            mechanics.deviatoric_part(x[self.c]))
+                           mechanics.deviatoric_part(x[self.c]))
 
+    def test_vectorize_eigenvalues(self):
+        x = np.random.random((self.n,3,3))
+        assert np.allclose(mechanics.eigenvalues(x)[self.c],
+                           mechanics.eigenvalues(x[self.c]))
 
-    def test_vectorize_spherical_part(self):
+    def test_vectorize_eigenvectors(self):
-         x = np.random.random((self.n,3,3))
+        x = np.random.random((self.n,3,3))
-         assert np.allclose(mechanics.spherical_part(x,True)[self.c],
+        assert np.allclose(mechanics.eigenvectors(x)[self.c],
-                            mechanics.spherical_part(x[self.c],True))
+                           mechanics.eigenvectors(x[self.c]))
-
-
-    def test_vectorize_Mises_stress(self):
-         sigma = np.random.random((self.n,3,3))
-         assert np.allclose(mechanics.Mises_stress(sigma)[self.c],
-                            mechanics.Mises_stress(sigma[self.c]))
-
-
-    def test_vectorize_Mises_strain(self):
-         epsilon = np.random.random((self.n,3,3))
-         assert np.allclose(mechanics.Mises_strain(epsilon)[self.c],
-                            mechanics.Mises_strain(epsilon[self.c]))
-
-
-    def test_vectorize_symmetric(self):
-         x = np.random.random((self.n,3,3))
-         assert np.allclose(mechanics.symmetric(x)[self.c],
-                            mechanics.symmetric(x[self.c]))
-
-
-    def test_vectorize_maximum_shear(self):
-         x = np.random.random((self.n,3,3))
-         assert np.allclose(mechanics.maximum_shear(x)[self.c],
-                            mechanics.maximum_shear(x[self.c]))
-
-
-    def test_vectorize_principal_components(self):
-         x = np.random.random((self.n,3,3))
-         assert np.allclose(mechanics.principal_components(x)[self.c],
-                            mechanics.principal_components(x[self.c]))
-
-
-    def test_vectorize_transpose(self):
-         x = np.random.random((self.n,3,3))
-         assert np.allclose(mechanics.transpose(x)[self.c],
-                            mechanics.transpose(x[self.c]))
-
-
-    def test_vectorize_rotational_part(self):
-         x = np.random.random((self.n,3,3))
-         assert np.allclose(mechanics.rotational_part(x)[self.c],
-                            mechanics.rotational_part(x[self.c]))
-
 
     def test_vectorize_left_stretch(self):
-         x = np.random.random((self.n,3,3))
+        x = np.random.random((self.n,3,3))
-         assert np.allclose(mechanics.left_stretch(x)[self.c],
+        assert np.allclose(mechanics.left_stretch(x)[self.c],
-                            mechanics.left_stretch(x[self.c]))
+                           mechanics.left_stretch(x[self.c]))
 
+    def test_vectorize_maximum_shear(self):
+        x = np.random.random((self.n,3,3))
+        assert np.allclose(mechanics.maximum_shear(x)[self.c],
+                           mechanics.maximum_shear(x[self.c]))
+
+    def test_vectorize_Mises_strain(self):
+        epsilon = np.random.random((self.n,3,3))
+        assert np.allclose(mechanics.Mises_strain(epsilon)[self.c],
+                           mechanics.Mises_strain(epsilon[self.c]))
+
+    def test_vectorize_Mises_stress(self):
+        sigma = np.random.random((self.n,3,3))
+        assert np.allclose(mechanics.Mises_stress(sigma)[self.c],
+                           mechanics.Mises_stress(sigma[self.c]))
+
+    def test_vectorize_PK2(self):
+        F = np.random.random((self.n,3,3))
+        P = np.random.random((self.n,3,3))
+        assert np.allclose(mechanics.PK2(P,F)[self.c],
+                           mechanics.PK2(P[self.c],F[self.c]))
 
     def test_vectorize_right_stretch(self):
-         x = np.random.random((self.n,3,3))
+        x = np.random.random((self.n,3,3))
-         assert np.allclose(mechanics.right_stretch(x)[self.c],
+        assert np.allclose(mechanics.right_stretch(x)[self.c],
-                            mechanics.right_stretch(x[self.c]))
+                           mechanics.right_stretch(x[self.c]))
+
+    def test_vectorize_rotational_part(self):
+        x = np.random.random((self.n,3,3))
+        assert np.allclose(mechanics.rotational_part(x)[self.c],
+                           mechanics.rotational_part(x[self.c]))
+
+    def test_vectorize_spherical_part(self):
+        x = np.random.random((self.n,3,3))
+        assert np.allclose(mechanics.spherical_part(x,True)[self.c],
+                           mechanics.spherical_part(x[self.c],True))
+
+    def test_vectorize_strain_tensor(self):
+        F = np.random.random((self.n,3,3))
+        t = ['V','U'][np.random.randint(0,2)]
+        m = np.random.random()*10. -5.0
+        assert np.allclose(mechanics.strain_tensor(F,t,m)[self.c],
+                           mechanics.strain_tensor(F[self.c],t,m))
+
+    def test_vectorize_symmetric(self):
+        x = np.random.random((self.n,3,3))
+        assert np.allclose(mechanics.symmetric(x)[self.c],
+                           mechanics.symmetric(x[self.c]))
+
+    def test_vectorize_transpose(self):
+        x = np.random.random((self.n,3,3))
+        assert np.allclose(mechanics.transpose(x)[self.c],
+                           mechanics.transpose(x[self.c]))
 
 
     def test_Cauchy(self):
-         """Ensure Cauchy stress is symmetrized 1. Piola-Kirchhoff stress for no deformation."""
+        """Ensure Cauchy stress is symmetrized 1. Piola-Kirchhoff stress for no deformation."""
-         P = np.random.random((self.n,3,3))
+        P = np.random.random((self.n,3,3))
-         assert np.allclose(mechanics.Cauchy(np.broadcast_to(np.eye(3),(self.n,3,3)),P),
+        assert np.allclose(mechanics.Cauchy(P,np.broadcast_to(np.eye(3),(self.n,3,3))),
-                            mechanics.symmetric(P))
+                           mechanics.symmetric(P))
+
 
     def test_polar_decomposition(self):
-         """F = RU = VR."""
+        """F = RU = VR."""
-         F = np.broadcast_to(np.eye(3),[self.n,3,3])*np.random.random((self.n,3,3))
+        F = np.broadcast_to(np.eye(3),[self.n,3,3])*np.random.random((self.n,3,3))
-         R = mechanics.rotational_part(F)
+        R = mechanics.rotational_part(F)
-         V = mechanics.left_stretch(F)
+        V = mechanics.left_stretch(F)
-         U = mechanics.right_stretch(F)
+        U = mechanics.right_stretch(F)
-         assert np.allclose(np.matmul(R,U),
+        assert np.allclose(np.matmul(R,U),
-                            np.matmul(V,R))
+                           np.matmul(V,R))
-   
+
+
+    def test_PK2(self):
+        """Ensure 2. Piola-Kirchhoff stress is symmetrized 1. Piola-Kirchhoff stress for no deformation."""
+        P = np.random.random((self.n,3,3))
+        assert np.allclose(mechanics.PK2(P,np.broadcast_to(np.eye(3),(self.n,3,3))),
+                           mechanics.symmetric(P))
+
 
     def test_strain_tensor_no_rotation(self):
-         """Ensure that left and right stretch give same results for no rotation."""
+        """Ensure that left and right stretch give same results for no rotation."""
-         F = np.broadcast_to(np.eye(3),[self.n,3,3])*np.random.random((self.n,3,3))
+        F = np.broadcast_to(np.eye(3),[self.n,3,3])*np.random.random((self.n,3,3))
-         m = np.random.random()*20.0-10.0
+        m = np.random.random()*20.0-10.0
-         assert np.allclose(mechanics.strain_tensor(F,'U',m),
+        assert np.allclose(mechanics.strain_tensor(F,'U',m),
-                            mechanics.strain_tensor(F,'V',m))
+                           mechanics.strain_tensor(F,'V',m))
-   
+
     def test_strain_tensor_rotation_equivalence(self):
-         """Ensure that left and right strain differ only by a rotation."""
+        """Ensure that left and right strain differ only by a rotation."""
-         F = np.broadcast_to(np.eye(3),[self.n,3,3]) + (np.random.random((self.n,3,3))*0.5 - 0.25)
+        F = np.broadcast_to(np.eye(3),[self.n,3,3]) + (np.random.random((self.n,3,3))*0.5 - 0.25)
-         m = np.random.random()*5.0-2.5
+        m = np.random.random()*5.0-2.5
-         assert np.allclose(np.linalg.det(mechanics.strain_tensor(F,'U',m)),
+        assert np.allclose(np.linalg.det(mechanics.strain_tensor(F,'U',m)),
-                            np.linalg.det(mechanics.strain_tensor(F,'V',m)))
+                           np.linalg.det(mechanics.strain_tensor(F,'V',m)))
 
     def test_strain_tensor_rotation(self):
-         """Ensure that pure rotation results in no strain."""
+        """Ensure that pure rotation results in no strain."""
-         F = mechanics.rotational_part(np.random.random((self.n,3,3)))
+        F = mechanics.rotational_part(np.random.random((self.n,3,3)))
-         t = ['V','U'][np.random.randint(0,2)]
+        t = ['V','U'][np.random.randint(0,2)]
-         m = np.random.random()*2.0 - 1.0
+        m = np.random.random()*2.0 - 1.0
-         assert np.allclose(mechanics.strain_tensor(F,t,m),
+        assert np.allclose(mechanics.strain_tensor(F,t,m),
-                            0.0)
+                           0.0)
-   
-    def test_rotation_determinant(self):
-         """
-         Ensure that the determinant of the rotational part is +- 1.
 
-         Should be +1, but random F might contain a reflection.
+    def test_rotation_determinant(self):
-         """
+        """
-         x = np.random.random((self.n,3,3))
+        Ensure that the determinant of the rotational part is +- 1.
-         assert np.allclose(np.abs(np.linalg.det(mechanics.rotational_part(x))),
+
-                            1.0)
+        Should be +1, but random F might contain a reflection.
+        """
+        x = np.random.random((self.n,3,3))
+        assert np.allclose(np.abs(np.linalg.det(mechanics.rotational_part(x))),
+                           1.0)
 
 
     def test_spherical_deviatoric_part(self):
-         """Ensure that full tensor is sum of spherical and deviatoric part."""
+        """Ensure that full tensor is sum of spherical and deviatoric part."""
-         x = np.random.random((self.n,3,3))
+        x = np.random.random((self.n,3,3))
-         sph = mechanics.spherical_part(x,True)
+        sph = mechanics.spherical_part(x,True)
-         assert np.allclose(sph + mechanics.deviatoric_part(x),
+        assert np.allclose(sph + mechanics.deviatoric_part(x),
-                            x)
+                           x)
 
     def test_deviatoric_Mises(self):
-         """Ensure that Mises equivalent stress depends only on deviatoric part."""
+        """Ensure that Mises equivalent stress depends only on deviatoric part."""
-         x = np.random.random((self.n,3,3))
+        x = np.random.random((self.n,3,3))
-         full = mechanics.Mises_stress(x)
+        full = mechanics.Mises_stress(x)
-         dev  = mechanics.Mises_stress(mechanics.deviatoric_part(x))
+        dev  = mechanics.Mises_stress(mechanics.deviatoric_part(x))
-         assert np.allclose(full,
+        assert np.allclose(full,
-                            dev)
+                           dev)
 
     def test_spherical_mapping(self):
-         """Ensure that mapping to tensor is correct."""
+        """Ensure that mapping to tensor is correct."""
-         x = np.random.random((self.n,3,3))
+        x = np.random.random((self.n,3,3))
-         tensor = mechanics.spherical_part(x,True)
+        tensor = mechanics.spherical_part(x,True)
-         scalar = mechanics.spherical_part(x)
+        scalar = mechanics.spherical_part(x)
-         assert np.allclose(np.linalg.det(tensor),
+        assert np.allclose(np.linalg.det(tensor),
-                            scalar**3.0)
+                           scalar**3.0)
 
     def test_spherical_Mises(self):
-         """Ensure that Mises equivalent strrain of spherical strain is 0."""
+        """Ensure that Mises equivalent strrain of spherical strain is 0."""
-         x = np.random.random((self.n,3,3))
+        x = np.random.random((self.n,3,3))
-         sph = mechanics.spherical_part(x,True)
+        sph = mechanics.spherical_part(x,True)
-         assert np.allclose(mechanics.Mises_strain(sph),
+        assert np.allclose(mechanics.Mises_strain(sph),
-                            0.0)
+                           0.0)
 
     def test_symmetric(self):
-         """Ensure that a symmetric tensor is half of the sum of a tensor and its transpose."""
+        """Ensure that a symmetric tensor is half of the sum of a tensor and its transpose."""
-         x = np.random.random((self.n,3,3))
+        x = np.random.random((self.n,3,3))
-         assert np.allclose(mechanics.symmetric(x)*2.0,
+        assert np.allclose(mechanics.symmetric(x)*2.0,
-                            mechanics.transpose(x)+x)
+                           mechanics.transpose(x)+x)
 
 
     def test_transpose(self):
-         """Ensure that a symmetric tensor equals its transpose."""
+        """Ensure that a symmetric tensor equals its transpose."""
-         x = mechanics.symmetric(np.random.random((self.n,3,3)))
+        x = mechanics.symmetric(np.random.random((self.n,3,3)))
-         assert np.allclose(mechanics.transpose(x),
+        assert np.allclose(mechanics.transpose(x),
-                            x)
+                           x)
 
 
     def test_Mises(self):
-         """Ensure that equivalent stress is 3/2 of equivalent strain."""
+        """Ensure that equivalent stress is 3/2 of equivalent strain."""
-         x = np.random.random((self.n,3,3))
+        x = np.random.random((self.n,3,3))
-         assert np.allclose(mechanics.Mises_stress(x)/mechanics.Mises_strain(x),
+        assert np.allclose(mechanics.Mises_stress(x)/mechanics.Mises_strain(x),
-                            1.5)
+                           1.5)
+
+
+    def test_eigenvalues(self):
+        """Ensure that the characteristic polynomial can be solved."""
+        A = mechanics.symmetric(np.random.random((self.n,3,3)))
+        lambd = mechanics.eigenvalues(A)
+        s = np.random.randint(self.n)
+        for i in range(3):
+           assert np.allclose(np.linalg.det(A[s]-lambd[s,i]*np.eye(3)),.0)
+
+    def test_eigenvalues_and_vectors(self):
+        """Ensure that eigenvalues and -vectors are the solution to the characteristic polynomial."""
+        A = mechanics.symmetric(np.random.random((self.n,3,3)))
+        lambd = mechanics.eigenvalues(A)
+        x     = mechanics.eigenvectors(A)
+        s = np.random.randint(self.n)
+        for i in range(3):
+           assert np.allclose(np.dot(A[s]-lambd[s,i]*np.eye(3),x[s,:,i]),.0)
+
+    def test_eigenvectors_RHS(self):
+        """Ensure that RHS coordinate system does only change sign of determinant."""
+        A = mechanics.symmetric(np.random.random((self.n,3,3)))
+        LRHS = np.linalg.det(mechanics.eigenvectors(A,RHS=False))
+        RHS  = np.linalg.det(mechanics.eigenvectors(A,RHS=True))
+        assert np.allclose(np.abs(LRHS),RHS)
+
+    def test_spherical_no_shear(self):
+        """Ensure that sherical stress has max shear of 0.0."""
+        A = mechanics.spherical_part(mechanics.symmetric(np.random.random((self.n,3,3))),True)
+        assert np.allclose(mechanics.maximum_shear(A),0.0)