From 3dfdbaff5ba132d75528466353236d0e17cb7cac Mon Sep 17 00:00:00 2001 From: Christoph Kords Date: Wed, 22 Jan 2014 08:38:13 +0000 Subject: [PATCH] Fixed wrong indices in tangents dT_dFe and dFe_dLp, which however luckily did not have any effect in the perturbed stiffness since they were transposed such that the double contraction of both remained unchanged. In contrast, the analytical jacobian will probably be affected by this change! @Pratheek: Can you check with me how this can be fixed? --- code/constitutive.f90 | 9 ++++++--- code/crystallite.f90 | 2 +- 2 files changed, 7 insertions(+), 4 deletions(-) diff --git a/code/constitutive.f90 b/code/constitutive.f90 index dee1e00f8..4ec21cb09 100644 --- a/code/constitutive.f90 +++ b/code/constitutive.f90 @@ -658,6 +658,7 @@ end subroutine constitutive_TandItsTangent !-------------------------------------------------------------------------------------------------- pure subroutine constitutive_hooke_TandItsTangent(T, dT_dFe, Fe, ipc, ip, el) use math, only : & + math_mul3x3, & math_mul33x33, & math_mul3333xx33, & math_Mandel66to3333, & @@ -676,7 +677,7 @@ use math, only : & real(pReal), intent(out), dimension(3,3,3,3) :: & dT_dFe !< dT/dFe - integer(pInt) :: o, p + integer(pInt) :: i, j, k, l real(pReal), dimension(3,3) :: FeT real(pReal), dimension(3,3,3,3) :: C @@ -684,9 +685,11 @@ use math, only : & C = math_Mandel66to3333(constitutive_homogenizedC(ipc,ip,el)) FeT = math_transpose33(Fe) - T = 0.5_pReal*math_mul3333xx33(C,math_mul33x33(FeT,Fe)-MATH_I3) + T = 0.5_pReal * math_mul3333xx33(C, math_mul33x33(FeT,Fe)-MATH_I3) - forall (o=1_pInt:3_pInt, p=1_pInt:3_pInt) dT_dFe(o,p,1:3,1:3) = math_mul33x33(C(o,p,1:3,1:3), FeT) ! dT*_ij/dFe_kl + dT_dFe = 0.0_pReal + forall (i=1_pInt:3_pInt, j=1_pInt:3_pInt, k=1_pInt:3_pInt, l=1_pInt:3_pInt) & + dT_dFe(i,j,k,l) = math_mul3x3(C(i,j,l,1:3),Fe(k,1:3)) ! dT*_ij/dFe_kl end subroutine constitutive_hooke_TandItsTangent diff --git a/code/crystallite.f90 b/code/crystallite.f90 index 40e85b532..6c8eda57d 100644 --- a/code/crystallite.f90 +++ b/code/crystallite.f90 @@ -3211,7 +3211,7 @@ logical function crystallite_integrateStress(& if (mod(jacoCounter, iJacoLpresiduum) == 0_pInt) then dFe_dLp3333 = 0.0_pReal do o=1_pInt,3_pInt; do p=1_pInt,3_pInt - dFe_dLp3333(p,o,1:3,p) = A(o,1:3) ! dFe_dLp(i,j,k,l) = -dt * A(i,k) delta(l,j) + dFe_dLp3333(o,p,1:3,p) = A(o,1:3) ! dFe_dLp(i,j,k,l) = -dt * A(i,k) delta(l,j) enddo; enddo dFe_dLp3333 = -dt * dFe_dLp3333 dFe_dLp = math_Plain3333to99(dFe_dLp3333)