diff --git a/trunk/constitutive_pheno.f90 b/trunk/constitutive_pheno.f90 index c6346ebe0..b3f3224b6 100644 --- a/trunk/constitutive_pheno.f90 +++ b/trunk/constitutive_pheno.f90 @@ -751,6 +751,8 @@ subroutine constitutive_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,state,Temperature, !********************************************************************* use prec, only: pReal,pInt use crystal, only: crystal_Sslip,crystal_Sslip_v +use math, only: math_Plain3333to99 + implicit none !* Definition of variables @@ -759,7 +761,8 @@ integer(pInt) matID,i,k,l,m,n real(pReal) Temperature real(pReal), dimension(6) :: Tstar_v real(pReal), dimension(3,3) :: Lp -real(pReal), dimension(3,3,3,3) :: dLp_dTstar +real(pReal), dimension(3,3,3,3) :: dLp +real(pReal), dimension(9,9) :: dLp_dTstar real(pReal), dimension(constitutive_Nstatevars(ipc,ip,el)) :: state real(pReal), dimension(material_Nslip(constitutive_matID(ipc,ip,el))) :: gdot_slip,dgdot_dtauslip,tau_slip @@ -776,17 +779,17 @@ do i=1,material_Nslip(matID) enddo !* Calculation of the tangent of Lp -dLp_dTstar=0.0_pReal +dLp = 0.0_pReal +dLp_dTstar = 0.0_pReal do i=1,material_Nslip(matID) dgdot_dtauslip(i) = material_gdot0_slip(matID)*(abs(tau_slip(i))/state(i))**& (material_n_slip(matID)-1.0_pReal)*material_n_slip(matID)/state(i) - forall (k=1:3,l=1:3,m=1:3,n=1:3) - dLp_dTstar(k,l,m,n) = dLp_dTstar(k,l,m,n)+ & - dgdot_dtauslip(i)*crystal_Sslip(k,l,i,material_CrystalStructure(matID))* & - (crystal_Sslip(m,n,i,material_CrystalStructure(matID))+ & - crystal_Sslip(n,m,i,material_CrystalStructure(matID)))/2.0_pReal ! force m,n symmetry - endforall + forall (k=1:3,l=1:3,m=1:3,n=1:3) & + dLp(k,l,m,n) = dLp(k,l,m,n) + & + dgdot_dtauslip(i)*crystal_Sslip(k,l,i,material_CrystalStructure(matID))* & + crystal_Sslip(m,n,i,material_CrystalStructure(matID)) enddo +dLp_dTstar = math_Plain3333to99(dLp) return end subroutine @@ -803,7 +806,7 @@ function constitutive_dotState(Tstar_v,state,Temperature,ipc,ip,el) !* - ip : current integration point * !* - el : current element * !* OUTPUT: * -!* - constitutive_DotState : evolution of state variable * +!* - constitutive_dotState : evolution of state variable * !********************************************************************* use prec, only: pReal,pInt use crystal, only: crystal_Sslip_v diff --git a/trunk/crystal.f90 b/trunk/crystal.f90 index 618f8df4b..ba02cbcd1 100644 --- a/trunk/crystal.f90 +++ b/trunk/crystal.f90 @@ -305,7 +305,7 @@ subroutine crystal_SchmidMatrices() !* Calculation of Schmid matrices * !************************************** use prec, only: pReal,pInt -use math, only: math_I3 +use math, only: math_I3,nrmMandel,mapMandel implicit none !* Definition of variables @@ -327,17 +327,10 @@ do l=1,crystal_MaxCrystalStructure crystal_st(:,k,l)=crystal_st(:,k,l)/norm_t crystal_sn(:,k,l)=crystal_sn(:,k,l)/norm_n !* Defintion of Schmid matrix - forall (i=1:3,j=1:3) - crystal_Sslip(i,j,k,l)=crystal_sd(i,k,l)*crystal_sn(j,k,l) - endforall + forall (i=1:3,j=1:3) crystal_Sslip(i,j,k,l)=crystal_sd(i,k,l)*crystal_sn(j,k,l) !* Vectorization of normalized Schmid matrix - crystal_Sslip_v(1,k,l)=crystal_Sslip(1,1,k,l) - crystal_Sslip_v(2,k,l)=crystal_Sslip(2,2,k,l) - crystal_Sslip_v(3,k,l)=crystal_Sslip(3,3,k,l) - !* be compatible with Mandel notation of Tstar - crystal_Sslip_v(4,k,l)=(crystal_Sslip(1,2,k,l)+crystal_Sslip(2,1,k,l))/dsqrt(2.0_pReal) - crystal_Sslip_v(5,k,l)=(crystal_Sslip(2,3,k,l)+crystal_Sslip(3,2,k,l))/dsqrt(2.0_pReal) - crystal_Sslip_v(6,k,l)=(crystal_Sslip(1,3,k,l)+crystal_Sslip(3,1,k,l))/dsqrt(2.0_pReal) + forall (i=1:6) crystal_Sslip_v(i,k,l) = nrmMandel(i)/2.0_pReal * & + (crystal_Sslip(mapMandel(1,i),mapMandel(2,i),k,l)+crystal_Sslip(mapMandel(2,i),mapMandel(1,i),k,l)) enddo !* Iteration over the twin systems diff --git a/trunk/debug.f90 b/trunk/debug.f90 index 2fd50fc39..4de13b164 100644 --- a/trunk/debug.f90 +++ b/trunk/debug.f90 @@ -2,49 +2,50 @@ !############################################################## MODULE debug !############################################################## - use prec + use prec + implicit none - integer(pInt), dimension(nCutback+1) :: debug_cutbackDistribution - integer(pInt), dimension(nInner) :: debug_innerLoopDistribution - integer(pInt), dimension(nOuter) :: debug_outerLoopDistribution - logical debugger - - CONTAINS - - -!******************************************************************** -! write debug statements to standard out -!******************************************************************** - SUBROUTINE debug_info() - - use prec - implicit none - - integer(pInt) i - - write(6,*) 'DEBUG Info' - write(6,*) 'distribution_cutback :' - do i=0,nCutback - if (debug_cutbackDistribution(i+1) > 0) write(6,*) i,debug_cutbackDistribution(i+1) - enddo - write(6,*) 'total',sum(debug_cutbackDistribution) - write(6,*) - - write(6,*) 'distribution_innerLoop :' - do i=1,nInner - if (debug_innerLoopDistribution(i) > 0) write(6,*) i,debug_innerLoopDistribution(i) - enddo - write(6,*) 'total',sum(debug_innerLoopDistribution) - write(6,*) - - write(6,*) 'distribution_outerLoop :' - do i=1,nOuter - if (debug_outerLoopDistribution(i) > 0) write(6,*) i,debug_outerLoopDistribution(i) - enddo - write(6,*) 'total',sum(debug_outerLoopDistribution) - write(6,*) - - END SUBROUTINE + integer(pInt), dimension(nCutback+1) :: debug_cutbackDistribution + integer(pInt), dimension(nInner) :: debug_InnerLoopDistribution + integer(pInt), dimension(nOuter) :: debug_OuterLoopDistribution + logical :: debugger = .false. + + CONTAINS + + +!******************************************************************** +! write debug statements to standard out +!******************************************************************** + SUBROUTINE debug_info() + + use prec + implicit none + + integer(pInt) i + + write(6,*) 'DEBUG Info' + write(6,*) 'distribution_cutback :' + do i=0,nCutback + if (debug_cutbackDistribution(i+1) /= 0) write(6,*) i,debug_cutbackDistribution(i+1) + enddo + write(6,*) 'total',sum(debug_cutbackDistribution) + write(6,*) + + write(6,*) 'distribution_InnerLoop :' + do i=1,nInner + if (debug_InnerLoopDistribution(i) /= 0) write(6,*) i,debug_InnerLoopDistribution(i) + enddo + write(6,*) 'total',sum(debug_InnerLoopDistribution) + write(6,*) + + write(6,*) 'distribution_OuterLoop :' + do i=1,nOuter + if (debug_OuterLoopDistribution(i) /= 0) write(6,*) i,debug_OuterLoopDistribution(i) + enddo + write(6,*) 'total',sum(debug_OuterLoopDistribution) + write(6,*) + + END SUBROUTINE END MODULE debug diff --git a/trunk/mpie_cpfem_marc2005r3.f90 b/trunk/mpie_cpfem_marc2005r3.f90 index ec425eaf0..7364ade17 100644 --- a/trunk/mpie_cpfem_marc2005r3.f90 +++ b/trunk/mpie_cpfem_marc2005r3.f90 @@ -4,7 +4,7 @@ ! written by F. Roters, P. Eisenlohr, L. Hantcherli, W.A. Counts ! MPI fuer Eisenforschung, Duesseldorf ! -! last modified: 18.02.2005 +! last modified: 08.11.2007 !******************************************************************** ! Usage: ! - choose material as hypela2 @@ -153,10 +153,6 @@ ! dimension e(*),de(*),t(*),dt(*),g(*),d(ngens,*),s(*), n(2),coord(ncrd,*),disp(ndeg,*),matus(2),dispt(ndeg,*),ffn(itel,*),& frotn(itel,*),strechn(itel),eigvn(itel,*),ffn1(itel,*),frotn1(itel,*),strechn1(itel),eigvn1(itel,*),kcus(2) - - logical stress_recovery - - stress_recovery = (lovl == 6) ! ! subroutine cpfem_general(mpie_ffn, mpie_ffn1, temperature, mpie_inc, mpie_subinc, mpie_cn, ! mpie_stress_recovery, mpie_tinc, mpie_en, mpie_in, mpie_s, mpie_d, mpie_ngens) @@ -178,10 +174,6 @@ ! mpie_ngens size of stress strain law !******************************************************************** call CPFEM_general(ffn, ffn1, t(1), inc, incsub, ncycle, stress_recovery, timinc, n(1), nn, s, d, ngens) -! Mandel: 11, 22, 33, SQRT(2)*12, SQRT(2)*23, SQRT(2)*13 -! Marc: 11, 22, 33, 12, 23, 13 - forall(i=1:ngens) d(1:ngens,i) = invnrmMandel(i)*d(1:ngens,i)*invnrmMandel(1:ngens) - s(1:ngens) = s(1:ngens)*invnrmMandel(1:ngens) return END SUBROUTINE diff --git a/trunk/mpie_cpfem_marc2007r1.f90 b/trunk/mpie_cpfem_marc2007r1.f90 index 003a6c740..6d249dace 100644 --- a/trunk/mpie_cpfem_marc2007r1.f90 +++ b/trunk/mpie_cpfem_marc2007r1.f90 @@ -28,7 +28,9 @@ !******************************************************************** ! include "prec.f90" + include "debug.f90" + include "math.f90" include "IO.f90" include "mesh.f90" @@ -195,6 +197,7 @@ ! Marc: 11, 22, 33, 12, 23, 13 forall(i=1:ngens) d(1:ngens,i) = invnrmMandel(i)*d(1:ngens,i)*invnrmMandel(1:ngens) s(1:ngens) = s(1:ngens)*invnrmMandel(1:ngens) + return END SUBROUTINE diff --git a/trunk/prec.f90 b/trunk/prec.f90 index d6c2ff6c7..371ab13b3 100644 --- a/trunk/prec.f90 +++ b/trunk/prec.f90 @@ -17,7 +17,7 @@ integer(pInt), parameter :: nReg = 1_pInt ! regularization attempts for Jacobi inversion real(pReal), parameter :: pert_Fg = 1.0e-5_pReal ! strain perturbation for FEM Jacobi integer(pInt), parameter :: nOuter = 10_pInt ! outer loop limit - integer(pInt), parameter :: nInner = 1000_pInt ! inner loop limit + integer(pInt), parameter :: nInner = 200_pInt ! inner loop limit real(pReal), parameter :: reltol_Outer = 1.0e-4_pReal ! relative tolerance in outer loop (state) real(pReal), parameter :: reltol_Inner = 1.0e-6_pReal ! relative tolerance in inner loop (Lp) real(pReal), parameter :: abstol_Inner = 1.0e-8_pReal ! absolute tolerance in inner loop (Lp)