diff --git a/README.md b/README.md
index cd6e7bd8..a68a95aa 100644
--- a/README.md
+++ b/README.md
@@ -48,8 +48,11 @@ All notable changes to this project will be documented in this file. The format
 
 ## [Unreleased]
 
+### Added
+- Add `SolidBody.evaluate.kirchhoff_stress()` method. Contrary to the Cauchy stress method, this gives correct results in incompressible plane stress.
+
 ### Fixed
-- Fix `tovoigt()` helper for data with more or less than two trailing axes.
+- Fix `tovoigt()` helper for data with more or less than two trailing axes and 2D tensors.
 
 ## [3.1.0] - 2022-05-02
 
diff --git a/felupe/math/_tensor.py b/felupe/math/_tensor.py
index 2c049bba..dd60343c 100644
--- a/felupe/math/_tensor.py
+++ b/felupe/math/_tensor.py
@@ -279,8 +279,12 @@ def ddot(A, B, n=2, parallel=False):
 
 def tovoigt(A):
     "Convert (3, 3) tensor to (6, ) voigt notation."
-    B = np.zeros((6, *A.shape[2:]))
-    ij = [(0, 0), (1, 1), (2, 2), (0, 1), (1, 2), (0, 2)]
+    if A.shape[0] * A.shape[1] == 4:
+        B = np.zeros((3, *A.shape[2:]))
+        ij = [(0, 0), (1, 1), (0, 1)]
+    else:
+        B = np.zeros((6, *A.shape[2:]))
+        ij = [(0, 0), (1, 1), (2, 2), (0, 1), (1, 2), (0, 2)]
     for i6, (i, j) in enumerate(ij):
         B[i6] = A[i, j]
     return B
diff --git a/felupe/mechanics/_helpers.py b/felupe/mechanics/_helpers.py
index 6c7f4339..45630e33 100644
--- a/felupe/mechanics/_helpers.py
+++ b/felupe/mechanics/_helpers.py
@@ -37,12 +37,13 @@ def __init__(self, vector, matrix):
 class Evaluate:
     "A class with evaluate methods of a SolidBody."
 
-    def __init__(self, gradient, hessian, cauchy_stress=None):
+    def __init__(self, gradient, hessian, cauchy_stress=None, kirchhoff_stress=None):
         self.gradient = gradient
         self.hessian = hessian
 
         if cauchy_stress is not None:
             self.cauchy_stress = cauchy_stress
+            self.kirchhoff_stress = kirchhoff_stress
 
 
 class Results:
@@ -51,6 +52,7 @@ class Results:
     def __init__(self, stress=False, elasticity=False):
 
         self.force = None
+        self._force = None
         self.stiffness = None
         self.kinematics = None
 
diff --git a/felupe/mechanics/_solidbody.py b/felupe/mechanics/_solidbody.py
index 707eec51..56d22353 100644
--- a/felupe/mechanics/_solidbody.py
+++ b/felupe/mechanics/_solidbody.py
@@ -28,7 +28,7 @@
 from .._field import Field, FieldMixed, FieldsMixed, FieldAxisymmetric
 from .._assembly import IntegralForm, IntegralFormMixed, IntegralFormAxisymmetric
 from ..constitution import AreaChange
-from ..math import inv, dot
+from ..math import inv, dot, transpose, det
 from ._helpers import Assemble, Evaluate, Results
 
 
@@ -54,6 +54,7 @@ def __init__(self, umat, field):
             gradient=self._gradient,
             hessian=self._hessian,
             cauchy_stress=self._cauchy_stress,
+            kirchhoff_stress=self._kirchhoff_stress,
         )
 
         self._area_change = AreaChange()
@@ -141,6 +142,19 @@ def _hessian(self, field=None, args=(), kwargs={}):
 
         return self.results.elasticity
 
+    def _kirchhoff_stress(self, field=None):
+
+        self._gradient(field)
+
+        if len(self.results.kinematics) > 1:
+            P = self.results.stress[0]
+        else:
+            P = self.results.stress
+
+        F = self.results.kinematics[0]
+
+        return dot(P, transpose(F))
+
     def _cauchy_stress(self, field=None):
 
         self._gradient(field)
@@ -150,6 +164,7 @@ def _cauchy_stress(self, field=None):
         else:
             P = self.results.stress
 
-        JiFT = self._area_change.function(self.results.kinematics[0])
+        F = self.results.kinematics[0]
+        J = det(F)
 
-        return dot(P, inv(JiFT))
+        return dot(P, transpose(F)) / J
diff --git a/tests/test_math.py b/tests/test_math.py
index 7b70cc42..f40dad5d 100644
--- a/tests/test_math.py
+++ b/tests/test_math.py
@@ -127,6 +127,7 @@ def test_math():
     fe.math.cdya(F, F)
 
     fe.math.tovoigt(C)
+    fe.math.tovoigt(C[:2, :2])
     fe.math.eigvals(C)
     fe.math.eigvals(C[:2, :2])
     fe.math.eigvals(C, shear=True)
diff --git a/tests/test_mechanics.py b/tests/test_mechanics.py
index 1d9b2c62..19671c29 100644
--- a/tests/test_mechanics.py
+++ b/tests/test_mechanics.py
@@ -46,6 +46,7 @@ def test_simple():
     F = b.results.kinematics[0]
     P = b.results.stress
     s = b.evaluate.cauchy_stress()
+    t = b.evaluate.kirchhoff_stress()
     C = b.results.elasticity
 
     assert K.shape == (81, 81)
@@ -53,6 +54,7 @@ def test_simple():
     assert F.shape == (3, 3, 8, 8)
     assert P.shape == (3, 3, 8, 8)
     assert s.shape == (3, 3, 8, 8)
+    assert t.shape == (3, 3, 8, 8)
     assert C.shape == (3, 3, 3, 3, 8, 8)
 
 
@@ -193,6 +195,10 @@ def test_solidbody():
             s2 = b.evaluate.cauchy_stress(u)
             assert np.allclose(s1, s2)
 
+            t1 = b.evaluate.kirchhoff_stress()
+            t2 = b.evaluate.kirchhoff_stress(u)
+            assert np.allclose(t1, t2)
+
 
 def test_solidbody_axi():
 
@@ -234,6 +240,10 @@ def test_solidbody_axi():
             s2 = b.evaluate.cauchy_stress(u)
             assert np.allclose(s1, s2)
 
+            t1 = b.evaluate.kirchhoff_stress()
+            t2 = b.evaluate.kirchhoff_stress(u)
+            assert np.allclose(t1, t2)
+
 
 def test_solidbody_mixed():
 
@@ -280,6 +290,10 @@ def test_solidbody_mixed():
             s2 = b.evaluate.cauchy_stress(u)
             assert np.allclose(s1, s2)
 
+            t1 = b.evaluate.kirchhoff_stress()
+            t2 = b.evaluate.kirchhoff_stress(u)
+            assert np.allclose(t1, t2)
+
 
 if __name__ == "__main__":
     test_simple()