From 7b8028b9b60c08147374fd93b681a67c77e61246 Mon Sep 17 00:00:00 2001
From: Andreas D <a.dutzler@gmail.com>
Date: Thu, 29 Jul 2021 11:01:19 +0200
Subject: [PATCH 1/2] tools: strain -> strain_voigt, new strain function

---
 felupe/tools.py     |  7 ++++++-
 tests/test_tools.py | 15 ++++++++++++---
 2 files changed, 18 insertions(+), 4 deletions(-)

diff --git a/felupe/tools.py b/felupe/tools.py
index 5e1be8f4..1f349041 100644
--- a/felupe/tools.py
+++ b/felupe/tools.py
@@ -32,7 +32,7 @@
 
 import meshio
 
-from .math import grad, identity, interpolate, norms, dot, transpose, det, eigvals
+from .math import grad, sym, identity, interpolate, norms, dot, transpose, det, eigvals
 from . import solve as solvetools
 from .doftools import partition as dofpartition, apply
 from .forms import IntegralFormMixed
@@ -55,6 +55,11 @@ def defgrad(field):
 
 def strain(field):
     "Calculate strains of a given displacement field."
+    return sym(grad(field))
+
+
+def strain_voigt(field):
+    "Calculate strains in voigt notation of a given displacement field."
     dudX = grad(field)
     dim = dudX.shape[1]
 
diff --git a/tests/test_tools.py b/tests/test_tools.py
index 9a0671f0..7857ab04 100644
--- a/tests/test_tools.py
+++ b/tests/test_tools.py
@@ -63,12 +63,20 @@ def test_defgrad():
 
 def test_strain():
 
+    r, fields = pre()
+    eps = fe.tools.strain(fields[0])
+
+    assert eps.shape == (3, 3, r.quadrature.npoints, r.mesh.ncells)
+
+
+def test_strain_voigt():
+
     m = fe.mesh.Line(n=3)
     e = fe.element.Line()
     q = fe.quadrature.GaussLegendre(1, 1)
     r = fe.Region(m, e, q)
     u = fe.Field(r, dim=1)
-    strain = fe.tools.strain(u)
+    strain = fe.tools.strain_voigt(u)
     assert strain.shape == (1, r.quadrature.npoints, r.mesh.ncells)
 
     m = fe.mesh.Rectangle(n=3)
@@ -76,7 +84,7 @@ def test_strain():
     q = fe.quadrature.GaussLegendre(1, 2)
     r = fe.Region(m, e, q)
     u = fe.Field(r, dim=2)
-    strain = fe.tools.strain(u)
+    strain = fe.tools.strain_voigt(u)
     assert strain.shape == (3, r.quadrature.npoints, r.mesh.ncells)
 
     m = fe.mesh.Cube(n=3)
@@ -84,7 +92,7 @@ def test_strain():
     q = fe.quadrature.GaussLegendre(1, 3)
     r = fe.Region(m, e, q)
     u = fe.Field(r, dim=3)
-    strain = fe.tools.strain(u)
+    strain = fe.tools.strain_voigt(u)
     assert strain.shape == (6, r.quadrature.npoints, r.mesh.ncells)
 
 
@@ -170,6 +178,7 @@ def jac(x):
     test_extract()
     test_defgrad()
     test_strain()
+    test_strain_voigt()
     test_update()
     test_solve_check()
     test_solve_mixed_check()

From 4a2b8e60d00b01a7ed386ea423d2284e7e492bba Mon Sep 17 00:00:00 2001
From: Andreas D <a.dutzler@gmail.com>
Date: Thu, 29 Jul 2021 11:03:25 +0200
Subject: [PATCH 2/2] constition: fix linearelastic elasticity

add initial stress elasticitiy if stress is passed
---
 felupe/constitution/models.py | 12 ++++++++++--
 1 file changed, 10 insertions(+), 2 deletions(-)

diff --git a/felupe/constitution/models.py b/felupe/constitution/models.py
index b3117029..eff29d69 100644
--- a/felupe/constitution/models.py
+++ b/felupe/constitution/models.py
@@ -54,9 +54,17 @@ def __init__(self, E, nu):
     def stress(self, strain):
         return 2 * self.mu * strain + self.gamma * trace(strain) * identity(strain)
 
-    def elasticity(self, strain):
+    def elasticity(self, strain, stress=None):
         I = identity(strain)
-        return 2 * self.mu * cdya(I, I) + self.gamma * dya(I, I)
+
+        elast = 2 * self.mu * cdya(I, I) + self.gamma * dya(I, I)
+
+        if stress is not None:
+            elast_stress = cdya_ik(stress, I)
+        else:
+            elast_stress = 0
+
+        return elast + elast_stress
 
     def lame(self, E, nu):
         mu = E / (2 * (1 + nu))