From 3336cfc3da46973afb62e1559ac8c5908d5e8105 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Fri, 18 Oct 2019 20:50:03 +0200
Subject: [PATCH] better have tensor functions available

allows in-memory evaluation of results
---
 python/damask/__init__.py  |  6 +++++-
 python/damask/dadf5.py     | 10 ++++------
 python/damask/mechanics.py | 23 +++++++++++++++++++++++
 3 files changed, 32 insertions(+), 7 deletions(-)
 create mode 100644 python/damask/mechanics.py

diff --git a/python/damask/__init__.py b/python/damask/__init__.py
index ed6747a93..8f1e64937 100644
--- a/python/damask/__init__.py
+++ b/python/damask/__init__.py
@@ -6,6 +6,7 @@ with open(os.path.join(os.path.dirname(__file__),'VERSION')) as f:
 
 name = 'damask'
 
+# classes
 from .environment import Environment      # noqa
 from .asciitable  import ASCIItable       # noqa
     
@@ -14,8 +15,11 @@ from .colormaps   import Colormap, Color  # noqa
 from .orientation import Symmetry, Lattice, Rotation, Orientation # noqa
 from .dadf5       import DADF5 # noqa
 
-#from .block       import Block           # only one class
 from .geom        import Geom             # noqa
 from .solver      import Solver           # noqa
 from .test        import Test             # noqa
 from .util        import extendableOption # noqa
+
+# functions in modules
+from . import mechanics
+
diff --git a/python/damask/dadf5.py b/python/damask/dadf5.py
index 54d92be40..ba81fa605 100644
--- a/python/damask/dadf5.py
+++ b/python/damask/dadf5.py
@@ -7,6 +7,7 @@ import numpy as np
 
 from . import util
 from . import version
+from . import mechanics
 
 # ------------------------------------------------------------------
 class DADF5():
@@ -379,10 +380,8 @@ class DADF5():
     Resulting tensor is symmetrized as the Cauchy stress should be symmetric.
     """
     def Cauchy(F,P):
-      sigma = np.einsum('i,ijk,ilk->ijl',1.0/np.linalg.det(F['data']),P['data'],F['data'])
-      sigma = (sigma + np.transpose(sigma,(0,2,1)))*0.5                                             # enforce symmetry
       return  {
-               'data'  : sigma,
+               'data'  : mechanics.Cauchy(F['data']),P['data']),
                'label' : 'sigma',
                'meta' : {
                         'Unit' :        P['meta']['Unit'],
@@ -529,13 +528,12 @@ class DADF5():
   def add_deviator(self,x):
     """Adds the deviator of a tensor."""
     def deviator(x):
-      d = x['data']
       
-      if not np.all(np.array(d.shape[1:]) == np.array([3,3])):
+      if not np.all(np.array(x['data'].shape[1:]) == np.array([3,3])):
         raise ValueError
       
       return  {
-               'data' :  d - np.einsum('ijk,i->ijk',np.broadcast_to(np.eye(3),[d.shape[0],3,3]),np.trace(d,axis1=1,axis2=2)/3.0), 
+               'data' :  mechanics.deviator(x['data']),
                'label' : 'dev({})'.format(x['label']),
                'meta' : {
                         'Unit' :        x['meta']['Unit'],
diff --git a/python/damask/mechanics.py b/python/damask/mechanics.py
new file mode 100644
index 000000000..5f8dbc02d
--- /dev/null
+++ b/python/damask/mechanics.py
@@ -0,0 +1,23 @@
+import numpy as np
+
+def Cauchy(F,P):
+  if np.shape(F) == np.shape(P) == (3,3):
+    sigma = 1.0/np.linalg.det(F) * np.dot(F,P)
+    return (sigma+sigma.T)*0.5
+  else:
+    sigma =  np.einsum('i,ijk,ilk->ijl',1.0/np.linalg.det(F),P,F)
+    return (sigma + np.transpose(sigma,(0,2,1)))*0.5
+
+
+def deviator(x):
+  if np.shape(x) == (3,3):
+    return x - np.eye(3)*np.trace(x)/3.0
+  else:
+    return d - np.einsum('ijk,i->ijk',np.broadcast_to(np.eye(3),[d.shape[0],3,3]),np.trace(d,axis1=1,axis2=2)/3.0) 
+
+
+def spherical(x):
+  if np.shape(x) == (3,3):
+    return np.trace(x)/3.0
+  else:
+    return np.trace(x,axis1=1,axis2=2)/3.0