From 9e7b291f8fe12528817ba9b6e5f800b5acb5d051 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Fri, 19 Aug 2022 11:50:01 +0200 Subject: [PATCH 1/3] field integration: fix field-trimming for plane strain fixes #259 --- felupe/_assembly/_base.py | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/felupe/_assembly/_base.py b/felupe/_assembly/_base.py index 9ec97e87..13f96115 100644 --- a/felupe/_assembly/_base.py +++ b/felupe/_assembly/_base.py @@ -163,8 +163,15 @@ def integrate(self, parallel=False, jit=False): dV = self.dV fun = self.fun - if len(fun) > v.dim: - fun = fun[tuple([slice(v.dim)] * (len(fun.shape) - 2))] + # plane strain + # trim 3d vector-valued functions to the dimension of the field + function_dimension = len(fun.shape) - 2 + function_is_vector = function_dimension >= 1 + function_is_3d = len(fun) == 3 + field_is_2d = v.dim == 2 + + if function_is_vector and function_is_3d and field_is_2d: + fun = fun[tuple([slice(2)] * function_dimension)] if parallel: einsum = einsumt From f86ef0090cfbe4efdfcd842caec8280def7e5d6c Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Fri, 19 Aug 2022 12:23:49 +0200 Subject: [PATCH 2/3] fix stress output of solid body tensor --- felupe/mechanics/_solidbody_tensor.py | 5 +- tests/test_mechanics_tensor.py | 125 ++++++++++++++++++-------- 2 files changed, 89 insertions(+), 41 deletions(-) diff --git a/felupe/mechanics/_solidbody_tensor.py b/felupe/mechanics/_solidbody_tensor.py index ea1a2983..059e1fb2 100644 --- a/felupe/mechanics/_solidbody_tensor.py +++ b/felupe/mechanics/_solidbody_tensor.py @@ -77,7 +77,6 @@ def _vector( self.field = field self.results.stress = self._gradient(field, args=args, kwargs=kwargs) - self.results.force = self._form( fun=self.results.stress[slice(items)], v=self.field, @@ -141,7 +140,7 @@ def _kirchhoff_stress(self, field=None): self._gradient(field) - P = self.results.stress[0][0] + P = self.results.stress[0] F = self.results.kinematics[0] return dot(P, transpose(F)) @@ -150,7 +149,7 @@ def _cauchy_stress(self, field=None): self._gradient(field) - P = self.results.stress[0][0] + P = self.results.stress[0] F = self.results.kinematics[0] J = det(F) diff --git a/tests/test_mechanics_tensor.py b/tests/test_mechanics_tensor.py index 9a3b523d..627bc61a 100644 --- a/tests/test_mechanics_tensor.py +++ b/tests/test_mechanics_tensor.py @@ -32,23 +32,43 @@ def pre_umat(): - LE = fe.LinearElastic(E=210000, nu=0.3) + NH = fe.NeoHooke(mu=1, bulk=5000) - class LETensor: - def __init__(self, LE): - self.LE = LE - self.x = [np.zeros((3, 3)), np.zeros((5, 16))] + class NHTensor: + def __init__(self, NH): + self.NH = NH + self.x = [np.zeros((3, 3)), np.zeros((1, 1))] def function(self, x): - F, statevars = x + F, statevars = x[:-1], x[-1] # return dummy state variables along with stress - return self.LE.stress(F), statevars + return [*self.NH.gradient(F), statevars] def gradient(self, x): - F, statevars = x - return self.LE.elasticity(F) + F, statevars = x[:-1], x[-1] + return self.NH.hessian(F) - return LETensor(LE) + return NHTensor(NH) + +def pre_umat_mixed(): + + NH = fe.ThreeFieldVariation(fe.NeoHooke(mu=1, bulk=5000)) + + class NHTensor: + def __init__(self, NH): + self.NH = NH + self.x = [np.zeros((3, 3)), np.zeros((1, 1)), np.zeros((1, 1)), np.zeros((1, 1))] + + def function(self, x): + F, statevars = x[:-1], x[-1] + # return dummy state variables along with stress + return [*self.NH.gradient(F), statevars] + + def gradient(self, x): + F, statevars = x[:-1], x[-1] + return self.NH.hessian(F) + + return NHTensor(NH) def test_solidbody_tensor(): @@ -57,35 +77,64 @@ def test_solidbody_tensor(): m = fe.Cube(n=3) r = fe.RegionHexahedron(m) - u = fe.Field(r, dim=3) - v = fe.FieldContainer([u]) - - sv = np.zeros((5, 16, r.quadrature.npoints, m.ncells)) - - for statevars in [sv, None]: - b = fe.SolidBodyTensor(umat, v, statevars) - r = b.assemble.vector() - - K = b.assemble.matrix(v) - K = b.assemble.matrix() - r = b.assemble.vector(v) - r = b.assemble.vector() - F = b.results.kinematics - P = b.results.stress - s = b.evaluate.cauchy_stress() - t = b.evaluate.kirchhoff_stress() - C = b.results.elasticity - z = b.results.statevars - - assert K.shape == (81, 81) - assert r.shape == (81, 1) - assert F[0].shape == (3, 3, 8, 8) - assert P[0][0].shape == (3, 3, 8, 8) - assert s.shape == (3, 3, 8, 8) - assert t.shape == (3, 3, 8, 8) - assert C[0].shape == (3, 3, 3, 3, 8, 8) - assert z.shape == (5, 16, 8, 8) + v = fe.FieldsMixed(r, n=1) + + b = fe.SolidBodyTensor(umat, v) + r = b.assemble.vector() + + K = b.assemble.matrix(v) + K = b.assemble.matrix() + r = b.assemble.vector(v) + r = b.assemble.vector() + F = b.results.kinematics + P = b.results.stress + s = b.evaluate.cauchy_stress() + t = b.evaluate.kirchhoff_stress() + C = b.results.elasticity + z = b.results.statevars + + assert K.shape == (81, 81) + assert r.shape == (81, 1) + assert F[0].shape == (3, 3, 8, 8) + assert P[0].shape == (3, 3, 8, 8) + assert s.shape == (3, 3, 8, 8) + assert t.shape == (3, 3, 8, 8) + assert C[0].shape == (3, 3, 3, 3, 8, 8) + assert z.shape == (1, 1, 8, 8) + + +def test_solidbody_tensor_mixed(): + + umat = pre_umat_mixed() + + m = fe.Cube(n=3) + r = fe.RegionHexahedron(m) + v = fe.FieldsMixed(r, n=3) + + b = fe.SolidBodyTensor(umat, v) + r = b.assemble.vector() + + K = b.assemble.matrix(v) + K = b.assemble.matrix() + r = b.assemble.vector(v) + r = b.assemble.vector() + F = b.results.kinematics + P = b.results.stress + s = b.evaluate.cauchy_stress() + t = b.evaluate.kirchhoff_stress() + C = b.results.elasticity + z = b.results.statevars + + assert K.shape == (97, 97) + assert r.shape == (97, 1) + assert F[0].shape == (3, 3, 8, 8) + assert P[0].shape == (3, 3, 8, 8) + assert s.shape == (3, 3, 8, 8) + assert t.shape == (3, 3, 8, 8) + assert C[0].shape == (3, 3, 3, 3, 8, 8) + assert z.shape == (1, 1, 8, 8) if __name__ == "__main__": test_solidbody_tensor() + test_solidbody_tensor_mixed() \ No newline at end of file From b43c371044ba2628aefb5ad83381038e53d8dc49 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Fri, 19 Aug 2022 12:27:55 +0200 Subject: [PATCH 3/3] solid body tensor: test user statevars --- tests/test_mechanics_tensor.py | 98 ++++++++++++++++++---------------- 1 file changed, 53 insertions(+), 45 deletions(-) diff --git a/tests/test_mechanics_tensor.py b/tests/test_mechanics_tensor.py index 627bc61a..f0688a64 100644 --- a/tests/test_mechanics_tensor.py +++ b/tests/test_mechanics_tensor.py @@ -78,29 +78,33 @@ def test_solidbody_tensor(): m = fe.Cube(n=3) r = fe.RegionHexahedron(m) v = fe.FieldsMixed(r, n=1) - - b = fe.SolidBodyTensor(umat, v) - r = b.assemble.vector() - - K = b.assemble.matrix(v) - K = b.assemble.matrix() - r = b.assemble.vector(v) - r = b.assemble.vector() - F = b.results.kinematics - P = b.results.stress - s = b.evaluate.cauchy_stress() - t = b.evaluate.kirchhoff_stress() - C = b.results.elasticity - z = b.results.statevars - - assert K.shape == (81, 81) - assert r.shape == (81, 1) - assert F[0].shape == (3, 3, 8, 8) - assert P[0].shape == (3, 3, 8, 8) - assert s.shape == (3, 3, 8, 8) - assert t.shape == (3, 3, 8, 8) - assert C[0].shape == (3, 3, 3, 3, 8, 8) - assert z.shape == (1, 1, 8, 8) + + sv = np.zeros((1, 1, r.quadrature.npoints, m.ncells)) + + for statevars in [sv, None]: + + b = fe.SolidBodyTensor(umat, v, statevars) + r = b.assemble.vector() + + K = b.assemble.matrix(v) + K = b.assemble.matrix() + r = b.assemble.vector(v) + r = b.assemble.vector() + F = b.results.kinematics + P = b.results.stress + s = b.evaluate.cauchy_stress() + t = b.evaluate.kirchhoff_stress() + C = b.results.elasticity + z = b.results.statevars + + assert K.shape == (81, 81) + assert r.shape == (81, 1) + assert F[0].shape == (3, 3, 8, 8) + assert P[0].shape == (3, 3, 8, 8) + assert s.shape == (3, 3, 8, 8) + assert t.shape == (3, 3, 8, 8) + assert C[0].shape == (3, 3, 3, 3, 8, 8) + assert z.shape == (1, 1, 8, 8) def test_solidbody_tensor_mixed(): @@ -111,28 +115,32 @@ def test_solidbody_tensor_mixed(): r = fe.RegionHexahedron(m) v = fe.FieldsMixed(r, n=3) - b = fe.SolidBodyTensor(umat, v) - r = b.assemble.vector() - - K = b.assemble.matrix(v) - K = b.assemble.matrix() - r = b.assemble.vector(v) - r = b.assemble.vector() - F = b.results.kinematics - P = b.results.stress - s = b.evaluate.cauchy_stress() - t = b.evaluate.kirchhoff_stress() - C = b.results.elasticity - z = b.results.statevars - - assert K.shape == (97, 97) - assert r.shape == (97, 1) - assert F[0].shape == (3, 3, 8, 8) - assert P[0].shape == (3, 3, 8, 8) - assert s.shape == (3, 3, 8, 8) - assert t.shape == (3, 3, 8, 8) - assert C[0].shape == (3, 3, 3, 3, 8, 8) - assert z.shape == (1, 1, 8, 8) + sv = np.zeros((1, 1, r.quadrature.npoints, m.ncells)) + + for statevars in [sv, None]: + + b = fe.SolidBodyTensor(umat, v, statevars) + r = b.assemble.vector() + + K = b.assemble.matrix(v) + K = b.assemble.matrix() + r = b.assemble.vector(v) + r = b.assemble.vector() + F = b.results.kinematics + P = b.results.stress + s = b.evaluate.cauchy_stress() + t = b.evaluate.kirchhoff_stress() + C = b.results.elasticity + z = b.results.statevars + + assert K.shape == (97, 97) + assert r.shape == (97, 1) + assert F[0].shape == (3, 3, 8, 8) + assert P[0].shape == (3, 3, 8, 8) + assert s.shape == (3, 3, 8, 8) + assert t.shape == (3, 3, 8, 8) + assert C[0].shape == (3, 3, 3, 3, 8, 8) + assert z.shape == (1, 1, 8, 8) if __name__ == "__main__":