Merge pull request #60 from adtzlr/fix-linearelastic

Fix linearelastic

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))

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]
 

tests/test_tools.py

@@ -63,28 +63,36 @@ 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)
     e = fe.element.Quad()
     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)
     e = fe.element.Hexahedron()
     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()