diff --git a/processing/post/addPK2.py b/processing/post/addPK2.py
index f38753619..45afdbcd1 100755
--- a/processing/post/addPK2.py
+++ b/processing/post/addPK2.py
@@ -2,10 +2,9 @@
 
 import os
 import sys
+from io import StringIO
 from optparse import OptionParser
 
-import numpy as np
-
 import damask
 
 
@@ -36,53 +35,16 @@ parser.set_defaults(defgrad = 'f',
                    )
 
 (options,filenames) = parser.parse_args()
-
-# --- loop over input files -------------------------------------------------------------------------
-
 if filenames == []: filenames = [None]
 
 for name in filenames:
-  try:
-    table = damask.ASCIItable(name = name,
-                              buffered = False)
-  except: continue
-  damask.util.report(scriptName,name)
+    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()
+    table.add_array('S',
+                    damask.mechanics.PK2(table.get_array(options.defgrad).reshape(-1,3,3),
+                                         table.get_array(options.stress).reshape(-1,3,3)).reshape(-1,9),
+                    scriptID+' '+' '.join(sys.argv[1:]))
 
-# ------------------------------------------ sanity checks ----------------------------------------
-
-  errors = []
-  column = {}
-  
-  for tensor in [options.defgrad,options.stress]:
-    dim = table.label_dimension(tensor)
-    if   dim <  0: errors.append('column {} not found.'.format(tensor))
-    elif dim != 9: errors.append('column {} is not a tensor.'.format(tensor))
-    else:
-      column[tensor] = table.label_index(tensor)
-
-  if errors != []:
-    damask.util.croak(errors)
-    table.close(dismiss = True)
-    continue
-
-# ------------------------------------------ assemble header --------------------------------------
-
-  table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:]))
-  table.labels_append(['{}_S'.format(i+1) for i in range(9)])                                       # extend ASCII header with new labels
-  table.head_write()
-
-# ------------------------------------------ process data ------------------------------------------
-  outputAlive = True
-  while outputAlive and table.data_read():                                                          # read next data line of ASCII table
-    F = np.array(list(map(float,table.data[column[options.defgrad]:column[options.defgrad]+9])),'d').reshape(3,3)
-    P = np.array(list(map(float,table.data[column[options.stress ]:column[options.stress ]+9])),'d').reshape(3,3)
-    table.data_append(list(np.dot(np.linalg.inv(F),P).reshape(9)))                                  # [S] =[P].[F-1]
-    outputAlive = table.data_write()                                                                # output processed line
-
-# ------------------------------------------ output finalization -----------------------------------
-
-  table.close()                                                                                     # close input ASCII table (works for stdin)
+    table.to_ASCII(sys.stdout if name is None else name)
diff --git a/processing/post/scaleData.py b/processing/post/scaleData.py
index 90eaeb0cc..ad57f4b7e 100755
--- a/processing/post/scaleData.py
+++ b/processing/post/scaleData.py
@@ -2,10 +2,9 @@
 
 import os
 import sys
+from io import StringIO
 from optparse import OptionParser
 
-import numpy as np
-
 import damask
 
 
diff --git a/python/damask/mechanics.py b/python/damask/mechanics.py
index 476682380..5503d7048 100644
--- a/python/damask/mechanics.py
+++ b/python/damask/mechanics.py
@@ -19,7 +19,26 @@ def Cauchy(F,P):
     else:
         sigma = np.einsum('i,ijk,ilk->ijl',1.0/np.linalg.det(F),P,F)
     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):
     """