diff --git a/trunk/CPFEM_GIA8.f90 b/trunk/CPFEM_GIA8.f90 new file mode 100644 index 000000000..52747fabe --- /dev/null +++ b/trunk/CPFEM_GIA8.f90 @@ -0,0 +1,798 @@ +!############################################################## + MODULE CPFEM +!############################################################## +! *** CPFEM engine *** +! + use prec, only: pReal,pInt + implicit none +! +! **************************************************************** +! *** General variables for the material behaviour calculation *** +! **************************************************************** + real(pReal), dimension (:,:), allocatable :: CPFEM_Temperature + real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_ffn_bar + real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_ffn1_bar + real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_PK1_bar + real(pReal), dimension (:,:,:,:,:,:),allocatable :: CPFEM_dPdF_bar + real(pReal), dimension (:,:,:), allocatable :: CPFEM_stress_bar + real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_jaco_bar + real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_jaco_knownGood + real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_results + real(pReal), dimension (:,:,:,:,:), allocatable :: CPFEM_Fp_old + real(pReal), dimension (:,:,:,:,:), allocatable :: CPFEM_Fp_new + real(pReal), parameter :: CPFEM_odd_stress = 1e15_pReal, CPFEM_odd_jacobian = 1e50_pReal + integer(pInt) :: CPFEM_Nresults = 4_pInt ! three Euler angles plus volume fraction + logical :: CPFEM_init_done = .false. ! remember if init has been done already + logical :: CPFEM_calc_done = .false. ! remember if first IP has already calced the results +! + real(pReal), dimension (:,:,:,:), allocatable :: GIA_rVect_new ! boundary relaxation vectors + real(pReal), dimension (:,:,:,:), allocatable :: GIA_rVect_old ! boundary relaxation vectors + real(pReal), dimension (:,:), allocatable :: GIA_bNorm ! grain boundary normals +! + CONTAINS +! +!********************************************************* +!*** allocate the arrays defined in module CPFEM *** +!*** and initialize them *** +!********************************************************* + SUBROUTINE CPFEM_init(Temperature) +! + use prec + use math, only: math_EulertoR, math_I3, math_identity2nd + use mesh + use constitutive +! + implicit none +! + real(pReal) Temperature + integer(pInt) e,i,g,b +! +! *** mpie.marc parameters *** + allocate(CPFEM_Temperature (mesh_maxNips,mesh_NcpElems)) ; CPFEM_Temperature = Temperature + allocate(CPFEM_ffn_bar (3,3,mesh_maxNips,mesh_NcpElems)) + forall(e=1:mesh_NcpElems,i=1:mesh_maxNips) CPFEM_ffn_bar(:,:,i,e) = math_I3 + allocate(CPFEM_ffn1_bar (3,3,mesh_maxNips,mesh_NcpElems)) ; CPFEM_ffn1_bar = CPFEM_ffn_bar + allocate(CPFEM_PK1_bar (3,3,mesh_maxNips,mesh_NcpElems)) ; CPFEM_PK1_bar = 0.0_pReal + allocate(CPFEM_dPdF_bar(3,3,3,3,mesh_maxNips,mesh_NcpElems)) ; CPFEM_dPdF_bar = 0.0_pReal + allocate(CPFEM_stress_bar(6,mesh_maxNips,mesh_NcpElems)) ; CPFEM_stress_bar = 0.0_pReal + allocate(CPFEM_jaco_bar(6,6,mesh_maxNips,mesh_NcpElems)) ; CPFEM_jaco_bar = 0.0_pReal + allocate(CPFEM_jaco_knownGood(6,6,mesh_maxNips,mesh_NcpElems)) ; CPFEM_jaco_knownGood = 0.0_pReal +! +! *** User defined results !!! MISSING incorporate consti_Nresults *** + allocate(CPFEM_results(CPFEM_Nresults+constitutive_maxNresults,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems)) + CPFEM_results = 0.0_pReal +! +! *** Plastic deformation gradient at (t=t0) and (t=t1) *** + allocate(CPFEM_Fp_new(3,3,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; CPFEM_Fp_new = 0.0_pReal + allocate(CPFEM_Fp_old(3,3,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems)) + forall (e=1:mesh_NcpElems,i=1:mesh_maxNips,g=1:constitutive_maxNgrains) & + CPFEM_Fp_old(:,:,g,i,e) = math_EulerToR(constitutive_EulerAngles(:,g,i,e)) ! plastic def gradient reflects init orientation +! + allocate(GIA_rVect_new(3,12,mesh_maxNips,mesh_NcpElems)) ; GIA_rVect_new = 0.0_pReal + allocate(GIA_rVect_old(3,12,mesh_maxNips,mesh_NcpElems)) ; GIA_rVect_old = 0.0_pReal + allocate(GIA_bNorm(3,12)) ; GIA_bNorm = 0.0_pReal + do b = 1,4 + GIA_bNorm(1,b) = 1.0_pReal + GIA_bNorm(2,b+4) = 1.0_pReal + GIA_bNorm(3,b+8) = 1.0_pReal + enddo +! +! *** Output to MARC output file *** + write(6,*) + write(6,*) 'CPFEM Initialization' + write(6,*) + write(6,*) 'CPFEM_Temperature: ', shape(CPFEM_Temperature) + write(6,*) 'CPFEM_ffn_bar: ', shape(CPFEM_ffn_bar) + write(6,*) 'CPFEM_ffn1_bar: ', shape(CPFEM_ffn1_bar) + write(6,*) 'CPFEM_PK1_bar: ', shape(CPFEM_PK1_bar) + write(6,*) 'CPFEM_dPdF_bar: ', shape(CPFEM_dPdF_bar) + write(6,*) 'CPFEM_stress_bar: ', shape(CPFEM_stress_bar) + write(6,*) 'CPFEM_jaco_bar: ', shape(CPFEM_jaco_bar) + write(6,*) 'CPFEM_jaco_knownGood: ', shape(CPFEM_jaco_knownGood) + write(6,*) 'CPFEM_results: ', shape(CPFEM_results) + write(6,*) 'CPFEM_Fp_old: ', shape(CPFEM_Fp_old) + write(6,*) 'CPFEM_Fp_new: ', shape(CPFEM_Fp_new) +! + write(6,*) 'GIA_rVect_new: ', shape(GIA_rVect_new) + write(6,*) 'GIA_rVect_old: ', shape(GIA_rVect_old) + write(6,*) 'GIA_bNorm: ', shape(GIA_bNorm) + write(6,*) + call flush(6) + return +! + END SUBROUTINE +! +! +!*********************************************************************** +!*** perform initialization at first call, update variables and *** +!*** call the actual material model *** +! +! CPFEM_mode computation mode (regular, collection, recycle) +! ffn deformation gradient for t=t0 +! ffn1 deformation gradient for t=t1 +! Temperature temperature +! CPFEM_dt time increment +! CPFEM_en element number +! CPFEM_in intergration point number +! CPFEM_stress stress vector in Mandel notation +! CPFEM_updateJaco flag to initiate computation of Jacobian +! CPFEM_jaco jacobian in Mandel notation +! CPFEM_ngens size of stress strain law +!*********************************************************************** + SUBROUTINE CPFEM_general(CPFEM_mode, ffn, ffn1, Temperature, CPFEM_dt,& + CPFEM_en, CPFEM_in, CPFEM_stress, CPFEM_updateJaco, CPFEM_jaco, CPFEM_ngens) +! note: CPFEM_stress = Cauchy stress cs(6) and CPFEM_jaco = Consistent tangent dcs/de +! + use prec, only: pReal,pInt + use FEsolving + use debug + use math, only: math_init, invnrmMandel, math_identity2nd, math_Mandel3333to66,math_Mandel33to6,math_Mandel6to33,math_det3x3,math_I3 + use mesh, only: mesh_init,mesh_FEasCP, mesh_NcpElems, FE_Nips, FE_mapElemtype, mesh_element + use lattice, only: lattice_init + use constitutive, only: constitutive_init,constitutive_state_old,constitutive_state_new,material_Cslip_66 + implicit none +! + integer(pInt) CPFEM_en, CPFEM_in, cp_en, CPFEM_ngens, i,j,k,l,m,n, e + real(pReal), dimension (3,3) :: ffn,ffn1,Kirchhoff_bar + real(pReal), dimension (3,3,3,3) :: H_bar + real(pReal), dimension(CPFEM_ngens) :: CPFEM_stress + real(pReal), dimension(CPFEM_ngens,CPFEM_ngens) :: CPFEM_jaco + real(pReal) Temperature,CPFEM_dt,J_inverse + integer(pInt) CPFEM_mode ! 1: regular computation with aged results& + ! 2: regular computation& + ! 3: collection of FEM data& + ! 4: recycling of former results (MARC speciality)& + ! 5: record tangent from former converged inc& + ! 6: restore tangent from former converged inc + logical CPFEM_updateJaco +! + if (.not. CPFEM_init_done) then ! initialization step (three dimensional stress state check missing?) + call math_init() + call mesh_init() + call lattice_init() + call constitutive_init() + call CPFEM_init(Temperature) + CPFEM_init_done = .true. + endif +! + cp_en = mesh_FEasCP('elem',CPFEM_en) + if (cp_en == 1 .and. CPFEM_in == 1) & + write(6,'(a6,x,i4,x,a4,x,i4,x,a10,x,f8.4,x,a10,x,i2,x,a10,x,i2,x,a10,x,i2,x,a10,x,i2)') & + 'elem',cp_en,'IP',CPFEM_in,& + 'theTime',theTime,'theInc',theInc,'theCycle',theCycle,'theLovl',theLovl,& + 'mode',CPFEM_mode +! + select case (CPFEM_mode) + case (2,1) ! regular computation (with aging of results) + if (.not. CPFEM_calc_done) then ! puuh, me needs doing all the work... + write (6,*) 'puuh me needs doing all the work', cp_en + if (CPFEM_mode == 1) then ! age results at start of new increment + CPFEM_Fp_old = CPFEM_Fp_new + constitutive_state_old = constitutive_state_new + GIA_rVect_old = GIA_rVect_new + write (6,*) '#### aged results' + endif + debug_cutbackDistribution = 0_pInt ! initialize debugging data + debug_InnerLoopDistribution = 0_pInt + debug_OuterLoopDistribution = 0_pInt +! + do e=1,mesh_NcpElems ! ## this shall be done in a parallel loop in the future ## + do i=1,FE_Nips(mesh_element(2,e)) ! iterate over all IPs of this element's type + debugger = (e==1 .and. i==1) ! switch on debugging for first IP in first element + call CPFEM_MaterialPoint(CPFEM_updateJaco, CPFEM_dt, i, e) + enddo + enddo + call debug_info() ! output of debugging/performance statistics + CPFEM_calc_done = .true. ! now calc is done + endif +! translate from P and dP/dF to CS and dCS/dE + Kirchhoff_bar = matmul(CPFEM_PK1_bar(:,:,CPFEM_in, cp_en),transpose(CPFEM_ffn1_bar(:,:,CPFEM_in, cp_en))) + J_inverse = 1.0_pReal/math_det3x3(CPFEM_ffn1_bar(:,:,CPFEM_in, cp_en)) + CPFEM_stress_bar(1:CPFEM_ngens,CPFEM_in,cp_en) = math_Mandel33to6(J_inverse*Kirchhoff_bar) +! + H_bar = 0.0_pReal + forall(i=1:3,j=1:3,k=1:3,l=1:3,m=1:3,n=1:3) & + H_bar(i,j,k,l) = H_bar(i,j,k,l) + & + (CPFEM_ffn1_bar(j,m,CPFEM_in,cp_en)*CPFEM_ffn1_bar(l,n,CPFEM_in,cp_en)*CPFEM_dPdF_bar(i,m,k,n,CPFEM_in,cp_en) - & + math_I3(j,l)*CPFEM_ffn1_bar(i,m,CPFEM_in,cp_en)*CPFEM_PK1_bar(k,m,CPFEM_in,cp_en)) + & + 0.5_pReal*(math_I3(i,k)*Kirchhoff_bar(j,l) + math_I3(j,l)*Kirchhoff_bar(i,k) + & + math_I3(i,l)*Kirchhoff_bar(j,k) + math_I3(j,k)*Kirchhoff_bar(i,l)) + CPFEM_jaco_bar(1:CPFEM_ngens,1:CPFEM_ngens,CPFEM_in,cp_en) = math_Mandel3333to66(J_inverse*H_bar) +! + case (3) ! collect and return odd result + CPFEM_Temperature(CPFEM_in,cp_en) = Temperature + CPFEM_ffn_bar(:,:,CPFEM_in,cp_en) = ffn + CPFEM_ffn1_bar(:,:,CPFEM_in,cp_en) = ffn1 + CPFEM_stress_bar(1:CPFEM_ngens,CPFEM_in,cp_en) = CPFEM_odd_stress + CPFEM_jaco_bar(1:CPFEM_ngens,1:CPFEM_ngens,CPFEM_in,cp_en) = CPFEM_odd_jacobian*math_identity2nd(CPFEM_ngens) + CPFEM_calc_done = .false. + + case (4) ! do nothing since we can recycle the former results (MARC specialty) + case (5) ! record consistent tangent at beginning of new increment + CPFEM_jaco_knownGood = CPFEM_jaco_bar + case (6) ! restore consistent tangent after cutback + CPFEM_jaco_bar = CPFEM_jaco_knownGood + end select +! +! return the local stress and the jacobian from storage + CPFEM_stress(1:CPFEM_ngens) = CPFEM_stress_bar(1:CPFEM_ngens,CPFEM_in,cp_en) + CPFEM_jaco(1:CPFEM_ngens,1:CPFEM_ngens) = CPFEM_jaco_bar(1:CPFEM_ngens,1:CPFEM_ngens,CPFEM_in,cp_en) + if (cp_en == 1 .and. CPFEM_in == 1) write (6,*) 'stress',CPFEM_stress + if (cp_en == 1 .and. CPFEM_in == 1 .and. CPFEM_updateJaco) write (6,*) 'stiffness',CPFEM_jaco +! + return +! + END SUBROUTINE +! +!********************************************************** +!*** calculate the material point behaviour *** +!********************************************************** + SUBROUTINE CPFEM_MaterialPoint(& + updateJaco,& ! flag to initiate Jacobian updating + CPFEM_dt,& ! Time increment (dt) + CPFEM_in,& ! Integration point number + cp_en) ! Element number +! + use prec + use debug + use math, only: math_pDecomposition,math_RtoEuler,inDeg,math_I3,math_invert3x3,math_permut,math_invert,math_delta + use IO, only: IO_error + use mesh, only: mesh_element + use crystallite + use constitutive + implicit none +! + character(len=128) msg + integer(pInt) cp_en,CPFEM_in,grain,max_cutbacks,i,j,k,l,m,n,iBoun,NRiter,dummy,ii,jj,kk,ll + logical updateJaco,error,NRconvergent,failed + real(pReal) CPFEM_dt,volfrac,dTime,shMod,C_kb,resNorm,resMax,subStep,subFrac,temp1,temp2 + real(pReal), dimension(3,3) :: F0_bar,F1_bar,dF_bar + real(pReal), dimension(3,3) :: U,R + real(pReal), dimension(3,3,8) :: PK1,Fp0,Fp1,Fe1,F1,F0,dFgrain,dFg_cor + real(pReal), dimension(3,3,12) :: GPK1,GF1,Nye,GRB1 + real(pReal), dimension(3,3,3,3,8) :: dPdF + real(pReal), dimension(3,3,3,3,12) :: dRdX1 + real(pReal), dimension(36) :: var,res + real(pReal), dimension(36,36) :: dresdvar,dvardres + real(pReal), dimension(3,12) :: rx,rVect + real(pReal), dimension(12) :: NyeNorm + real(pReal), dimension(constitutive_maxNstatevars,8) :: state0,state1 +! + if (texture_Ngrains(mesh_element(4,cp_en)) /= 8_pInt) then + call IO_error(800) + return + endif +! + CPFEM_PK1_bar(:,:,CPFEM_in,cp_en) = 0.0_pReal ! zero out average first PK stress + if (updateJaco) CPFEM_dPdF_bar(:,:,:,:,CPFEM_in,cp_en) = 0.0_pReal ! zero out average consistent tangent +! +! ------------- GIA loop -------------------- +! +! collect information + shMod = 0.2_pReal*(material_C11(1) - material_C12(1)) + 0.3_pReal*material_C44(1) ! equivalent shear modulus + C_kb = material_bg(1)*shMod/material_GrainSize(1) ! equivalent boundary stiffness +! + F0_bar = CPFEM_ffn_bar(:,:,CPFEM_in,cp_en) ! effective deformation gradient at t_n + state0 = constitutive_state_old(:,:,CPFEM_in,cp_en) ! state variables at t_n + Fp0 = CPFEM_Fp_old(:,:,:,CPFEM_in,cp_en) ! grain plastic def. gradient at t_n + rVect = GIA_rVect_old(:,:,CPFEM_in,cp_en) ! relaxation vectors from previous convergent step +! + dF_bar = CPFEM_ffn1_bar(:,:,CPFEM_in,cp_en) - CPFEM_ffn_bar(:,:,CPFEM_in,cp_en) ! deformation gradient increment + subFrac = 0.0_pReal + subStep = 1.0_pReal +! +! Substepping procedure to improve N-R iteration +SubStepping: do + dTime = subStep*CPFEM_dt + call GIA_RelaxedDeformation(F0,F0_bar,rVect) ! def. gradient of indiv. grains at t_n + F1_bar = F0_bar + subStep*dF_bar ! effective def. gradient at t_n+1 + forall (iBoun=1:12,i=1:3) var(3_pInt*(iBoun-1_pInt)+i) = rVect(i,iBoun) ! primary variable: relaxation vector +! +! Newton-Raphson iteration block + NRiter = 1_pInt +NRIteration: do + forall (iBoun=1:12,i=1:3) rx(i,iBoun) = var(3_pInt*(iBoun-1_pInt)+i) ! relaxation vectors (guess) +! +! deformation gradients of grains at t_n+1 (guess) + call GIA_RelaxedDeformation(F1,F1_bar,rx) +! +! -------------- grain loop ----------------- + do grain = 1,8 + call SingleCrystallite(msg,PK1(:,:,grain),dPdF(:,:,:,:,grain),& + CPFEM_results(5:4+constitutive_Nresults(grain,CPFEM_in,cp_en),grain,CPFEM_in,cp_en),& + Fp1(:,:,grain),Fe1(:,:,grain),state1(:,grain),& ! output up to here + dTime,cp_en,CPFEM_in,grain,.true.,& + CPFEM_Temperature(CPFEM_in,cp_en),F1(:,:,grain),F0(:,:,grain),Fp0(:,:,grain),state0(:,grain)) + if (msg /= 'ok') then ! solution not reached --> exit NRIteration + write(6,*) 'GIA: grain loop failed to converge within allowable step-size' + NRconvergent = .false. + exit NRiteration + endif + enddo ! grain loop +! +! calculate the deformation jump and stress jump across the boundaries + dFgrain = F1 - F0 + call GIA_BoundaryJump(GF1,F1) + call GIA_BoundaryJump(GPK1,PK1) +! +! compute the Nye tensor at the boundary + Nye = 0.0_pReal + NyeNorm = 0.0_pReal + do iBoun = 1,12 + do i = 1,3 + do j = 1,3 + do k = 1,3 + do l = 1,3 + Nye(i,j,iBoun) = Nye(i,j,iBoun) - 0.5_pReal*math_permut(j,k,l)*GIA_bNorm(k,iBoun)*GF1(i,l,iBoun) + enddo + enddo + NyeNorm(iBoun) = NyeNorm(iBoun) + Nye(i,j,iBoun)*Nye(i,j,iBoun) + enddo + enddo + NyeNorm(iBoun) = sqrt(NyeNorm(iBoun)) + if (NyeNorm(iBoun) > 1.0e-8_pReal) Nye(:,:,iBoun) = Nye(:,:,iBoun)/NyeNorm(iBoun) + enddo +! +! compute the stress-like penalty at the boundary + GRB1 = 0.0_pReal + do iBoun = 1,12 + do i = 1,3 + do j = 1,3 + do k = 1,3 + do l = 1,3 + GRB1(i,j,iBoun) = GRB1(i,j,iBoun) + Nye(i,k,iBoun)*GIA_bNorm(l,iBoun)*math_permut(k,l,j) + enddo + enddo + enddo + enddo + GRB1(:,:,iBoun) = 0.5_pReal*(C_kb + C_kb)*GRB1(:,:,iBoun) + enddo +! +! compute the resiudal of stress at the boundary + res = 0.0_pReal + resNorm = 0.0_pReal + do iBoun = 1,12 + do j = 1,3 + do i = 1,3 + res(3_pInt*(iBoun-1_pInt)+j) = res(3_pInt*(iBoun-1_pInt)+j) - & + GIA_bNorm(i,iBoun)*(GPK1(i,j,iBoun) - GRB1(i,j,iBoun)) + enddo + resNorm = resNorm + res(3_pInt*(iBoun-1_pInt)+j)*res(3_pInt*(iBoun-1_pInt)+j) + enddo + enddo + resNorm = sqrt(resNorm) +! +! write(6,'(x,a,i3,a,i3,a,i3,a,e10.4)')'EL = ',cp_en,':IP = ',CPFEM_in,':Iter = ',NRiter,':RNorm = ',resNorm + if (NRiter == 1_pInt) resMax = resNorm + if ((resNorm < resToler*resMax) .or. (resNorm < resAbsol)) then ! resNorm < tolerance ===> convergent + NRconvergent = .true. + exit NRiteration + elseif ((NRiter > NRiterMax) .or. (resNorm > resBound*resMax)) then ! resNorm > up. bound ===> substepping + NRconvergent = .false. + exit NRiteration + else ! update the residual + dRdX1 = 0.0_pReal + do iBoun = 1,12 + if (NyeNorm(iBoun) < 1.0e-8_pReal) NyeNorm(iBoun) = 1.0e-8_pReal + do i = 1,3 + do j = 1,3 + do k = 1,3 + do l = 1,3 + temp1 = 0.0_pReal + temp2 = 0.0_pReal + do ii = 1,3 + do jj = 1,3 + do kk = 1,3 + temp1 = temp1 + GIA_bNorm(jj,iBoun)*math_permut(ii,jj,j)*math_delta(i,k)* & + GIA_bNorm(kk,iBoun)*math_permut(ii,kk,l) + do ll = 1,3 + temp2 = temp2 + Nye(i,ii,iBoun)*GIA_bNorm(jj,iBoun)*math_permut(ii,jj,j)* & + Nye(k,kk,iBoun)*GIA_bNorm(ll,iBoun)*math_permut(kk,ll,l) + enddo + enddo + enddo + enddo + dRdX1(i,j,k,l,iBoun) = 0.25_pReal*(C_kb + C_kb)*(temp1 - temp2)/NyeNorm(iBoun) + enddo + enddo + enddo + enddo + enddo + call GIA_JacobianMatrix(dresdvar,dPdF,dRdX1) + dvardres = 0.0_pReal + call math_invert(36,dresdvar,dvardres,dummy,failed) + if (failed) then + write(6,*) 'GIA: failed to invert the Jacobian' + NRconvergent = .false. + exit NRiteration + endif + forall (i=1:36,j=1:36) var(i) = var(i) - dvardres(i,j)*res(j) + endif +! + NRiter = NRiter + 1_pInt + enddo NRIteration ! End of N-R iteration blok +! + if (.not. NRconvergent) then + subStep = 0.5_pReal*subStep + else + subFrac = subFrac + subStep + subStep = 1.0_pReal - subFrac + Fp0 = Fp1 + F0_bar = F1_bar + state0 = state1 + rVect = rx + endif +! + if (subStep < subStepMin) exit SubStepping + enddo SubStepping ! End of substepping blok +! +! ------------- GIA loop (end) -------------- +! +! return to the general subroutine when convergence is not reached + if (.not. NRconvergent) then + write(6,'(x,a)') 'GIA: convergence is not reached within allowable step-size' + write(6,'(x,a,i3,a,i3)') 'Element: ',cp_en,' @ IP: ',CPFEM_in + call IO_error(600) + return + endif +! +! updates all variables, deformation gradients, and vectors + GIA_rVect_new(:,:,CPFEM_in,cp_en) = rVect + CPFEM_Fp_new(:,:,:,CPFEM_in,cp_en) = Fp1 + constitutive_state_new(:,:,CPFEM_in,cp_en) = state1 +! +! compute the effective stress and consistent tangent + call GIA_TangentCorrection(dFgrain,dFg_cor) + do grain = 1,texture_Ngrains(mesh_element(4,cp_en)) + volfrac = constitutive_matVolFrac(grain,CPFEM_in,cp_en)*constitutive_texVolFrac(grain,CPFEM_in,cp_en) + CPFEM_PK1_bar(:,:,CPFEM_in,cp_en) = CPFEM_PK1_bar(:,:,CPFEM_in,cp_en) + & + volfrac*PK1(:,:,grain) ! average Cauchy stress + if (updateJaco) then ! consistent tangent + do i = 1,3 + do j = 1,3 + CPFEM_dPdF_bar(:,:,i,j,CPFEM_in,cp_en) = CPFEM_dPdF_bar(:,:,i,j,CPFEM_in,cp_en) + & + volfrac*dPdF(:,:,i,j,grain)*dFg_cor(i,j,grain) + enddo + enddo + endif +! +! update results plotted in MENTAT + call math_pDecomposition(Fe1(:,:,grain),U,R,error) ! polar decomposition + if (error) then + write(6,*) Fe1(:,:,grain) + write(6,*) 'polar decomposition' + write(6,*) 'Grain: ',grain + write(6,*) 'Integration point: ',CPFEM_in + write(6,*) 'Element: ',mesh_element(1,cp_en) + call IO_error(650) + return + endif + CPFEM_results(1:3,grain,CPFEM_in,cp_en) = math_RtoEuler(transpose(R))*inDeg ! orientation + CPFEM_results(4 ,grain,CPFEM_in,cp_en) = volfrac ! volume fraction of orientation + enddo +! + return +! + END SUBROUTINE +! +! +!******************************************************************** +! Calculates the relaxed deformation gradients of grains +!******************************************************************** + subroutine GIA_RelaxedDeformation(& + F,& ! relaxed deformation gradient of grains + F_bar,& ! effective deformation gradient + r) ! relaxation vectors at boundary +! + implicit none +! + real(pReal), dimension(3,3) :: F_bar + real(pReal), dimension(3,3,8) :: F + real(pReal), dimension(3,12) :: r,n + integer(pInt) i,j,iBoun,grain +! + n = GIA_bNorm + do i = 1,3 + do j = 1,3 + F(i,j,1) = F_bar(i,j) + n(i, 1)*r(j, 1) + n(i, 5)*r(j, 5) + n(i, 9)*r(j, 9) + F(i,j,2) = F_bar(i,j) - n(i, 1)*r(j, 1) + n(i, 6)*r(j, 6) + n(i,10)*r(j,10) + F(i,j,3) = F_bar(i,j) + n(i, 2)*r(j, 2) - n(i, 5)*r(j, 5) + n(i,11)*r(j,11) + F(i,j,4) = F_bar(i,j) - n(i, 2)*r(j, 2) - n(i, 6)*r(j, 6) + n(i,12)*r(j,12) + F(i,j,5) = F_bar(i,j) + n(i, 3)*r(j, 3) + n(i, 7)*r(j, 7) - n(i, 9)*r(j, 9) + F(i,j,6) = F_bar(i,j) - n(i, 3)*r(j, 3) + n(i, 8)*r(j, 8) - n(i,10)*r(j,10) + F(i,j,7) = F_bar(i,j) + n(i, 4)*r(j, 4) - n(i, 7)*r(j, 7) - n(i,11)*r(j,11) + F(i,j,8) = F_bar(i,j) - n(i, 4)*r(j, 4) - n(i, 8)*r(j, 8) - n(i,12)*r(j,12) + enddo + enddo +! + return +! + END SUBROUTINE +! +! +!******************************************************************** +! Calculates the jump of tensors across the grain boundary +!******************************************************************** + subroutine GIA_BoundaryJump(& + F_boun,& ! tensor jump across the boundary + F_bulk) ! bulk tensor +! + implicit none +! + real(pReal), dimension(3,3,12) :: F_boun + real(pReal), dimension(3,3,8) :: F_bulk + integer(pInt) i,j,iBoun,grain +! + F_boun(:,:, 1) = F_bulk(:,:,2) - F_bulk(:,:,1) + F_boun(:,:, 2) = F_bulk(:,:,4) - F_bulk(:,:,3) + F_boun(:,:, 3) = F_bulk(:,:,6) - F_bulk(:,:,5) + F_boun(:,:, 4) = F_bulk(:,:,8) - F_bulk(:,:,7) + F_boun(:,:, 5) = F_bulk(:,:,3) - F_bulk(:,:,1) + F_boun(:,:, 6) = F_bulk(:,:,4) - F_bulk(:,:,2) + F_boun(:,:, 7) = F_bulk(:,:,7) - F_bulk(:,:,5) + F_boun(:,:, 8) = F_bulk(:,:,8) - F_bulk(:,:,6) + F_boun(:,:, 9) = F_bulk(:,:,5) - F_bulk(:,:,1) + F_boun(:,:,10) = F_bulk(:,:,6) - F_bulk(:,:,2) + F_boun(:,:,11) = F_bulk(:,:,7) - F_bulk(:,:,3) + F_boun(:,:,12) = F_bulk(:,:,8) - F_bulk(:,:,4) +! + return +! + END SUBROUTINE +! +! +!******************************************************************** +! Calculates the jump of tensors across the grain boundary +!******************************************************************** + subroutine GIA_JacobianMatrix(& + dresdvar,& ! Jacobian matrix + dPdF,& ! stress consistent tangent of bulk + dRdX) ! stress-like penalty tangent at boundary +! + implicit none +! + real(pReal), dimension(3,3,3,3,8) :: dPdF + real(pReal), dimension(3,3,3,3,12) :: dRdX + real(pReal), dimension(36,36) :: dresdvar + real(pReal), dimension(3,12) :: n + integer(pInt) i,j,k,l +! + n = GIA_bNorm + dresdvar = 0.0_pReal + do i = 1,3 + do k = 1,3 + do l = 1,3 + do j = 1,3 +! +! at boundary 1, influenced by boundary +5, -6, +9, -10 + dresdvar(( 1-1)*3 + j,( 1-1)*3 + l) = dresdvar(( 1-1)*3 + j,( 1-1)*3 + l) & + + (dPdF(i,j,k,l, 1) + dPdF(i,j,k,l, 2))*n(i, 1)*n(k, 1) & + + (dRdX(i,j,k,l, 1) + dRdX(i,j,k,l, 1))*n(i, 1)*n(k, 1) + dresdvar(( 1-1)*3 + j,( 5-1)*3 + l) = dresdvar(( 1-1)*3 + j,( 5-1)*3 + l) + dPdF(i,j,k,l, 1)*n(i, 1)*n(k, 5) & + + dRdX(i,j,k,l, 1)*n(i, 1)*n(k, 5) + dresdvar(( 1-1)*3 + j,( 6-1)*3 + l) = dresdvar(( 1-1)*3 + j,( 6-1)*3 + l) - dPdF(i,j,k,l, 2)*n(i, 1)*n(k, 6) & + - dRdX(i,j,k,l, 1)*n(i, 1)*n(k, 6) + dresdvar(( 1-1)*3 + j,( 9-1)*3 + l) = dresdvar(( 1-1)*3 + j,( 9-1)*3 + l) + dPdF(i,j,k,l, 1)*n(i, 1)*n(k, 9) & + + dRdX(i,j,k,l, 1)*n(i, 1)*n(k, 9) + dresdvar(( 1-1)*3 + j,(10-1)*3 + l) = dresdvar(( 1-1)*3 + j,(10-1)*3 + l) - dPdF(i,j,k,l, 2)*n(i, 1)*n(k,10) & + - dRdX(i,j,k,l, 1)*n(i, 1)*n(k,10) +! +! at boundary 2, influenced by boundary -5, +6, +11, -12 + dresdvar(( 2-1)*3 + j,( 2-1)*3 + l) = dresdvar(( 2-1)*3 + j,( 2-1)*3 + l) & + + (dPdF(i,j,k,l, 3) + dPdF(i,j,k,l, 4))*n(i, 2)*n(k, 2) & + + (dRdX(i,j,k,l, 2) + dRdX(i,j,k,l, 2))*n(i, 2)*n(k, 2) + dresdvar(( 2-1)*3 + j,( 5-1)*3 + l) = dresdvar(( 2-1)*3 + j,( 5-1)*3 + l) - dPdF(i,j,k,l, 3)*n(i, 2)*n(k, 5) & + - dRdX(i,j,k,l, 2)*n(i, 2)*n(k, 5) + dresdvar(( 2-1)*3 + j,( 6-1)*3 + l) = dresdvar(( 2-1)*3 + j,( 6-1)*3 + l) + dPdF(i,j,k,l, 4)*n(i, 2)*n(k, 6) & + + dRdX(i,j,k,l, 2)*n(i, 2)*n(k, 6) + dresdvar(( 2-1)*3 + j,(11-1)*3 + l) = dresdvar(( 2-1)*3 + j,(11-1)*3 + l) + dPdF(i,j,k,l, 3)*n(i, 2)*n(k,11) & + + dRdX(i,j,k,l, 2)*n(i, 2)*n(k,11) + dresdvar(( 2-1)*3 + j,(12-1)*3 + l) = dresdvar(( 2-1)*3 + j,(12-1)*3 + l) - dPdF(i,j,k,l, 4)*n(i, 2)*n(k,12) & + - dRdX(i,j,k,l, 2)*n(i, 2)*n(k,12) +! +! at boundary 3, influenced by boundary +7, -8, -9, +10 + dresdvar(( 3-1)*3 + j,( 3-1)*3 + l) = dresdvar(( 3-1)*3 + j,( 3-1)*3 + l) & + + (dPdF(i,j,k,l, 5) + dPdF(i,j,k,l, 6))*n(i, 3)*n(k, 3) & + + (dRdX(i,j,k,l, 3) + dRdX(i,j,k,l, 3))*n(i, 3)*n(k, 3) + dresdvar(( 3-1)*3 + j,( 7-1)*3 + l) = dresdvar(( 3-1)*3 + j,( 7-1)*3 + l) + dPdF(i,j,k,l, 5)*n(i, 3)*n(k, 7) & + + dRdX(i,j,k,l, 3)*n(i, 3)*n(k, 7) + dresdvar(( 3-1)*3 + j,( 8-1)*3 + l) = dresdvar(( 3-1)*3 + j,( 8-1)*3 + l) - dPdF(i,j,k,l, 6)*n(i, 3)*n(k, 8) & + - dRdX(i,j,k,l, 3)*n(i, 3)*n(k, 8) + dresdvar(( 3-1)*3 + j,( 9-1)*3 + l) = dresdvar(( 3-1)*3 + j,( 9-1)*3 + l) - dPdF(i,j,k,l, 5)*n(i, 3)*n(k, 9) & + - dRdX(i,j,k,l, 3)*n(i, 3)*n(k, 9) + dresdvar(( 3-1)*3 + j,(10-1)*3 + l) = dresdvar(( 3-1)*3 + j,(10-1)*3 + l) + dPdF(i,j,k,l, 6)*n(i, 3)*n(k,10) & + + dRdX(i,j,k,l, 3)*n(i, 3)*n(k,10) +! +! at boundary 4, influenced by boundary -7, +8, -11, +12 + dresdvar(( 4-1)*3 + j,( 4-1)*3 + l) = dresdvar(( 4-1)*3 + j,( 4-1)*3 + l) & + + (dPdF(i,j,k,l, 7) + dPdF(i,j,k,l, 8))*n(i, 4)*n(k, 4) & + + (dRdX(i,j,k,l, 4) + dRdX(i,j,k,l, 4))*n(i, 4)*n(k, 4) + dresdvar(( 4-1)*3 + j,( 7-1)*3 + l) = dresdvar(( 4-1)*3 + j,( 7-1)*3 + l) - dPdF(i,j,k,l, 7)*n(i, 4)*n(k, 7) & + - dRdX(i,j,k,l, 4)*n(i, 4)*n(k, 7) + dresdvar(( 4-1)*3 + j,( 8-1)*3 + l) = dresdvar(( 4-1)*3 + j,( 8-1)*3 + l) + dPdF(i,j,k,l, 8)*n(i, 4)*n(k, 8) & + + dRdX(i,j,k,l, 4)*n(i, 4)*n(k, 8) + dresdvar(( 4-1)*3 + j,(11-1)*3 + l) = dresdvar(( 4-1)*3 + j,(11-1)*3 + l) - dPdF(i,j,k,l, 7)*n(i, 4)*n(k,11) & + - dRdX(i,j,k,l, 4)*n(i, 4)*n(k,11) + dresdvar(( 4-1)*3 + j,(12-1)*3 + l) = dresdvar(( 4-1)*3 + j,(12-1)*3 + l) + dPdF(i,j,k,l, 8)*n(i, 4)*n(k,12) & + + dRdX(i,j,k,l, 4)*n(i, 4)*n(k,12) +! +! at boundary 5, influenced by boundary +1, -2, +9, -11 + dresdvar(( 5-1)*3 + j,( 5-1)*3 + l) = dresdvar(( 5-1)*3 + j,( 5-1)*3 + l) & + + (dPdF(i,j,k,l, 1) + dPdF(i,j,k,l, 3))*n(i, 5)*n(k, 5) & + + (dRdX(i,j,k,l, 5) + dRdX(i,j,k,l, 5))*n(i, 5)*n(k, 5) + dresdvar(( 5-1)*3 + j,( 1-1)*3 + l) = dresdvar(( 5-1)*3 + j,( 1-1)*3 + l) + dPdF(i,j,k,l, 1)*n(i, 5)*n(k, 1) & + + dRdX(i,j,k,l, 5)*n(i, 5)*n(k, 1) + dresdvar(( 5-1)*3 + j,( 2-1)*3 + l) = dresdvar(( 5-1)*3 + j,( 2-1)*3 + l) - dPdF(i,j,k,l, 3)*n(i, 5)*n(k, 2) & + - dRdX(i,j,k,l, 5)*n(i, 5)*n(k, 2) + dresdvar(( 5-1)*3 + j,( 9-1)*3 + l) = dresdvar(( 5-1)*3 + j,( 9-1)*3 + l) + dPdF(i,j,k,l, 1)*n(i, 5)*n(k, 9) & + + dRdX(i,j,k,l, 5)*n(i, 5)*n(k, 9) + dresdvar(( 5-1)*3 + j,(11-1)*3 + l) = dresdvar(( 5-1)*3 + j,(11-1)*3 + l) - dPdF(i,j,k,l, 3)*n(i, 5)*n(k,11) & + - dRdX(i,j,k,l, 5)*n(i, 5)*n(k,11) +! +! at boundary 6, influenced by boundary -1, +2, +10, -12 + dresdvar(( 6-1)*3 + j,( 6-1)*3 + l) = dresdvar(( 6-1)*3 + j,( 6-1)*3 + l) & + + (dPdF(i,j,k,l, 2) + dPdF(i,j,k,l, 4))*n(i, 6)*n(k, 6) & + + (dRdX(i,j,k,l, 6) + dRdX(i,j,k,l, 6))*n(i, 6)*n(k, 6) + dresdvar(( 6-1)*3 + j,( 1-1)*3 + l) = dresdvar(( 6-1)*3 + j,( 1-1)*3 + l) - dPdF(i,j,k,l, 2)*n(i, 6)*n(k, 1) & + - dRdX(i,j,k,l, 6)*n(i, 6)*n(k, 1) + dresdvar(( 6-1)*3 + j,( 2-1)*3 + l) = dresdvar(( 6-1)*3 + j,( 2-1)*3 + l) + dPdF(i,j,k,l, 4)*n(i, 6)*n(k, 2) & + + dRdX(i,j,k,l, 6)*n(i, 6)*n(k, 2) + dresdvar(( 6-1)*3 + j,(10-1)*3 + l) = dresdvar(( 6-1)*3 + j,(10-1)*3 + l) + dPdF(i,j,k,l, 2)*n(i, 6)*n(k,10) & + + dRdX(i,j,k,l, 6)*n(i, 6)*n(k,10) + dresdvar(( 6-1)*3 + j,(12-1)*3 + l) = dresdvar(( 6-1)*3 + j,(12-1)*3 + l) - dPdF(i,j,k,l, 4)*n(i, 6)*n(k,12) & + - dRdX(i,j,k,l, 6)*n(i, 6)*n(k,12) +! +! at boundary 7, influenced by boundary +3, -4, -9, +11 + dresdvar(( 7-1)*3 + j,( 7-1)*3 + l) = dresdvar(( 7-1)*3 + j,( 7-1)*3 + l) & + + (dPdF(i,j,k,l, 5) + dPdF(i,j,k,l, 7))*n(i, 7)*n(k, 7) & + + (dRdX(i,j,k,l, 7) + dRdX(i,j,k,l, 7))*n(i, 7)*n(k, 7) + dresdvar(( 7-1)*3 + j,( 3-1)*3 + l) = dresdvar(( 7-1)*3 + j,( 3-1)*3 + l) + dPdF(i,j,k,l, 5)*n(i, 7)*n(k, 3) & + + dRdX(i,j,k,l, 7)*n(i, 7)*n(k, 3) + dresdvar(( 7-1)*3 + j,( 4-1)*3 + l) = dresdvar(( 7-1)*3 + j,( 4-1)*3 + l) - dPdF(i,j,k,l, 7)*n(i, 7)*n(k, 4) & + - dRdX(i,j,k,l, 7)*n(i, 7)*n(k, 4) + dresdvar(( 7-1)*3 + j,( 9-1)*3 + l) = dresdvar(( 7-1)*3 + j,( 9-1)*3 + l) - dPdF(i,j,k,l, 5)*n(i, 7)*n(k, 9) & + - dRdX(i,j,k,l, 7)*n(i, 7)*n(k, 9) + dresdvar(( 7-1)*3 + j,(11-1)*3 + l) = dresdvar(( 7-1)*3 + j,(11-1)*3 + l) + dPdF(i,j,k,l, 7)*n(i, 7)*n(k,11) & + + dRdX(i,j,k,l, 7)*n(i, 7)*n(k,11) +! +! at boundary 8, influenced by boundary -3, +4, -10, +12 + dresdvar(( 8-1)*3 + j,( 8-1)*3 + l) = dresdvar(( 8-1)*3 + j,( 8-1)*3 + l) & + + (dPdF(i,j,k,l, 6) + dPdF(i,j,k,l, 8))*n(i, 8)*n(k, 8) & + + (dRdX(i,j,k,l, 8) + dRdX(i,j,k,l, 8))*n(i, 8)*n(k, 8) + dresdvar(( 8-1)*3 + j,( 3-1)*3 + l) = dresdvar(( 8-1)*3 + j,( 3-1)*3 + l) - dPdF(i,j,k,l, 6)*n(i, 8)*n(k, 3) & + - dRdX(i,j,k,l, 8)*n(i, 8)*n(k, 3) + dresdvar(( 8-1)*3 + j,( 4-1)*3 + l) = dresdvar(( 8-1)*3 + j,( 4-1)*3 + l) + dPdF(i,j,k,l, 8)*n(i, 8)*n(k, 4) & + + dRdX(i,j,k,l, 8)*n(i, 8)*n(k, 4) + dresdvar(( 8-1)*3 + j,(10-1)*3 + l) = dresdvar(( 8-1)*3 + j,(10-1)*3 + l) - dPdF(i,j,k,l, 6)*n(i, 8)*n(k,10) & + - dRdX(i,j,k,l, 8)*n(i, 8)*n(k,10) + dresdvar(( 8-1)*3 + j,(12-1)*3 + l) = dresdvar(( 8-1)*3 + j,(12-1)*3 + l) + dPdF(i,j,k,l, 8)*n(i, 8)*n(k,12) & + + dRdX(i,j,k,l, 8)*n(i, 8)*n(k,12) +! +! at boundary 9, influenced by boundary +1, -3, +5, -7 + dresdvar(( 9-1)*3 + j,( 9-1)*3 + l) = dresdvar(( 9-1)*3 + j,( 9-1)*3 + l) & + + (dPdF(i,j,k,l, 1) + dPdF(i,j,k,l, 5))*n(i, 9)*n(k, 9) & + + (dRdX(i,j,k,l, 9) + dRdX(i,j,k,l, 9))*n(i, 9)*n(k, 9) + dresdvar(( 9-1)*3 + j,( 1-1)*3 + l) = dresdvar(( 9-1)*3 + j,( 1-1)*3 + l) + dPdF(i,j,k,l, 1)*n(i, 9)*n(k, 1) & + + dRdX(i,j,k,l, 9)*n(i, 9)*n(k, 1) + dresdvar(( 9-1)*3 + j,( 3-1)*3 + l) = dresdvar(( 9-1)*3 + j,( 3-1)*3 + l) - dPdF(i,j,k,l, 5)*n(i, 9)*n(k, 3) & + - dRdX(i,j,k,l, 9)*n(i, 9)*n(k, 3) + dresdvar(( 9-1)*3 + j,( 5-1)*3 + l) = dresdvar(( 9-1)*3 + j,( 5-1)*3 + l) + dPdF(i,j,k,l, 1)*n(i, 9)*n(k, 5) & + + dRdX(i,j,k,l, 9)*n(i, 9)*n(k, 5) + dresdvar(( 9-1)*3 + j,( 7-1)*3 + l) = dresdvar(( 9-1)*3 + j,( 7-1)*3 + l) - dPdF(i,j,k,l, 5)*n(i, 9)*n(k, 7) & + - dRdX(i,j,k,l, 9)*n(i, 9)*n(k, 7) +! +! at boundary 10, influenced by boundary -1, +3, +6, -8 + dresdvar((10-1)*3 + j,(10-1)*3 + l) = dresdvar((10-1)*3 + j,(10-1)*3 + l) & + + (dPdF(i,j,k,l, 2) + dPdF(i,j,k,l, 6))*n(i,10)*n(k,10) & + + (dRdX(i,j,k,l,10) + dRdX(i,j,k,l,10))*n(i,10)*n(k,10) + dresdvar((10-1)*3 + j,( 1-1)*3 + l) = dresdvar((10-1)*3 + j,( 1-1)*3 + l) - dPdF(i,j,k,l, 2)*n(i,10)*n(k, 1) & + - dRdX(i,j,k,l,10)*n(i,10)*n(k, 1) + dresdvar((10-1)*3 + j,( 3-1)*3 + l) = dresdvar((10-1)*3 + j,( 3-1)*3 + l) + dPdF(i,j,k,l, 6)*n(i,10)*n(k, 3) & + + dRdX(i,j,k,l,10)*n(i,10)*n(k, 3) + dresdvar((10-1)*3 + j,( 6-1)*3 + l) = dresdvar((10-1)*3 + j,( 6-1)*3 + l) + dPdF(i,j,k,l, 2)*n(i,10)*n(k, 6) & + + dRdX(i,j,k,l,10)*n(i,10)*n(k, 6) + dresdvar((10-1)*3 + j,( 8-1)*3 + l) = dresdvar((10-1)*3 + j,( 8-1)*3 + l) - dPdF(i,j,k,l, 6)*n(i,10)*n(k, 8) & + - dRdX(i,j,k,l,10)*n(i,10)*n(k, 8) +! +! at boundary 11, influenced by boundary +2, -4, -5, +7 + dresdvar((11-1)*3 + j,(11-1)*3 + l) = dresdvar((11-1)*3 + j,(11-1)*3 + l) & + + (dPdF(i,j,k,l, 3) + dPdF(i,j,k,l, 7))*n(i,11)*n(k,11) & + + (dRdX(i,j,k,l,11) + dRdX(i,j,k,l,11))*n(i,11)*n(k,11) + dresdvar((11-1)*3 + j,( 2-1)*3 + l) = dresdvar((11-1)*3 + j,( 2-1)*3 + l) + dPdF(i,j,k,l, 3)*n(i,11)*n(k, 2) & + + dRdX(i,j,k,l,11)*n(i,11)*n(k, 2) + dresdvar((11-1)*3 + j,( 4-1)*3 + l) = dresdvar((11-1)*3 + j,( 4-1)*3 + l) - dPdF(i,j,k,l, 7)*n(i,11)*n(k, 4) & + - dRdX(i,j,k,l,11)*n(i,11)*n(k, 4) + dresdvar((11-1)*3 + j,( 5-1)*3 + l) = dresdvar((11-1)*3 + j,( 5-1)*3 + l) - dPdF(i,j,k,l, 3)*n(i,11)*n(k, 5) & + - dRdX(i,j,k,l,11)*n(i,11)*n(k, 5) + dresdvar((11-1)*3 + j,( 7-1)*3 + l) = dresdvar((11-1)*3 + j,( 7-1)*3 + l) + dPdF(i,j,k,l, 7)*n(i,11)*n(k, 7) & + + dRdX(i,j,k,l,11)*n(i,11)*n(k, 7) +! +! at boundary 12, influenced by boundary -2, +4, -6, +8 + dresdvar((12-1)*3 + j,(12-1)*3 + l) = dresdvar((12-1)*3 + j,(12-1)*3 + l) & + + (dPdF(i,j,k,l, 4) + dPdF(i,j,k,l, 8))*n(i,12)*n(k,12) & + + (dRdX(i,j,k,l,12) + dRdX(i,j,k,l,12))*n(i,12)*n(k,12) + dresdvar((12-1)*3 + j,( 2-1)*3 + l) = dresdvar((12-1)*3 + j,( 2-1)*3 + l) - dPdF(i,j,k,l, 4)*n(i,12)*n(k, 2) & + - dRdX(i,j,k,l,12)*n(i,12)*n(k, 2) + dresdvar((12-1)*3 + j,( 4-1)*3 + l) = dresdvar((12-1)*3 + j,( 4-1)*3 + l) + dPdF(i,j,k,l, 8)*n(i,12)*n(k, 4) & + + dRdX(i,j,k,l,12)*n(i,12)*n(k, 4) + dresdvar((12-1)*3 + j,( 6-1)*3 + l) = dresdvar((12-1)*3 + j,( 6-1)*3 + l) - dPdF(i,j,k,l, 4)*n(i,12)*n(k, 6) & + - dRdX(i,j,k,l,12)*n(i,12)*n(k, 6) + dresdvar((12-1)*3 + j,( 8-1)*3 + l) = dresdvar((12-1)*3 + j,( 8-1)*3 + l) + dPdF(i,j,k,l, 8)*n(i,12)*n(k, 8) & + + dRdX(i,j,k,l,12)*n(i,12)*n(k, 8) +! + enddo + enddo + enddo + enddo +! + return +! + END SUBROUTINE +! +! +!******************************************************************** +! Calculate the correction for the effective consisten tangent +!******************************************************************** + subroutine GIA_TangentCorrection(& + dF_grain,& ! deformation gradient increment of grains + ddF_corr) ! +! + implicit none +! + real(pReal), dimension(3,3,8) :: dF_grain,ddF_corr + integer(pInt) i,j,iBoun,grain +! + ddF_corr = 0.0_pReal + do j = 1,3 +! +! first relaxation direction + ddF_corr(1,j,1) = 2.0_pReal*abs(dF_grain(1,j,1))/(abs(dF_grain(1,j,1)) + abs(dF_grain(1,j,2))) + if (abs(dF_grain(1,j,1)) < 1.0e-8_pReal) ddF_corr(1,j,1) = 1.0_pReal + ddF_corr(1,j,2) = 2.0_pReal - ddF_corr(1,j,1) + ddF_corr(1,j,3) = 2.0_pReal*abs(dF_grain(1,j,3))/(abs(dF_grain(1,j,3)) + abs(dF_grain(1,j,4))) + if (abs(dF_grain(1,j,3)) < 1.0e-8_pReal) ddF_corr(1,j,3) = 1.0_pReal + ddF_corr(1,j,4) = 2.0_pReal - ddF_corr(1,j,3) + ddF_corr(1,j,5) = 2.0_pReal*abs(dF_grain(1,j,5))/(abs(dF_grain(1,j,5)) + abs(dF_grain(1,j,6))) + if (abs(dF_grain(1,j,5)) < 1.0e-8_pReal) ddF_corr(1,j,5) = 1.0_pReal + ddF_corr(1,j,6) = 2.0_pReal - ddF_corr(1,j,5) + ddF_corr(1,j,7) = 2.0_pReal*abs(dF_grain(1,j,7))/(abs(dF_grain(1,j,7)) + abs(dF_grain(1,j,8))) + if (abs(dF_grain(1,j,7)) < 1.0e-8_pReal) ddF_corr(1,j,7) = 1.0_pReal + ddF_corr(1,j,8) = 2.0_pReal - ddF_corr(1,j,7) +! +! second relaxation direction + ddF_corr(2,j,1) = 2.0_pReal*abs(dF_grain(2,j,1))/(abs(dF_grain(2,j,1)) + abs(dF_grain(2,j,3))) + if (abs(dF_grain(2,j,1)) < 1.0e-8_pReal) ddF_corr(2,j,1) = 1.0_pReal + ddF_corr(2,j,2) = 2.0_pReal*abs(dF_grain(2,j,2))/(abs(dF_grain(2,j,2)) + abs(dF_grain(2,j,4))) + if (abs(dF_grain(2,j,2)) < 1.0e-8_pReal) ddF_corr(2,j,2) = 1.0_pReal + ddF_corr(2,j,3) = 2.0_pReal - ddF_corr(2,j,1) + ddF_corr(2,j,4) = 2.0_pReal - ddF_corr(2,j,2) + ddF_corr(2,j,5) = 2.0_pReal*abs(dF_grain(2,j,5))/(abs(dF_grain(2,j,5)) + abs(dF_grain(2,j,7))) + if (abs(dF_grain(2,j,5)) < 1.0e-8_pReal) ddF_corr(2,j,5) = 1.0_pReal + ddF_corr(2,j,6) = 2.0_pReal*abs(dF_grain(2,j,6))/(abs(dF_grain(2,j,6)) + abs(dF_grain(2,j,8))) + if (abs(dF_grain(2,j,6)) < 1.0e-8_pReal) ddF_corr(2,j,6) = 1.0_pReal + ddF_corr(2,j,7) = 2.0_pReal - ddF_corr(2,j,5) + ddF_corr(2,j,8) = 2.0_pReal - ddF_corr(2,j,6) +! +! third relaxation direction + ddF_corr(3,j,1) = 2.0_pReal*abs(dF_grain(3,j,1))/(abs(dF_grain(3,j,1)) + abs(dF_grain(3,j,5))) + if (abs(dF_grain(3,j,1)) < 1.0e-8_pReal) ddF_corr(3,j,1) = 1.0_pReal + ddF_corr(3,j,2) = 2.0_pReal*abs(dF_grain(3,j,2))/(abs(dF_grain(3,j,2)) + abs(dF_grain(3,j,6))) + if (abs(dF_grain(3,j,2)) < 1.0e-8_pReal) ddF_corr(3,j,2) = 1.0_pReal + ddF_corr(3,j,3) = 2.0_pReal*abs(dF_grain(3,j,3))/(abs(dF_grain(3,j,3)) + abs(dF_grain(3,j,7))) + if (abs(dF_grain(3,j,3)) < 1.0e-8_pReal) ddF_corr(3,j,3) = 1.0_pReal + ddF_corr(3,j,4) = 2.0_pReal*abs(dF_grain(3,j,4))/(abs(dF_grain(3,j,4)) + abs(dF_grain(3,j,8))) + if (abs(dF_grain(3,j,4)) < 1.0e-8_pReal) ddF_corr(3,j,4) = 1.0_pReal + ddF_corr(3,j,5) = 2.0_pReal - ddF_corr(3,j,1) + ddF_corr(3,j,6) = 2.0_pReal - ddF_corr(3,j,2) + ddF_corr(3,j,7) = 2.0_pReal - ddF_corr(3,j,3) + ddF_corr(3,j,8) = 2.0_pReal - ddF_corr(3,j,4) +! + enddo +! + return +! + END SUBROUTINE +! + END MODULE +!############################################################## + diff --git a/trunk/CPFEM_Taylor.f90 b/trunk/CPFEM_Taylor.f90 new file mode 100644 index 000000000..b6f92a883 --- /dev/null +++ b/trunk/CPFEM_Taylor.f90 @@ -0,0 +1,290 @@ +!############################################################## + MODULE CPFEM +!############################################################## +! *** CPFEM engine *** +! + use prec, only: pReal,pInt + implicit none +! +! **************************************************************** +! *** General variables for the material behaviour calculation *** +! **************************************************************** + real(pReal), dimension (:,:), allocatable :: CPFEM_Temperature + real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_ffn_bar + real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_ffn1_bar + real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_PK1_bar + real(pReal), dimension (:,:,:,:,:,:),allocatable :: CPFEM_dPdF_bar + real(pReal), dimension (:,:,:), allocatable :: CPFEM_stress_bar + real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_jaco_bar + real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_jaco_knownGood + real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_results + real(pReal), dimension (:,:,:,:,:), allocatable :: CPFEM_Fp_old + real(pReal), dimension (:,:,:,:,:), allocatable :: CPFEM_Fp_new + real(pReal), parameter :: CPFEM_odd_stress = 1e15_pReal, CPFEM_odd_jacobian = 1e50_pReal + integer(pInt) :: CPFEM_Nresults = 4_pInt ! three Euler angles plus volume fraction + logical :: CPFEM_init_done = .false. ! remember if init has been done already + logical :: CPFEM_calc_done = .false. ! remember if first IP has already calced the results +! + CONTAINS +! +!********************************************************* +!*** allocate the arrays defined in module CPFEM *** +!*** and initialize them *** +!********************************************************* + SUBROUTINE CPFEM_init(Temperature) +! + use prec + use math, only: math_EulertoR, math_I3, math_identity2nd + use mesh + use constitutive +! + implicit none +! + real(pReal) Temperature + integer(pInt) e,i,g +! +! *** mpie.marc parameters *** + allocate(CPFEM_Temperature (mesh_maxNips,mesh_NcpElems)) ; CPFEM_Temperature = Temperature + allocate(CPFEM_ffn_bar (3,3,mesh_maxNips,mesh_NcpElems)) + forall(e=1:mesh_NcpElems,i=1:mesh_maxNips) CPFEM_ffn_bar(:,:,i,e) = math_I3 + allocate(CPFEM_ffn1_bar (3,3,mesh_maxNips,mesh_NcpElems)) ; CPFEM_ffn1_bar = CPFEM_ffn_bar + allocate(CPFEM_PK1_bar (3,3,mesh_maxNips,mesh_NcpElems)) ; CPFEM_PK1_bar = 0.0_pReal + allocate(CPFEM_dPdF_bar(3,3,3,3,mesh_maxNips,mesh_NcpElems)) ; CPFEM_dPdF_bar = 0.0_pReal + allocate(CPFEM_stress_bar(6,mesh_maxNips,mesh_NcpElems)) ; CPFEM_stress_bar = 0.0_pReal + allocate(CPFEM_jaco_bar(6,6,mesh_maxNips,mesh_NcpElems)) ; CPFEM_jaco_bar = 0.0_pReal + allocate(CPFEM_jaco_knownGood(6,6,mesh_maxNips,mesh_NcpElems)) ; CPFEM_jaco_knownGood = 0.0_pReal +! +! *** User defined results !!! MISSING incorporate consti_Nresults *** + allocate(CPFEM_results(CPFEM_Nresults+constitutive_maxNresults,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems)) + CPFEM_results = 0.0_pReal +! +! *** Plastic deformation gradient at (t=t0) and (t=t1) *** + allocate(CPFEM_Fp_new(3,3,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; CPFEM_Fp_new = 0.0_pReal + allocate(CPFEM_Fp_old(3,3,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems)) + forall (e=1:mesh_NcpElems,i=1:mesh_maxNips,g=1:constitutive_maxNgrains) & + CPFEM_Fp_old(:,:,g,i,e) = math_EulerToR(constitutive_EulerAngles(:,g,i,e)) ! plastic def gradient reflects init orientation +! +! *** Output to MARC output file *** + write(6,*) + write(6,*) 'CPFEM Initialization' + write(6,*) + write(6,*) 'CPFEM_Temperature: ', shape(CPFEM_Temperature) + write(6,*) 'CPFEM_ffn_bar: ', shape(CPFEM_ffn_bar) + write(6,*) 'CPFEM_ffn1_bar: ', shape(CPFEM_ffn1_bar) + write(6,*) 'CPFEM_PK1_bar: ', shape(CPFEM_PK1_bar) + write(6,*) 'CPFEM_dPdF_bar: ', shape(CPFEM_dPdF_bar) + write(6,*) 'CPFEM_stress_bar: ', shape(CPFEM_stress_bar) + write(6,*) 'CPFEM_jaco_bar: ', shape(CPFEM_jaco_bar) + write(6,*) 'CPFEM_jaco_knownGood: ', shape(CPFEM_jaco_knownGood) + write(6,*) 'CPFEM_results: ', shape(CPFEM_results) + write(6,*) 'CPFEM_Fp_old: ', shape(CPFEM_Fp_old) + write(6,*) 'CPFEM_Fp_new: ', shape(CPFEM_Fp_new) + write(6,*) + call flush(6) + return +! + END SUBROUTINE +! +! +!*********************************************************************** +!*** perform initialization at first call, update variables and *** +!*** call the actual material model *** +! +! CPFEM_mode computation mode (regular, collection, recycle) +! ffn deformation gradient for t=t0 +! ffn1 deformation gradient for t=t1 +! Temperature temperature +! CPFEM_dt time increment +! CPFEM_en element number +! CPFEM_in intergration point number +! CPFEM_stress stress vector in Mandel notation +! CPFEM_updateJaco flag to initiate computation of Jacobian +! CPFEM_jaco jacobian in Mandel notation +! CPFEM_ngens size of stress strain law +!*********************************************************************** + SUBROUTINE CPFEM_general(CPFEM_mode, ffn, ffn1, Temperature, CPFEM_dt,& + CPFEM_en, CPFEM_in, CPFEM_stress, CPFEM_updateJaco, CPFEM_jaco, CPFEM_ngens) +! note: CPFEM_stress = Cauchy stress cs(6) and CPFEM_jaco = Consistent tangent dcs/de +! + use prec, only: pReal,pInt + use FEsolving + use debug + use math, only: math_init, invnrmMandel, math_identity2nd, math_Mandel3333to66,math_Mandel33to6,math_Mandel6to33,math_det3x3,math_I3 + use mesh, only: mesh_init,mesh_FEasCP, mesh_NcpElems, FE_Nips, FE_mapElemtype, mesh_element + use lattice, only: lattice_init + use constitutive, only: constitutive_init,constitutive_state_old,constitutive_state_new,material_Cslip_66 + implicit none +! + integer(pInt) CPFEM_en, CPFEM_in, cp_en, CPFEM_ngens, i,j,k,l,m,n, e + real(pReal), dimension (3,3) :: ffn,ffn1,Kirchhoff_bar + real(pReal), dimension (3,3,3,3) :: H_bar + real(pReal), dimension(CPFEM_ngens) :: CPFEM_stress + real(pReal), dimension(CPFEM_ngens,CPFEM_ngens) :: CPFEM_jaco + real(pReal) Temperature,CPFEM_dt,J_inverse + integer(pInt) CPFEM_mode ! 1: regular computation with aged results& + ! 2: regular computation& + ! 3: collection of FEM data& + ! 4: recycling of former results (MARC speciality)& + ! 5: record tangent from former converged inc& + ! 6: restore tangent from former converged inc + logical CPFEM_updateJaco +! + if (.not. CPFEM_init_done) then ! initialization step (three dimensional stress state check missing?) + call math_init() + call mesh_init() + call lattice_init() + call constitutive_init() + call CPFEM_init(Temperature) + CPFEM_init_done = .true. + endif +! + cp_en = mesh_FEasCP('elem',CPFEM_en) + if (cp_en == 1 .and. CPFEM_in == 1) & + write(6,'(a6,x,i4,x,a4,x,i4,x,a10,x,f8.4,x,a10,x,i2,x,a10,x,i2,x,a10,x,i2,x,a10,x,i2)') & + 'elem',cp_en,'IP',CPFEM_in,& + 'theTime',theTime,'theInc',theInc,'theCycle',theCycle,'theLovl',theLovl,& + 'mode',CPFEM_mode +! + select case (CPFEM_mode) + case (2,1) ! regular computation (with aging of results) + if (.not. CPFEM_calc_done) then ! puuh, me needs doing all the work... + write (6,*) 'puuh me needs doing all the work', cp_en + if (CPFEM_mode == 1) then ! age results at start of new increment + CPFEM_Fp_old = CPFEM_Fp_new + constitutive_state_old = constitutive_state_new + write (6,*) '#### aged results' + endif + debug_cutbackDistribution = 0_pInt ! initialize debugging data + debug_InnerLoopDistribution = 0_pInt + debug_OuterLoopDistribution = 0_pInt +! + do e=1,mesh_NcpElems ! ## this shall be done in a parallel loop in the future ## + do i=1,FE_Nips(mesh_element(2,e)) ! iterate over all IPs of this element's type + debugger = (e==1 .and. i==1) ! switch on debugging for first IP in first element + call CPFEM_MaterialPoint(CPFEM_updateJaco, CPFEM_dt, i, e) + enddo + enddo + call debug_info() ! output of debugging/performance statistics + CPFEM_calc_done = .true. ! now calc is done + endif +! translate from P and dP/dF to CS and dCS/dE + Kirchhoff_bar = matmul(CPFEM_PK1_bar(:,:,CPFEM_in, cp_en),transpose(CPFEM_ffn1_bar(:,:,CPFEM_in, cp_en))) + J_inverse = 1.0_pReal/math_det3x3(CPFEM_ffn1_bar(:,:,CPFEM_in, cp_en)) + CPFEM_stress_bar(1:CPFEM_ngens,CPFEM_in,cp_en) = math_Mandel33to6(J_inverse*Kirchhoff_bar) +! + H_bar = 0.0_pReal + forall(i=1:3,j=1:3,k=1:3,l=1:3,m=1:3,n=1:3) & + H_bar(i,j,k,l) = H_bar(i,j,k,l) + & + (CPFEM_ffn1_bar(j,m,CPFEM_in,cp_en)*CPFEM_ffn1_bar(l,n,CPFEM_in,cp_en)*CPFEM_dPdF_bar(i,m,k,n,CPFEM_in,cp_en) - & + math_I3(j,l)*CPFEM_ffn1_bar(i,m,CPFEM_in,cp_en)*CPFEM_PK1_bar(k,m,CPFEM_in,cp_en)) + & + 0.5_pReal*(math_I3(i,k)*Kirchhoff_bar(j,l) + math_I3(j,l)*Kirchhoff_bar(i,k) + & + math_I3(i,l)*Kirchhoff_bar(j,k) + math_I3(j,k)*Kirchhoff_bar(i,l)) + CPFEM_jaco_bar(1:CPFEM_ngens,1:CPFEM_ngens,CPFEM_in,cp_en) = math_Mandel3333to66(J_inverse*H_bar) +! + case (3) ! collect and return odd result + CPFEM_Temperature(CPFEM_in,cp_en) = Temperature + CPFEM_ffn_bar(:,:,CPFEM_in,cp_en) = ffn + CPFEM_ffn1_bar(:,:,CPFEM_in,cp_en) = ffn1 + CPFEM_stress_bar(1:CPFEM_ngens,CPFEM_in,cp_en) = CPFEM_odd_stress + CPFEM_jaco_bar(1:CPFEM_ngens,1:CPFEM_ngens,CPFEM_in,cp_en) = CPFEM_odd_jacobian*math_identity2nd(CPFEM_ngens) + CPFEM_calc_done = .false. + + case (4) ! do nothing since we can recycle the former results (MARC specialty) + case (5) ! record consistent tangent at beginning of new increment + CPFEM_jaco_knownGood = CPFEM_jaco_bar + case (6) ! restore consistent tangent after cutback + CPFEM_jaco_bar = CPFEM_jaco_knownGood + end select +! +! return the local stress and the jacobian from storage + CPFEM_stress(1:CPFEM_ngens) = CPFEM_stress_bar(1:CPFEM_ngens,CPFEM_in,cp_en) + CPFEM_jaco(1:CPFEM_ngens,1:CPFEM_ngens) = CPFEM_jaco_bar(1:CPFEM_ngens,1:CPFEM_ngens,CPFEM_in,cp_en) + if (cp_en == 1 .and. CPFEM_in == 1) write (6,*) 'stress',CPFEM_stress + if (cp_en == 1 .and. CPFEM_in == 1 .and. CPFEM_updateJaco) write (6,*) 'stiffness',CPFEM_jaco +! + return +! + END SUBROUTINE +! +!********************************************************** +!*** calculate the material point behaviour *** +!********************************************************** + SUBROUTINE CPFEM_MaterialPoint(& + updateJaco,& ! flag to initiate Jacobian updating + CPFEM_dt,& ! Time increment (dt) + CPFEM_in,& ! Integration point number + cp_en) ! Element number +! + use prec + use debug + use math, only: math_pDecomposition,math_RtoEuler,inDeg,math_I3,math_invert3x3,math_permut,math_invert,math_delta + use IO, only: IO_error + use mesh, only: mesh_element + use crystallite + use constitutive + implicit none +! + character(len=128) msg + integer(pInt) cp_en,CPFEM_in,grain,max_cutbacks,i,j,k,l,m,n + logical updateJaco,error + real(pReal) CPFEM_dt,volfrac + real(pReal), dimension(3,3) :: F0_bar,F1_bar,dF_bar + real(pReal), dimension(3,3) :: U,R + real(pReal), dimension(3,3,8) :: PK1,Fp0,Fp1,Fe1 + real(pReal), dimension(3,3,3,3,8) :: dPdF + real(pReal), dimension(constitutive_maxNstatevars,8) :: state0,state1 +! + CPFEM_PK1_bar(:,:,CPFEM_in,cp_en) = 0.0_pReal ! zero out average first PK stress + if (updateJaco) CPFEM_dPdF_bar(:,:,:,:,CPFEM_in,cp_en) = 0.0_pReal ! zero out average consistent tangent +! + F0_bar = CPFEM_ffn_bar(:,:,CPFEM_in,cp_en) ! effective deformation gradient at t_n + F1_bar = CPFEM_ffn1_bar(:,:,CPFEM_in,cp_en) ! effective deformation gradient at t_n+1 + state0 = constitutive_state_old(:,:,CPFEM_in,cp_en) ! state variables at t_n + Fp0 = CPFEM_Fp_old(:,:,:,CPFEM_in,cp_en) ! grain plastic def. gradient at t_n +! +! -------------- grain loop ----------------- + do grain = 1,texture_Ngrains(mesh_element(4,cp_en)) + call SingleCrystallite(msg,PK1(:,:,grain),dPdF(:,:,:,:,grain),& + CPFEM_results(5:4+constitutive_Nresults(grain,CPFEM_in,cp_en),grain,CPFEM_in,cp_en),& + Fp1(:,:,grain),Fe1(:,:,grain),state1(:,grain),& ! output up to here + CPFEM_dt,cp_en,CPFEM_in,grain,.true.,& + CPFEM_Temperature(CPFEM_in,cp_en),F1_bar,F0_bar,Fp0(:,:,grain),state0(:,grain)) + if (msg /= 'ok') then ! solution not reached + call IO_error(600) + return + endif + +! update results plotted in MENTAT + call math_pDecomposition(Fe1(:,:,grain),U,R,error) ! polar decomposition + if (error) then + write(6,*) Fe1(:,:,grain) + write(6,*) 'polar decomposition' + write(6,*) 'Grain: ',grain + write(6,*) 'Integration point: ',CPFEM_in + write(6,*) 'Element: ',mesh_element(1,cp_en) + call IO_error(650) + return + endif +! + volfrac = constitutive_matVolFrac(grain,CPFEM_in,cp_en)*constitutive_texVolFrac(grain,CPFEM_in,cp_en) + CPFEM_PK1_bar(:,:,CPFEM_in,cp_en) = CPFEM_PK1_bar(:,:,CPFEM_in,cp_en) + & + volfrac*PK1(:,:,grain) ! average Cauchy stress + if (updateJaco) then ! consistent tangent + CPFEM_dPdF_bar(:,:,:,:,CPFEM_in,cp_en) = CPFEM_dPdF_bar(:,:,:,:,CPFEM_in,cp_en) + & + volfrac*dPdF(:,:,:,:,grain) + endif + CPFEM_results(1:3,grain,CPFEM_in,cp_en) = math_RtoEuler(transpose(R))*inDeg ! orientation + CPFEM_results(4 ,grain,CPFEM_in,cp_en) = volfrac ! volume fraction of orientation + enddo ! grain loop +! +! updates all variables, deformation gradients, and vectors + CPFEM_Fp_new(:,:,:,CPFEM_in,cp_en) = Fp1 + constitutive_state_new(:,:,CPFEM_in,cp_en) = state1 +! + return +! + END SUBROUTINE +! + END MODULE +!############################################################## + diff --git a/trunk/FEsolving.f90 b/trunk/FEsolving.f90 index aacbd802a..63fcb33e7 100644 --- a/trunk/FEsolving.f90 +++ b/trunk/FEsolving.f90 @@ -3,13 +3,12 @@ MODULE FEsolving !############################################################## - use prec, only: pInt + use prec, only: pInt,pReal implicit none integer(pInt) cycleCounter - integer(pInt) theInc, theCycle, theLovl - logical :: outdatedByNewInc = .false. - - + integer(pInt) theInc,theCycle,theLovl + real(pReal) theTime + logical :: lastIncConverged = .false.,outdatedByNewInc = .false. END MODULE FEsolving diff --git a/trunk/IO.f90 b/trunk/IO.f90 index 85eb8ff7c..406fecdad 100644 --- a/trunk/IO.f90 +++ b/trunk/IO.f90 @@ -595,6 +595,8 @@ msg='Polar decomposition failed' case (700) msg='Singular matrix in stress iteration' + case (800) + msg='GIA requires 8 grains per IP (bonehead, you!)' case default msg='Unknown error number' end select diff --git a/trunk/constitutive_pheno.f90 b/trunk/constitutive_pheno.f90 index 8314cefd8..1d004e96d 100644 --- a/trunk/constitutive_pheno.f90 +++ b/trunk/constitutive_pheno.f90 @@ -271,7 +271,7 @@ do while(.true.) else if (section>0) then select case(tag) - case ('crystal_structure') + case ('lattice_structure') material_CrystalStructure(section)=IO_intValue(line,positions,2) case ('nslip') material_Nslip(section)=IO_intValue(line,positions,2) @@ -413,7 +413,7 @@ subroutine constitutive_Parse_MatTexDat(filename) use prec, only: pReal,pInt use IO, only: IO_error, IO_open_file use math, only: math_Mandel3333to66, math_Voigt66to3333 -use crystal, only: crystal_SlipIntType +use lattice, only: lattice_SlipIntType implicit none !* Definition of variables @@ -454,7 +454,7 @@ allocate(material_n_slip(material_maxN)) ; material_n_slip=0.0_pReal allocate(material_h0(material_maxN)) ; material_h0=0.0_pReal allocate(material_s_sat(material_maxN)) ; material_s_sat=0.0_pReal allocate(material_w0(material_maxN)) ; material_w0=0.0_pReal -allocate(material_SlipIntCoeff(maxval(crystal_SlipIntType),material_maxN)) ; material_SlipIntCoeff=0.0_pReal +allocate(material_SlipIntCoeff(maxval(lattice_SlipIntType),material_maxN)) ; material_SlipIntCoeff=0.0_pReal allocate(material_GrainSize(material_maxN)) ; material_GrainSize=0.0_pReal allocate(material_bg(material_maxN)) ; material_bg=0.0_pReal allocate(texture_ODFfile(texture_maxN)) ; texture_ODFfile='' @@ -546,7 +546,7 @@ use prec, only: pReal,pInt use math, only: math_sampleGaussOri,math_sampleFiberOri,math_sampleRandomOri,math_symmetricEulers,math_EulerToR use mesh, only: mesh_NcpElems,FE_Nips,mesh_maxNips,mesh_element use IO, only: IO_hybridIA -use crystal, only: crystal_SlipIntType +use lattice, only: lattice_SlipIntType implicit none @@ -688,7 +688,7 @@ do i=1,material_maxN do j=1,material_Nslip(i) do k=1,material_Nslip(i) !* min function is used to distinguish self hardening from latent hardening - constitutive_HardeningMatrix(k,j,i) = material_SlipIntCoeff(max(2,min(3,crystal_SlipIntType(k,j,i)))-1,i) ! 1,2,3,4,5 --> 1,1,2,2,2 + constitutive_HardeningMatrix(k,j,i) = material_SlipIntCoeff(max(2,min(3,lattice_SlipIntType(k,j,i)))-1,i) ! 1,2,3,4,5 --> 1,1,2,2,2 enddo enddo enddo @@ -759,7 +759,7 @@ subroutine constitutive_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,state,Temperature, !* - dLp_dTstar : derivative of Lp (4th-order tensor) * !********************************************************************* use prec, only: pReal,pInt -use crystal, only: crystal_Sslip,crystal_Sslip_v +use lattice, only: lattice_Sslip,lattice_Sslip_v use math, only: math_Plain3333to99 use debug @@ -782,10 +782,10 @@ matID = constitutive_matID(ipc,ip,el) !* Calculation of Lp Lp = 0.0_pReal do i=1,material_Nslip(matID) - tau_slip(i)=dot_product(Tstar_v,crystal_Sslip_v(:,i,material_CrystalStructure(matID))) + tau_slip(i)=dot_product(Tstar_v,lattice_Sslip_v(:,i,material_CrystalStructure(matID))) gdot_slip(i)=material_gdot0_slip(matID)*(abs(tau_slip(i))/state(i))**& material_n_slip(matID)*sign(1.0_pReal,tau_slip(i)) - Lp=Lp+gdot_slip(i)*crystal_Sslip(:,:,i,material_CrystalStructure(matID)) + Lp=Lp+gdot_slip(i)*lattice_Sslip(:,:,i,material_CrystalStructure(matID)) enddo !* Calculation of the tangent of Lp @@ -796,8 +796,8 @@ do i=1,material_Nslip(matID) (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_dTstar3333(k,l,m,n) = dLp_dTstar3333(k,l,m,n) + & - dgdot_dtauslip(i)*crystal_Sslip(k,l,i,material_CrystalStructure(matID))* & - crystal_Sslip(m,n,i,material_CrystalStructure(matID)) + dgdot_dtauslip(i)*lattice_Sslip(k,l,i,material_CrystalStructure(matID))* & + lattice_Sslip(m,n,i,material_CrystalStructure(matID)) enddo dLp_dTstar = math_Plain3333to99(dLp_dTstar3333) @@ -819,7 +819,7 @@ function constitutive_dotState(Tstar_v,state,Temperature,ipc,ip,el) !* - constitutive_dotState : evolution of state variable * !********************************************************************* use prec, only: pReal,pInt -use crystal, only: crystal_Sslip_v +use lattice, only: lattice_Sslip_v implicit none !* Definition of variables @@ -834,7 +834,7 @@ matID = constitutive_matID(ipc,ip,el) !* Self-Hardening of each system do i=1,constitutive_Nstatevars(ipc,ip,el) - tau_slip = dot_product(Tstar_v,crystal_Sslip_v(:,i,material_CrystalStructure(matID))) + tau_slip = dot_product(Tstar_v,lattice_Sslip_v(:,i,material_CrystalStructure(matID))) gdot_slip = material_gdot0_slip(matID)*(abs(tau_slip)/state(i))**& material_n_slip(matID)*sign(1.0_pReal,tau_slip) self_hardening(i)=material_h0(matID)*(1.0_pReal-state(i)/& @@ -848,7 +848,7 @@ return end function -function constitutive_post_results(Tstar_v,state,dt,Temperature,ipc,ip,el) +function constitutive_post_results(Tstar_v,state,Temperature,dt,ipc,ip,el) !********************************************************************* !* return array of constitutive results * !* INPUT: * @@ -860,7 +860,7 @@ function constitutive_post_results(Tstar_v,state,dt,Temperature,ipc,ip,el) !* - el : current element * !********************************************************************* use prec, only: pReal,pInt -use crystal, only: crystal_Sslip_v +use lattice, only: lattice_Sslip_v implicit none !* Definition of variables @@ -878,7 +878,7 @@ if(constitutive_Nresults(ipc,ip,el)==0) return do i=1,material_Nslip(matID) constitutive_post_results(i) = state(i) - tau_slip=dot_product(Tstar_v,crystal_Sslip_v(:,i,material_CrystalStructure(matID))) + tau_slip=dot_product(Tstar_v,lattice_Sslip_v(:,i,material_CrystalStructure(matID))) constitutive_post_results(i+material_Nslip(matID)) = & dt*material_gdot0_slip(matID)*(abs(tau_slip)/state(i))**material_n_slip(matID)*sign(1.0_pReal,tau_slip) enddo diff --git a/trunk/crystallite.f90 b/trunk/crystallite.f90 new file mode 100644 index 000000000..7c8b37d2e --- /dev/null +++ b/trunk/crystallite.f90 @@ -0,0 +1,307 @@ +!############################################################## + MODULE crystallite +!############################################################## +! *** Solution at single crystallite level *** +! +CONTAINS +! +! +!******************************************************************** +! Calculates the stress for a single component +!******************************************************************** +!*********************************************************************** +!*** calculation of stress (P), stiffness (dPdF), *** +!*** and announcment of any *** +!*** acceleration of the Newton-Raphson correction *** +!*********************************************************************** + subroutine SingleCrystallite(& + msg,& ! return message + P,& ! first PK stress + dPdF,& ! consistent tangent + post_results,& ! plot results from constitutive model + Fp_new,& ! new plastic deformation gradient + Fe_new,& ! new "elastic" deformation gradient + state_new,& ! new state variable array +! + dt,& ! time increment + cp_en,& ! element number + ip,& ! integration point number + grain,& ! grain number + updateJaco,& ! update of Jacobian required + Temperature,& ! temperature of crystallite + Fg_new,& ! new global deformation gradient + Fg_old,& ! old global deformation gradient + Fp_old,& ! old plastic deformation gradient + state_old) ! old state variable array +! + use prec, only: pReal,pInt,pert_Fg,subStepMin + use debug + use constitutive, only: constitutive_Nstatevars,constitutive_Nresults + use mesh, only: mesh_element + use math, only: math_Mandel6to33,math_Mandel33to6,math_Mandel3333to66,& + math_I3,math_det3x3,math_invert3x3 + implicit none +! + character(len=*) msg + logical updateJaco,error,cuttedBack,guessNew + integer(pInt) cp_en,ip,grain,i,j,k,l,m,n, nCutbacks + real(pReal) Temperature + real(pReal) dt,dt_aim,subFrac,subStep,invJ,det + real(pReal), dimension(3,3) :: Lp,Lp_pert,inv + real(pReal), dimension(3,3) :: Fg_old,Fg_current,Fg_aim,Fg_new,Fg_pert,deltaFg + real(pReal), dimension(3,3) :: Fp_old,Fp_current,Fp_new,Fp_pert + real(pReal), dimension(3,3) :: Fe_old,Fe_current,Fe_new,Fe_pert + real(pReal), dimension(3,3) :: Tstar,tau,P,P_pert + real(pReal), dimension(3,3,3,3) :: dPdF + real(pReal), dimension(constitutive_Nstatevars(grain,ip,cp_en)) :: state_old,state_current,state_new,state_pert + real(pReal), dimension(constitutive_Nresults(grain,ip,cp_en)) :: post_results +! + deltaFg = Fg_new-Fg_old + subFrac = 0.0_pReal + subStep = 1.0_pReal + nCutbacks = 0_pInt +! + Fg_aim = Fg_old ! make "new", "aim" a synonym for "old" + Fp_new = Fp_old + call math_invert3x3(Fp_old,inv,det,error) + Fe_new = matmul(Fg_old,inv) + state_new = state_old +! + cuttedBack = .false. + guessNew = .true. +! +! begin the cutback loop + do while (subStep > subStepMin) ! continue until finished or too much cut backing + if (.not. cuttedBack) then + Fg_current = Fg_aim ! wind forward + Fp_current = Fp_new + Fe_current = Fe_new + state_current = state_new + endif +! + Fg_aim = Fg_current + subStep*deltaFg ! aim for Fg + dt_aim = subStep*dt ! aim for dt + msg = '' ! error free so far + if (guessNew) then ! calculate new Lp guess when cutted back + if (dt_aim /= 0.0_pReal) then + call math_invert3x3(Fg_aim,inv,det,error) + Lp = (math_I3-matmul(Fp_current,matmul(inv,Fe_current)))/dt ! fully plastic initial guess + else + Lp = 0.0_pReal ! fully elastic guess + endif + endif + call TimeIntegration(msg,Lp,Fp_new,Fe_new,P,state_new,post_results,.true., & ! def gradients and PK2 at end of time step + dt_aim,cp_en,ip,grain,Temperature,Fg_aim,Fp_current,state_current) +! + if (msg == 'ok') then + cuttedBack = .false. ! no cut back required + guessNew = .false. ! keep the Lp + subFrac = subFrac + subStep + subStep = 1.0_pReal - subFrac ! try one go + else + nCutbacks = nCutbacks + 1 ! record additional cutback + cuttedBack = .true. ! encountered problems --> + guessNew = .true. ! redo plastic Lp guess + subStep = subStep / 2.0_pReal ! cut time step in half + endif + enddo +! + debug_cutbackDistribution(nCutbacks+1) = debug_cutbackDistribution(nCutbacks+1)+1 +! + if (msg /= 'ok') return ! solution not reached --> report back + if (updateJaco) then ! consistent tangent using + do k=1,3 + do l=1,3 + Fg_pert = Fg_new ! initialize perturbed Fg + Fg_pert(k,l) = Fg_pert(k,l) + pert_Fg ! perturb single component + Lp_pert = Lp + state_pert = state_new ! initial guess from end of time step + call TimeIntegration(msg,Lp,Fp_pert,Fe_pert,P_pert,state_pert,post_results,.false., & ! def gradients and PK2 at end of time step + dt_aim,cp_en,ip,grain,Temperature,Fg_pert,Fp_current,state_current) + if (msg /= 'ok') then + msg = 'consistent tangent --> '//msg + return + endif + dPdF(:,:,k,l) = (P_pert-P)/pert_Fg ! constructing the tangent dP_ij/dFg_kl from forward differences + enddo + enddo + endif +! + return +! + END SUBROUTINE +! +!*********************************************************************** +!*** fully-implicit two-level time integration *** +!*** based on a residuum in Lp and intermediate *** +!*** acceleration of the Newton-Raphson correction *** +!*********************************************************************** + SUBROUTINE TimeIntegration(& + msg,& ! return message + Lpguess,& ! guess of plastic velocity gradient + Fp_new,& ! new plastic deformation gradient + Fe_new,& ! new "elastic" deformation gradient + P,& ! 1nd PK stress (taken as initial guess if /= 0) + state,& ! current microstructure at end of time inc (taken as guess if /= 0) + results,& ! post results from constitutive + wantsConstitutiveResults,& ! its flag +! + dt,& ! time increment + cp_en,& ! element number + ip,& ! integration point number + grain,& ! grain number + Temperature,& ! temperature + Fg_new,& ! new total def gradient + Fp_old,& ! former plastic def gradient + state_old) ! former microstructure +! + use prec + use debug + use mesh, only: mesh_element + use constitutive, only: constitutive_Nstatevars,& + constitutive_homogenizedC,constitutive_dotState,constitutive_LpAndItsTangent,& + constitutive_Nresults,constitutive_Microstructure,constitutive_post_results + use math + implicit none +! + character(len=*) msg + logical failed,wantsConstitutiveResults + integer(pInt) cp_en, ip, grain + integer(pInt) iOuter,iInner,dummy, i,j,k,l,m,n + real(pReal) dt, Temperature, det, p_hydro, leapfrog,maxleap + real(pReal), dimension(6) :: Tstar_v + real(pReal), dimension(9,9) :: dLp,dTdLp,dRdLp,invdRdLp,eye2 + real(pReal), dimension(6,6) :: C_66 + real(pReal), dimension(3,3) :: Fg_new,Fp_new,invFp_new,Fp_old,invFp_old,Fe_new,Fe_old + real(pReal), dimension(3,3) :: P,Tstar + real(pReal), dimension(3,3) :: Lp,Lpguess,Lpguess_old,Rinner,Rinner_old,A,B,BT,AB,BTA + real(pReal), dimension(3,3,3,3) :: C + real(pReal), dimension(constitutive_Nstatevars(grain, ip, cp_en)) :: state_old,state,ROuter + real(pReal), dimension(constitutive_Nresults(grain,ip,cp_en)) :: results +! + msg = 'ok' ! error-free so far + eye2 = math_identity2nd(9) + + call math_invert3x3(Fp_old,invFp_old,det,failed) ! inversion of Fp_old + if (failed) then + msg = 'inversion Fp_old' + return + endif + + A = matmul(transpose(invFp_old), matmul(transpose(Fg_new),matmul(Fg_new,invFp_old))) +! + if (all(state == 0.0_pReal)) state = state_old ! former state guessed, if none specified + iOuter = 0_pInt ! outer counter +! +! +Outer: do ! outer iteration: State + iOuter = iOuter+1 + if (iOuter > nOuter) then + msg = 'limit Outer iteration' + debug_OuterLoopDistribution(nOuter) = debug_OuterLoopDistribution(nOuter)+1 + return + endif + call constitutive_Microstructure(state,Temperature,grain,ip,cp_en) + C_66 = constitutive_HomogenizedC(state, grain, ip, cp_en) + C = math_Mandel66to3333(C_66) ! 4th rank elasticity tensor +! + iInner = 0_pInt + leapfrog = 1.0_pReal ! correction as suggested by invdRdLp-step + maxleap = 1024.0_pReal ! preassign maximum acceleration level +! +Inner: do ! inner iteration: Lp + iInner = iInner+1 + if (iInner > nInner) then ! too many loops required + msg = 'limit Inner iteration' + debug_InnerLoopDistribution(nInner) = debug_InnerLoopDistribution(nInner)+1 + return + endif + B = math_i3 - dt*Lpguess + BT = transpose(B) + AB = matmul(A,B) + BTA = matmul(BT,A) + Tstar_v = 0.5_pReal*matmul(C_66,math_mandel33to6(matmul(BT,AB)-math_I3)) + Tstar = math_Mandel6to33(Tstar_v) + p_hydro=(Tstar_v(1)+Tstar_v(2)+Tstar_v(3))/3.0_pReal + forall(i=1:3) Tstar_v(i) = Tstar_v(i)-p_hydro ! subtract hydrostatic pressure + call constitutive_LpAndItsTangent(Lp,dLp, & + Tstar_v,state,Temperature,grain,ip,cp_en) + Rinner = Lpguess - Lp ! update current residuum + if ((maxval(abs(Rinner)) < abstol_Inner) .or. & + (any(abs(dt*Lpguess) > relevantStrain) .and. & + maxval(abs(Rinner/Lpguess),abs(dt*Lpguess) > relevantStrain) < reltol_Inner)) & + exit Inner +! +! check for acceleration/deceleration in Newton--Raphson correction + if (leapfrog > 1.0_pReal .and. & + (sum(Rinner*Rinner) > sum(Rinner_old*Rinner_old) .or. & ! worse residuum + sum(Rinner*Rinner_old) < 0.0_pReal)) then ! residuum changed sign (overshoot) + maxleap = 0.5_pReal * leapfrog ! limit next acceleration + leapfrog = 1.0_pReal ! grinding halt + else ! better residuum + dTdLp = 0.0_pReal ! calc dT/dLp + forall (i=1:3,j=1:3,k=1:3,l=1:3,m=1:3,n=1:3) & + dTdLp(3*(i-1)+j,3*(k-1)+l) = dTdLp(3*(i-1)+j,3*(k-1)+l) + & + C(i,j,l,n)*AB(k,n)+C(i,j,m,l)*BTA(m,k) + dTdLp = -0.5_pReal*dt*dTdLp + dRdLp = eye2 - matmul(dLp,dTdLp) ! calc dR/dLp + invdRdLp = 0.0_pReal + call math_invert(9,dRdLp,invdRdLp,dummy,failed) ! invert dR/dLp --> dLp/dR + if (failed) then + msg = 'inversion dR/dLp' + if (debugger) then + write (6,*) msg + write (6,'(a,/,9(9(e9.3,x)/))') 'dRdLp', dRdLp(1:9,:) + write (6,*) 'state',state + write (6,'(a,/,3(3(f12.7,x)/))') 'Lpguess',Lpguess(1:3,:) + write (6,*) 'Tstar',Tstar_v + endif + return + endif +! + Rinner_old = Rinner ! remember current residuum + Lpguess_old = Lpguess ! remember current Lp guess + if (iInner > 1 .and. leapfrog < maxleap) leapfrog = 2.0_pReal * leapfrog ! accelerate + endif +! + Lpguess = Lpguess_old ! start from current guess + Rinner = Rinner_old ! use current residuum + forall (i=1:3,j=1:3,k=1:3,l=1:3) & ! leapfrog to updated Lpguess + Lpguess(i,j) = Lpguess(i,j) - leapfrog*invdRdLp(3*(i-1)+j,3*(k-1)+l)*Rinner(k,l) + enddo Inner +! + debug_InnerLoopDistribution(iInner) = debug_InnerLoopDistribution(iInner)+1 + ROuter = state - state_old - & + dt*constitutive_dotState(Tstar_v,state,Temperature,& + grain,ip,cp_en) ! residuum from evolution of microstructure + state = state - ROuter ! update of microstructure + if (maxval(abs(Router/state),state /= 0.0_pReal) < reltol_Outer) exit Outer + enddo Outer +! + debug_OuterLoopDistribution(iOuter) = debug_OuterLoopDistribution(iOuter)+1 + + + invFp_new = matmul(invFp_old,B) + call math_invert3x3(invFp_new,Fp_new,det,failed) + if (failed) then + msg = 'inversion Fp_new^-1' + return + endif +! + if (wantsConstitutiveResults) then ! get the post_results upon request + results = 0.0_pReal + results = constitutive_post_results(Tstar_v,state,Temperature,dt,grain,ip,cp_en) + endif +! + Fp_new = Fp_new*det**(1.0_pReal/3.0_pReal) ! regularize Fp by det = det(InvFp_new) !! + Fe_new = matmul(Fg_new,invFp_new) ! calc resulting Fe + forall (i=1:3) Tstar_v(i) = Tstar_v(i)+p_hydro ! add hydrostatic component back + P = matmul(Fe_new,matmul(Tstar,transpose(invFp_new))) ! first PK stress +! + return +! + END SUBROUTINE +!! + END MODULE +!############################################################## + diff --git a/trunk/lattice.f90 b/trunk/lattice.f90 new file mode 100644 index 000000000..6463556c2 --- /dev/null +++ b/trunk/lattice.f90 @@ -0,0 +1,370 @@ + +!************************************ +!* Module: LATTICE * +!************************************ +!* contains: * +!* - Lattice structure definition * +!* - Slip system definition * +!* - Schmid matrices calculation * +!************************************ + +MODULE lattice + +!*** Include other modules *** +use prec, only: pReal,pInt +implicit none + +!************************************ +!* Lattice structures * +!************************************ +!* Number of lattice structures (1-FCC,2-BCC,3-HCP) +integer(pInt), parameter :: lattice_MaxLatticeStructure = 3 +!* Total number of slip systems per lattice structure +!* (has to be changed according the definition of slip systems) +integer(pInt), dimension(lattice_MaxLatticeStructure), parameter :: lattice_MaxNslipOfStructure = & +reshape((/12,48,12/),(/lattice_MaxLatticeStructure/)) +!* Total number of twin systems per lattice structure +!* (has to be changed according the definition of twin systems) +integer(pInt), dimension(lattice_MaxLatticeStructure), parameter :: lattice_MaxNtwinOfStructure = & +reshape((/12,12,6/),(/lattice_MaxLatticeStructure/)) +!* Maximum number of slip systems over lattice structures +integer(pInt), parameter :: lattice_MaxMaxNslipOfStructure = 48 +!* Maximum number of twin systems over lattice structures +integer(pInt), parameter :: lattice_MaxMaxNtwinOfStructure = 12 +!* Slip direction, slip normales and Schmid matrices +real(pReal), dimension(3,3,lattice_MaxMaxNslipOfStructure,lattice_MaxLatticeStructure) :: lattice_Sslip +real(pReal), dimension(6,lattice_MaxMaxNslipOfStructure,lattice_MaxLatticeStructure) :: lattice_Sslip_v +real(pReal), dimension(3,lattice_MaxMaxNslipOfStructure,lattice_MaxLatticeStructure) :: lattice_sn +real(pReal), dimension(3,lattice_MaxMaxNslipOfStructure,lattice_MaxLatticeStructure) :: lattice_sd +real(pReal), dimension(3,lattice_MaxMaxNslipOfStructure,lattice_MaxLatticeStructure) :: lattice_st +!* twin direction, twin normales, Schmid matrices and transformation matrices +real(pReal), dimension(3,3,lattice_MaxMaxNtwinOfStructure,lattice_MaxLatticeStructure) :: lattice_Stwin +real(pReal), dimension(6,lattice_MaxMaxNtwinOfStructure,lattice_MaxLatticeStructure) :: lattice_Stwin_v +real(pReal), dimension(3,lattice_MaxMaxNtwinOfStructure,lattice_MaxLatticeStructure) :: lattice_tn +real(pReal), dimension(3,lattice_MaxMaxNtwinOfStructure,lattice_MaxLatticeStructure) :: lattice_td +real(pReal), dimension(3,lattice_MaxMaxNtwinOfStructure,lattice_MaxLatticeStructure) :: lattice_tt +real(pReal), dimension(3,3,lattice_MaxMaxNtwinOfStructure,lattice_MaxLatticeStructure) :: lattice_Qtwin +real(pReal), dimension(lattice_MaxLatticeStructure), parameter :: lattice_TwinShear = & +reshape((/0.7071067812,0.7071067812,0.7071067812/),(/lattice_MaxLatticeStructure/)) ! Depends surely on c/a ratio for HCP +!* Slip_slip interaction matrices +integer(pInt), dimension(lattice_MaxMaxNslipOfStructure,lattice_MaxMaxNslipOfStructure,lattice_MaxLatticeStructure) :: & +lattice_SlipIntType +!* Twin-twin interaction matrices +integer(pInt), dimension(lattice_MaxMaxNtwinOfStructure,lattice_MaxMaxNtwinOfStructure,lattice_MaxLatticeStructure) :: & +lattice_TwinIntType + +!*** Slip systems for FCC structures (1) *** +!* System {111}<110> Sort according Eisenlohr&Hantcherli +data lattice_sd(:, 1,1)/ 0, 1,-1/ ; data lattice_sn(:, 1,1)/ 1, 1, 1/ +data lattice_sd(:, 2,1)/-1, 0, 1/ ; data lattice_sn(:, 2,1)/ 1, 1, 1/ +data lattice_sd(:, 3,1)/ 1,-1, 0/ ; data lattice_sn(:, 3,1)/ 1, 1, 1/ +data lattice_sd(:, 4,1)/ 0,-1,-1/ ; data lattice_sn(:, 4,1)/-1,-1, 1/ +data lattice_sd(:, 5,1)/ 1, 0, 1/ ; data lattice_sn(:, 5,1)/-1,-1, 1/ +data lattice_sd(:, 6,1)/-1, 1, 0/ ; data lattice_sn(:, 6,1)/-1,-1, 1/ +data lattice_sd(:, 7,1)/ 0,-1, 1/ ; data lattice_sn(:, 7,1)/ 1,-1,-1/ +data lattice_sd(:, 8,1)/-1, 0,-1/ ; data lattice_sn(:, 8,1)/ 1,-1,-1/ +data lattice_sd(:, 9,1)/ 1, 1, 0/ ; data lattice_sn(:, 9,1)/ 1,-1,-1/ +data lattice_sd(:,10,1)/ 0, 1, 1/ ; data lattice_sn(:,10,1)/-1, 1,-1/ +data lattice_sd(:,11,1)/ 1, 0,-1/ ; data lattice_sn(:,11,1)/-1, 1,-1/ +data lattice_sd(:,12,1)/-1,-1, 0/ ; data lattice_sn(:,12,1)/-1, 1,-1/ + +!*** Twin systems for FCC structures (1) *** +!* System {111}<112> Sort according Eisenlohr&Hantcherli +data lattice_td(:, 1,1)/-2, 1, 1/ ; data lattice_tn(:, 1,1)/ 1, 1, 1/ +data lattice_td(:, 2,1)/ 1,-2, 1/ ; data lattice_tn(:, 2,1)/ 1, 1, 1/ +data lattice_td(:, 3,1)/ 1, 1,-2/ ; data lattice_tn(:, 3,1)/ 1, 1, 1/ +data lattice_td(:, 4,1)/ 2,-1, 1/ ; data lattice_tn(:, 4,1)/-1,-1, 1/ +data lattice_td(:, 5,1)/-1, 2, 1/ ; data lattice_tn(:, 5,1)/-1,-1, 1/ +data lattice_td(:, 6,1)/-1,-1,-2/ ; data lattice_tn(:, 6,1)/-1,-1, 1/ +data lattice_td(:, 7,1)/-2,-1,-1/ ; data lattice_tn(:, 7,1)/ 1,-1,-1/ +data lattice_td(:, 8,1)/ 1, 2,-1/ ; data lattice_tn(:, 8,1)/ 1,-1,-1/ +data lattice_td(:, 9,1)/ 1,-1, 2/ ; data lattice_tn(:, 9,1)/ 1,-1,-1/ +data lattice_td(:,10,1)/ 2, 1,-1/ ; data lattice_tn(:,10,1)/-1, 1,-1/ +data lattice_td(:,11,1)/-1,-2,-1/ ; data lattice_tn(:,11,1)/-1, 1,-1/ +data lattice_td(:,12,1)/-1, 1, 2/ ; data lattice_tn(:,12,1)/-1, 1,-1/ + +!*** Slip-Slip interactions for FCC structures (1) *** +data lattice_SlipIntType( 1,1:lattice_MaxNslipOfStructure(1),1)/1,2,2,4,6,5,3,5,5,4,5,6/ +data lattice_SlipIntType( 2,1:lattice_MaxNslipOfStructure(1),1)/2,1,2,6,4,5,5,4,6,5,3,5/ +data lattice_SlipIntType( 3,1:lattice_MaxNslipOfStructure(1),1)/2,2,1,5,5,3,5,6,4,6,5,4/ +data lattice_SlipIntType( 4,1:lattice_MaxNslipOfStructure(1),1)/4,6,5,1,2,2,4,5,6,3,5,5/ +data lattice_SlipIntType( 5,1:lattice_MaxNslipOfStructure(1),1)/6,4,5,2,1,2,5,3,5,5,4,6/ +data lattice_SlipIntType( 6,1:lattice_MaxNslipOfStructure(1),1)/5,5,3,2,2,1,6,5,4,5,6,4/ +data lattice_SlipIntType( 7,1:lattice_MaxNslipOfStructure(1),1)/3,5,5,4,5,6,1,2,2,4,6,5/ +data lattice_SlipIntType( 8,1:lattice_MaxNslipOfStructure(1),1)/5,4,6,5,3,5,2,1,2,6,4,5/ +data lattice_SlipIntType( 9,1:lattice_MaxNslipOfStructure(1),1)/5,6,4,6,5,4,2,2,1,5,5,3/ +data lattice_SlipIntType(10,1:lattice_MaxNslipOfStructure(1),1)/4,5,6,3,5,5,4,6,5,1,2,2/ +data lattice_SlipIntType(11,1:lattice_MaxNslipOfStructure(1),1)/5,3,5,5,4,6,6,4,5,2,1,2/ +data lattice_SlipIntType(12,1:lattice_MaxNslipOfStructure(1),1)/6,5,4,5,6,4,5,5,3,2,2,1/ + +!*** Twin-Twin interactions for FCC structures (1) *** +data lattice_TwinIntType( 1,1:lattice_MaxNtwinOfStructure(1),1)/0,0,0,1,1,1,1,1,1,1,1,1/ +data lattice_TwinIntType( 2,1:lattice_MaxNtwinOfStructure(1),1)/0,0,0,1,1,1,1,1,1,1,1,1/ +data lattice_TwinIntType( 3,1:lattice_MaxNtwinOfStructure(1),1)/0,0,0,1,1,1,1,1,1,1,1,1/ +data lattice_TwinIntType( 4,1:lattice_MaxNtwinOfStructure(1),1)/1,1,1,0,0,0,1,1,1,1,1,1/ +data lattice_TwinIntType( 5,1:lattice_MaxNtwinOfStructure(1),1)/1,1,1,0,0,0,1,1,1,1,1,1/ +data lattice_TwinIntType( 6,1:lattice_MaxNtwinOfStructure(1),1)/1,1,1,0,0,0,1,1,1,1,1,1/ +data lattice_TwinIntType( 7,1:lattice_MaxNtwinOfStructure(1),1)/1,1,1,1,1,1,0,0,0,1,1,1/ +data lattice_TwinIntType( 8,1:lattice_MaxNtwinOfStructure(1),1)/1,1,1,1,1,1,0,0,0,1,1,1/ +data lattice_TwinIntType( 9,1:lattice_MaxNtwinOfStructure(1),1)/1,1,1,1,1,1,0,0,0,1,1,1/ +data lattice_TwinIntType(10,1:lattice_MaxNtwinOfStructure(1),1)/1,1,1,1,1,1,1,1,1,0,0,0/ +data lattice_TwinIntType(11,1:lattice_MaxNtwinOfStructure(1),1)/1,1,1,1,1,1,1,1,1,0,0,0/ +data lattice_TwinIntType(12,1:lattice_MaxNtwinOfStructure(1),1)/1,1,1,1,1,1,1,1,1,0,0,0/ + +!*** Slip systems for BCC structures (2) *** +!* System {110}<111> +!* Sort? +data lattice_sd(:, 1,2)/ 1,-1, 1/ ; data lattice_sn(:, 1,2)/ 0, 1, 1/ +data lattice_sd(:, 2,2)/-1,-1, 1/ ; data lattice_sn(:, 2,2)/ 0, 1, 1/ +data lattice_sd(:, 3,2)/ 1, 1, 1/ ; data lattice_sn(:, 3,2)/ 0,-1, 1/ +data lattice_sd(:, 4,2)/-1, 1, 1/ ; data lattice_sn(:, 4,2)/ 0,-1, 1/ +data lattice_sd(:, 5,2)/-1, 1, 1/ ; data lattice_sn(:, 5,2)/ 1, 0, 1/ +data lattice_sd(:, 6,2)/-1,-1, 1/ ; data lattice_sn(:, 6,2)/ 1, 0, 1/ +data lattice_sd(:, 7,2)/ 1, 1, 1/ ; data lattice_sn(:, 7,2)/-1, 0, 1/ +data lattice_sd(:, 8,2)/ 1,-1, 1/ ; data lattice_sn(:, 8,2)/-1, 0, 1/ +data lattice_sd(:, 9,2)/-1, 1, 1/ ; data lattice_sn(:, 9,2)/ 1, 1, 0/ +data lattice_sd(:,10,2)/-1, 1,-1/ ; data lattice_sn(:,10,2)/ 1, 1, 0/ +data lattice_sd(:,11,2)/ 1, 1, 1/ ; data lattice_sn(:,11,2)/-1, 1, 0/ +data lattice_sd(:,12,2)/ 1, 1,-1/ ; data lattice_sn(:,12,2)/-1, 1, 0/ +!* System {112}<111> +!* Sort? +data lattice_sd(:,13,2)/-1, 1, 1/ ; data lattice_sn(:,13,2)/ 2, 1, 1/ +data lattice_sd(:,14,2)/ 1, 1, 1/ ; data lattice_sn(:,14,2)/-2, 1, 1/ +data lattice_sd(:,15,2)/ 1, 1,-1/ ; data lattice_sn(:,15,2)/ 2,-1, 1/ +data lattice_sd(:,16,2)/ 1,-1, 1/ ; data lattice_sn(:,16,2)/ 2, 1,-1/ +data lattice_sd(:,17,2)/ 1,-1, 1/ ; data lattice_sn(:,17,2)/ 1, 2, 1/ +data lattice_sd(:,18,2)/ 1, 1,-1/ ; data lattice_sn(:,18,2)/-1, 2, 1/ +data lattice_sd(:,19,2)/ 1, 1, 1/ ; data lattice_sn(:,19,2)/ 1,-2, 1/ +data lattice_sd(:,20,2)/-1, 1, 1/ ; data lattice_sn(:,20,2)/ 1, 2,-1/ +data lattice_sd(:,21,2)/ 1, 1,-1/ ; data lattice_sn(:,21,2)/ 1, 1, 2/ +data lattice_sd(:,22,2)/ 1,-1, 1/ ; data lattice_sn(:,22,2)/-1, 1, 2/ +data lattice_sd(:,23,2)/-1, 1, 1/ ; data lattice_sn(:,23,2)/ 1,-1, 2/ +data lattice_sd(:,24,2)/ 1, 1, 1/ ; data lattice_sn(:,24,2)/ 1, 1,-2/ +!* System {123}<111> +!* Sort? +data lattice_sd(:,25,2)/ 1, 1,-1/ ; data lattice_sn(:,25,2)/ 1, 2, 3/ +data lattice_sd(:,26,2)/ 1,-1, 1/ ; data lattice_sn(:,26,2)/-1, 2, 3/ +data lattice_sd(:,27,2)/-1, 1, 1/ ; data lattice_sn(:,27,2)/ 1,-2, 3/ +data lattice_sd(:,28,2)/ 1, 1, 1/ ; data lattice_sn(:,28,2)/ 1, 2,-3/ +data lattice_sd(:,29,2)/ 1,-1, 1/ ; data lattice_sn(:,29,2)/ 1, 3, 2/ +data lattice_sd(:,30,2)/ 1, 1,-1/ ; data lattice_sn(:,30,2)/-1, 3, 2/ +data lattice_sd(:,31,2)/ 1, 1, 1/ ; data lattice_sn(:,31,2)/ 1,-3, 2/ +data lattice_sd(:,32,2)/-1, 1, 1/ ; data lattice_sn(:,32,2)/ 1, 3,-2/ +data lattice_sd(:,33,2)/ 1, 1,-1/ ; data lattice_sn(:,33,2)/ 2, 1, 3/ +data lattice_sd(:,34,2)/ 1,-1, 1/ ; data lattice_sn(:,34,2)/-2, 1, 3/ +data lattice_sd(:,35,2)/-1, 1, 1/ ; data lattice_sn(:,35,2)/ 2,-1, 3/ +data lattice_sd(:,36,2)/ 1, 1, 1/ ; data lattice_sn(:,36,2)/ 2, 1,-3/ +data lattice_sd(:,37,2)/ 1,-1, 1/ ; data lattice_sn(:,37,2)/ 2, 3, 1/ +data lattice_sd(:,38,2)/ 1, 1,-1/ ; data lattice_sn(:,38,2)/-2, 3, 1/ +data lattice_sd(:,39,2)/ 1, 1, 1/ ; data lattice_sn(:,39,2)/ 2,-3, 1/ +data lattice_sd(:,40,2)/-1, 1, 1/ ; data lattice_sn(:,40,2)/ 2, 3,-1/ +data lattice_sd(:,41,2)/-1, 1, 1/ ; data lattice_sn(:,41,2)/ 3, 1, 2/ +data lattice_sd(:,42,2)/ 1, 1, 1/ ; data lattice_sn(:,42,2)/-3, 1, 2/ +data lattice_sd(:,43,2)/ 1, 1,-1/ ; data lattice_sn(:,43,2)/ 3,-1, 2/ +data lattice_sd(:,44,2)/ 1,-1, 1/ ; data lattice_sn(:,44,2)/ 3, 1,-2/ +data lattice_sd(:,45,2)/-1, 1, 1/ ; data lattice_sn(:,45,2)/ 3, 2, 1/ +data lattice_sd(:,46,2)/ 1, 1, 1/ ; data lattice_sn(:,46,2)/-3, 2, 1/ +data lattice_sd(:,47,2)/ 1, 1,-1/ ; data lattice_sn(:,47,2)/ 3,-2, 1/ +data lattice_sd(:,48,2)/ 1,-1, 1/ ; data lattice_sn(:,48,2)/ 3, 2,-1/ + +!*** Twin systems for BCC structures (2) *** +!* System {112}<111> +!* Sort? +!* MISSING: not implemented yet + +!*** Slip-Slip interactions for BCC structures (2) *** +data lattice_SlipIntType( 1,:,2)/1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/ +data lattice_SlipIntType( 2,:,2)/2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/ +data lattice_SlipIntType( 3,:,2)/2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/ +data lattice_SlipIntType( 4,:,2)/2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/ +data lattice_SlipIntType( 5,:,2)/2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/ +data lattice_SlipIntType( 6,:,2)/2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/ +data lattice_SlipIntType( 7,:,2)/2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/ +data lattice_SlipIntType( 8,:,2)/2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/ +data lattice_SlipIntType( 9,:,2)/2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/ +data lattice_SlipIntType(10,:,2)/2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/ +data lattice_SlipIntType(11,:,2)/2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/ +data lattice_SlipIntType(12,:,2)/2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/ +data lattice_SlipIntType(13,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/ +data lattice_SlipIntType(14,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/ +data lattice_SlipIntType(15,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/ +data lattice_SlipIntType(16,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/ +data lattice_SlipIntType(17,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/ +data lattice_SlipIntType(18,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/ +data lattice_SlipIntType(19,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/ +data lattice_SlipIntType(20,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/ +data lattice_SlipIntType(21,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/ +data lattice_SlipIntType(22,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/ +data lattice_SlipIntType(23,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/ +data lattice_SlipIntType(24,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/ +data lattice_SlipIntType(25,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/ +data lattice_SlipIntType(26,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/ +data lattice_SlipIntType(27,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/ +data lattice_SlipIntType(28,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/ +data lattice_SlipIntType(29,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/ +data lattice_SlipIntType(30,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/ +data lattice_SlipIntType(31,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/ +data lattice_SlipIntType(32,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/ +data lattice_SlipIntType(33,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/ +data lattice_SlipIntType(34,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2/ +data lattice_SlipIntType(35,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2/ +data lattice_SlipIntType(36,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2/ +data lattice_SlipIntType(37,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2/ +data lattice_SlipIntType(38,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2/ +data lattice_SlipIntType(39,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2/ +data lattice_SlipIntType(40,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2/ +data lattice_SlipIntType(41,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2/ +data lattice_SlipIntType(42,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2/ +data lattice_SlipIntType(43,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2/ +data lattice_SlipIntType(44,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2/ +data lattice_SlipIntType(45,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2/ +data lattice_SlipIntType(46,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2/ +data lattice_SlipIntType(47,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2/ +data lattice_SlipIntType(48,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1/ + +!*** Twin-twin interactions for BCC structures (2) *** +! MISSING: not implemented yet + +!*** Slip systems for HCP structures (3) *** +!* Basal systems {0001}<1120> (independent of c/a-ratio) +!* 1- (0 0 0 1)[-2 1 1 0] +!* 2- (0 0 0 1)[ 1 -2 1 0] +!* 3- (0 0 0 1)[ 1 1 -2 0] +!* Plane (hkil)->(hkl) +!* Direction [uvtw]->[(u-t) (v-t) w] +!* Automatical transformation from Bravais to Miller +!* not done for the moment +!* Sort? +data lattice_sd(:, 1,3)/-1, 0, 0/ ; data lattice_sn(:, 1,3)/ 0, 0, 1/ +data lattice_sd(:, 2,3)/ 0,-1, 0/ ; data lattice_sn(:, 2,3)/ 0, 0, 1/ +data lattice_sd(:, 3,3)/ 1, 1, 0/ ; data lattice_sn(:, 3,3)/ 0, 0, 1/ +!* 1st type prismatic systems {1010}<1120> (independent of c/a-ratio) +!* 1- ( 0 1 -1 0)[-2 1 1 0] +!* 2- ( 1 0 -1 0)[ 1 -2 1 0] +!* 3- (-1 1 0 0)[ 1 1 -2 0] +!* Sort? +data lattice_sd(:, 4,3)/-1, 0, 0/ ; data lattice_sn(:, 4,3)/ 0, 1, 0/ +data lattice_sd(:, 5,3)/ 0,-1, 0/ ; data lattice_sn(:, 5,3)/ 1, 0, 0/ +data lattice_sd(:, 6,3)/ 1, 1, 0/ ; data lattice_sn(:, 6,3)/-1, 1, 0/ +!* 1st type 1st order pyramidal systems {1011}<1120> +!* plane normales depend on the c/a-ratio +!* 1- ( 0 -1 1 1)[-2 1 1 0] +!* 2- ( 0 1 -1 1)[-2 1 1 0] +!* 3- (-1 0 1 1)[ 1 -2 1 0] +!* 4- ( 1 0 -1 1)[ 1 -2 1 0] +!* 5- (-1 1 0 1)[ 1 1 -2 0] +!* 6- ( 1 -1 0 1)[ 1 1 -2 0] +!* Sort? +data lattice_sd(:, 7,3)/-1, 0, 0/ ; data lattice_sn(:, 7,3)/ 0,-1, 1/ +data lattice_sd(:, 8,3)/ 0,-1, 0/ ; data lattice_sn(:, 8,3)/ 0, 1, 1/ +data lattice_sd(:, 9,3)/ 1, 1, 0/ ; data lattice_sn(:, 9,3)/-1, 0, 1/ +data lattice_sd(:,10,3)/-1, 0, 0/ ; data lattice_sn(:,10,3)/ 1, 0, 1/ +data lattice_sd(:,11,3)/ 0,-1, 0/ ; data lattice_sn(:,11,3)/-1, 1, 1/ +data lattice_sd(:,12,3)/ 1, 1, 0/ ; data lattice_sn(:,12,3)/ 1,-1, 1/ + +!*** Twin systems for HCP structures (3) *** +!* System {1012}<1011> +!* Sort? +!* MISSING: not implemented yet + +!*** Slip-Slip interactions for HCP structures (3) *** +data lattice_SlipIntType( 1,1:lattice_MaxNslipOfStructure(3),3)/1,2,2,2,2,2,2,2,2,2,2,2/ +data lattice_SlipIntType( 2,1:lattice_MaxNslipOfStructure(3),3)/2,1,2,2,2,2,2,2,2,2,2,2/ +data lattice_SlipIntType( 3,1:lattice_MaxNslipOfStructure(3),3)/2,2,1,2,2,2,2,2,2,2,2,2/ +data lattice_SlipIntType( 4,1:lattice_MaxNslipOfStructure(3),3)/2,2,2,1,2,2,2,2,2,2,2,2/ +data lattice_SlipIntType( 5,1:lattice_MaxNslipOfStructure(3),3)/2,2,2,2,1,2,2,2,2,2,2,2/ +data lattice_SlipIntType( 6,1:lattice_MaxNslipOfStructure(3),3)/2,2,2,2,2,1,2,2,2,2,2,2/ +data lattice_SlipIntType( 7,1:lattice_MaxNslipOfStructure(3),3)/2,2,2,2,2,2,1,2,2,2,2,2/ +data lattice_SlipIntType( 8,1:lattice_MaxNslipOfStructure(3),3)/2,2,2,2,2,2,2,1,2,2,2,2/ +data lattice_SlipIntType( 9,1:lattice_MaxNslipOfStructure(3),3)/2,2,2,2,2,2,2,2,1,2,2,2/ +data lattice_SlipIntType(10,1:lattice_MaxNslipOfStructure(3),3)/2,2,2,2,2,2,2,2,2,1,2,2/ +data lattice_SlipIntType(11,1:lattice_MaxNslipOfStructure(3),3)/2,2,2,2,2,2,2,2,2,2,1,2/ +data lattice_SlipIntType(12,1:lattice_MaxNslipOfStructure(3),3)/2,2,2,2,2,2,2,2,2,2,2,1/ + +!*** Twin-twin interactions for HCP structures (3) *** +! MISSING: not implemented yet + + +CONTAINS +!**************************************** +!* - lattice_Init +!* - lattice_SchmidMatrices +!**************************************** + + +subroutine lattice_init() +!************************************** +!* Module initialization * +!************************************** +call lattice_SchmidMatrices() +end subroutine + + +subroutine lattice_SchmidMatrices() +!************************************** +!* Calculation of Schmid matrices * +!************************************** +use prec, only: pReal,pInt +use math, only: math_I3,nrmMandel,mapMandel +implicit none + +!* Definition of variables +integer(pInt) i,j,k,l +real(pReal) norm_d,norm_t,norm_n + +!* Iteration over the lattice structures +do l=1,lattice_MaxLatticeStructure +!* Iteration over the slip systems + do k=1,lattice_MaxNslipOfStructure(l) +!* Definition of transverse direction st for the frame (sd,st,sn) + lattice_st(1,k,l)=lattice_sn(2,k,l)*lattice_sd(3,k,l)-lattice_sn(3,k,l)*lattice_sd(2,k,l) + lattice_st(2,k,l)=lattice_sn(3,k,l)*lattice_sd(1,k,l)-lattice_sn(1,k,l)*lattice_sd(3,k,l) + lattice_st(3,k,l)=lattice_sn(1,k,l)*lattice_sd(2,k,l)-lattice_sn(2,k,l)*lattice_sd(1,k,l) + norm_d=dsqrt(lattice_sd(1,k,l)**2+lattice_sd(2,k,l)**2+lattice_sd(3,k,l)**2) + norm_t=dsqrt(lattice_st(1,k,l)**2+lattice_st(2,k,l)**2+lattice_st(3,k,l)**2) + norm_n=dsqrt(lattice_sn(1,k,l)**2+lattice_sn(2,k,l)**2+lattice_sn(3,k,l)**2) + lattice_sd(:,k,l)=lattice_sd(:,k,l)/norm_d + lattice_st(:,k,l)=lattice_st(:,k,l)/norm_t + lattice_sn(:,k,l)=lattice_sn(:,k,l)/norm_n +!* Defintion of Schmid matrix + forall (i=1:3,j=1:3) lattice_Sslip(i,j,k,l)=lattice_sd(i,k,l)*lattice_sn(j,k,l) +!* Vectorization of normalized Schmid matrix + forall (i=1:6) lattice_Sslip_v(i,k,l) = nrmMandel(i)/2.0_pReal * & + (lattice_Sslip(mapMandel(1,i),mapMandel(2,i),k,l)+lattice_Sslip(mapMandel(2,i),mapMandel(1,i),k,l)) + enddo + +!* Iteration over the twin systems + do k=1,lattice_MaxNtwinOfStructure(l) +!* Definition of transverse direction tt for the frame (td,tt,tn) + lattice_tt(1,k,l)=lattice_tn(2,k,l)*lattice_td(3,k,l)-lattice_tn(3,k,l)*lattice_td(2,k,l) + lattice_tt(2,k,l)=lattice_tn(3,k,l)*lattice_td(1,k,l)-lattice_tn(1,k,l)*lattice_td(3,k,l) + lattice_tt(3,k,l)=lattice_tn(1,k,l)*lattice_td(2,k,l)-lattice_tn(2,k,l)*lattice_td(1,k,l) + norm_d=dsqrt(lattice_td(1,k,l)**2+lattice_td(2,k,l)**2+lattice_td(3,k,l)**2) + norm_t=dsqrt(lattice_tt(1,k,l)**2+lattice_tt(2,k,l)**2+lattice_tt(3,k,l)**2) + norm_n=dsqrt(lattice_tn(1,k,l)**2+lattice_tn(2,k,l)**2+lattice_tn(3,k,l)**2) + lattice_td(:,k,l)=lattice_td(:,k,l)/norm_d + lattice_tt(:,k,l)=lattice_tt(:,k,l)/norm_t + lattice_tn(:,k,l)=lattice_tn(:,k,l)/norm_n +!* Defintion of Schmid matrix and transformation matrices + lattice_Qtwin(:,:,k,l)=-math_I3 + forall (i=1:3,j=1:3) + lattice_Stwin(i,j,k,l)=lattice_td(i,k,l)*lattice_tn(j,k,l) + lattice_Qtwin(i,j,k,l)=lattice_Qtwin(i,j,k,l)+2*lattice_tn(i,k,l)*lattice_tn(j,k,l) + endforall +!* Vectorization of normalized Schmid matrix + lattice_Stwin_v(1,k,l)=lattice_Stwin(1,1,k,l) + lattice_Stwin_v(2,k,l)=lattice_Stwin(2,2,k,l) + lattice_Stwin_v(3,k,l)=lattice_Stwin(3,3,k,l) + !* be compatible with Mandel notation of Tstar + lattice_Stwin_v(4,k,l)=(lattice_Stwin(1,2,k,l)+lattice_Stwin(2,1,k,l))/dsqrt(2.0_pReal) + lattice_Stwin_v(5,k,l)=(lattice_Stwin(2,3,k,l)+lattice_Stwin(3,2,k,l))/dsqrt(2.0_pReal) + lattice_Stwin_v(6,k,l)=(lattice_Stwin(1,3,k,l)+lattice_Stwin(3,1,k,l))/dsqrt(2.0_pReal) + enddo +enddo + +end subroutine + +END MODULE + + + diff --git a/trunk/math.f90 b/trunk/math.f90 index e38162e46..024de791a 100644 --- a/trunk/math.f90 +++ b/trunk/math.f90 @@ -79,6 +79,7 @@ call random_seed() call get_seed(seed) + seed = 1 call halton_seed_set(seed) call halton_ndim_set(3) @@ -306,7 +307,6 @@ InvA = math_identity2nd(dimen) B = A CALL Gauss(dimen,B,InvA,LogAbsDetA,AnzNegEW,error) - RETURN END SUBROUTINE math_invert @@ -1878,6 +1878,5 @@ endif END FUNCTION - END MODULE math diff --git a/trunk/mpie_cpfem_marc2005r3.f90 b/trunk/mpie_cpfem_marc2005r3.f90 index dfdf7f643..ac6b9c545 100644 --- a/trunk/mpie_cpfem_marc2005r3.f90 +++ b/trunk/mpie_cpfem_marc2005r3.f90 @@ -33,8 +33,9 @@ include "math.f90" include "IO.f90" include "mesh.f90" - include "crystal.f90" + include "lattice.f90" include "constitutive.f90" + include "crystallite.f90" include "CPFEM.f90" ! @@ -74,7 +75,37 @@ ! s stress - should be updated by user ! t state variables (comes in at t=n, must be updated ! to have state variables at t=n+1) -! dt increment of state variables +! dt ! Marc common blocks are in fixed format so they have to be pasted in here +! +! Marc common blocks are in fixed format so they have to be pasted in here +! Beware of changes in newer Marc versions -- these are from 2005r3 +! concom is needed for inc, subinc, ncycle, lovl +! include 'concom' + common/concom/ & + iacous, iasmbl, iautth, ibear, icompl, iconj, icreep, ideva(50), idyn, idynt,& + ielas, ielcma, ielect, iform, ifour, iharm, ihcps, iheat, iheatt, ihresp,& + ijoule, ilem, ilnmom, iloren, inc, incext, incsub, ipass, iplres, ipois,& + ipoist, irpflo, ismall, ismalt, isoil, ispect, ispnow, istore, iswep, ithcrp,& + itherm, iupblg, iupdat, jacflg, jel, jparks, largst, lfond, loadup, loaduq,& + lodcor, lovl, lsub, magnet, ncycle, newtnt, newton, noshr, linear, ivscpl,& + icrpim, iradrt, ipshft, itshr, iangin, iupmdr, iconjf, jincfl, jpermg, jhour,& + isolvr, jritz, jtable, jshell, jdoubl, jform, jcentr, imini, kautth, iautof,& + ibukty, iassum, icnstd, icnstt, kmakmas, imethvp,iradrte,iradrtp, iupdate,iupdatp,& + ncycnt, marmen ,idynme, ihavca, ispf, kmini, imixed, largtt, kdoela, iautofg,& + ipshftp,idntrc, ipore, jtablm, jtablc, isnecma,itrnspo,imsdif, jtrnspo,mcnear,& + imech, imecht, ielcmat, ielectt,magnett, imsdift,noplas, jtabls, jactch, jtablth,& + kgmsto ,jpzo, ifricsh, iremkin,iremfor, ishearp,jspf, machining, jlshell,icompsol,& + iupblgfo,jcondir,nstcrp, nactive,ipassref, nstspnt,ibeart,icheckmpc, noline, icuring,& + ishrink,ioffsflg,isetoff, iharmt, inc_incdat, iautspc,ibrake + +! creeps is needed for timinc (time increment) +! include 'creeps' + common/marc_creeps/ & + cptim,timinc,timinc_p,timinc_s,timincm,timinc_a,timinc_b,creept(33),icptim,icfte,icfst,& + icfeq,icftm,icetem,mcreep,jcreep,icpa,icftmp,icfstr,icfqcp,icfcpm,icrppr,icrcha,icpb,iicpmt,iicpa + + + increment of state variables ! ngens size of stress - strain law ! n element number ! nn integration point number @@ -132,30 +163,34 @@ integer(pInt) computationMode - dimension e(*),de(*),t(*),dt(*),g(*),d(ngens,*),s(*), n(2),coord(ncrd,*),disp(ndeg,*),matus(2),dispt(ndeg,*),ffn(itel,*),& + + + 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) + ! Marc common blocks are in fixed format so they have to be pasted in here ! Beware of changes in newer Marc versions -- these are from 2005r3 ! concom is needed for inc, subinc, ncycle, lovl ! include 'concom' - common/concom/ & - iacous, iasmbl, iautth, ibear, icompl, iconj, icreep, ideva(50), idyn, idynt,& - ielas, ielcma, ielect, iform, ifour, iharm, ihcps, iheat, iheatt, ihresp,& - ijoule, ilem, ilnmom, iloren, inc, incext, incsub, ipass, iplres, ipois,& - ipoist, irpflo, ismall, ismalt, isoil, ispect, ispnow, istore, iswep, ithcrp,& - itherm, iupblg, iupdat, jacflg, jel, jparks, largst, lfond, loadup, loaduq,& - lodcor, lovl, lsub, magnet, ncycle, newtnt, newton, noshr, linear, ivscpl,& - icrpim, iradrt, ipshft, itshr, iangin, iupmdr, iconjf, jincfl, jpermg, jhour,& - isolvr, jritz, jtable, jshell, jdoubl, jform, jcentr, imini, kautth, iautof,& - ibukty, iassum, icnstd, icnstt, kmakmas, imethvp,iradrte,iradrtp, iupdate,iupdatp,& - ncycnt, marmen ,idynme, ihavca, ispf, kmini, imixed, largtt, kdoela, iautofg,& - ipshftp,idntrc, ipore, jtablm, jtablc, isnecma,itrnspo,imsdif, jtrnspo,mcnear,& - imech, imecht, ielcmat, ielectt,magnett, imsdift,noplas, jtabls, jactch, jtablth,& - kgmsto ,jpzo, ifricsh, iremkin,iremfor, ishearp,jspf, machining, jlshell,icompsol,& - iupblgfo,jcondir,nstcrp, nactive,ipassref, nstspnt,ibeart,icheckmpc, noline, icuring,& - ishrink,ioffsflg,isetoff, iharmt, inc_incdat, iautspc,ibrake + common/marc_concom/ & + iacous, iasmbl, iautth, ibear, icompl, iconj, icreep, ideva(50), idyn, idynt,& + ielas, ielcma, ielect, iform, ifour, iharm, ihcps, iheat, iheatt, ihresp,& + ijoule, ilem, ilnmom, iloren, inc, incext, incsub, ipass, iplres, ipois,& + ipoist, irpflo, ismall, ismalt, isoil, ispect, ispnow, istore, iswep, ithcrp,& + itherm, iupblg, iupdat, jacflg, jel, jparks, largst, lfond, loadup, loaduq,& + lodcor, lovl, lsub, magnet, ncycle, newtnt, newton, noshr, linear, ivscpl,& + icrpim, iradrt, ipshft, itshr, iangin, iupmdr, iconjf, jincfl, jpermg, jhour,& + isolvr, jritz, jtable, jshell, jdoubl, jform, jcentr, imini, kautth, iautof,& + ibukty, iassum, icnstd, icnstt, kmakmas, imethvp,iradrte,iradrtp, iupdate,iupdatp,& + ncycnt, marmen ,idynme, ihavca, ispf, kmini, imixed, largtt, kdoela, iautofg,& + ipshftp,idntrc, ipore, jtablm, jtablc, isnecma,itrnspo,imsdif, jtrnspo,mcnear,& + imech, imecht, ielcmat, ielectt,magnett, imsdift,noplas, jtabls, jactch, jtablth,& + kgmsto ,jpzo, ifricsh, iremkin,iremfor, ishearp,jspf, machining, jlshell,icompsol,& + iupblgfo,jcondir,nstcrp, nactive,ipassref, nstspnt,ibeart,icheckmpc, noline, icuring,& + ishrink,ioffsflg,isetoff, ioffsetm,iharmt, inc_incdat,iautspc,ibrake, icbush ,istream_input,& + iprsinp,ivlsinp,ifirst_time,ipin_m,jgnstr_glb, imarc_return,iqvcinp,nqvceid,istpnx,imicro1 ! creeps is needed for timinc (time increment) ! include 'creeps' @@ -167,22 +202,32 @@ if (inc == 0) then cycleCounter = 0 else - if (theInc /= inc .or. theCycle /= ncycle .or. theLovl /= lovl) cycleCounter = cycleCounter+1 + if (theCycle > ncycle) cycleCounter = 0 ! reset counter for each cutback + if (theCycle /= ncycle .or. theLovl /= lovl) cycleCounter = cycleCounter+1 ! ping pong + endif + if (cptim > theTime .or. theInc /= inc) then ! reached convergence + lastIncConverged = .true. + outdatedByNewInc = .true. endif - - if (theInc /= inc) outdatedByNewInc = .true. if (mod(cycleCounter,2) /= 0) computationMode = 4 ! recycle if (mod(cycleCounter,4) == 2) computationMode = 3 ! collect if (mod(cycleCounter,4) == 0) computationMode = 2 ! compute + if (computationMode == 4 .and. ncycle == 0 .and. .not. lastIncConverged) & + computationMode = 6 ! recycle but restore known good consistent tangent + if (computationMode == 4 .and. lastIncConverged) then + computationMode = 5 ! recycle and record former consistent tangent + lastIncConverged = .false. + endif if (computationMode == 2 .and. outdatedByNewInc) then + computationMode = 1 ! compute and age former results outdatedByNewInc = .false. - computationMode = 1 ! compute and age former results endif - theInc = inc - theCycle = ncycle - theLovl = lovl + theTime = cptim ! record current starting time + theInc = inc ! record current increment number + theCycle = ncycle ! record current cycle count + theLovl = lovl ! record current lovl call CPFEM_general(computationMode,ffn,ffn1,t(1),timinc,n(1),nn,s,mod(theCycle,2_pInt*ijaco)==0,d,ngens) @@ -196,7 +241,7 @@ END SUBROUTINE ! -! + SUBROUTINE plotv(v,s,sp,etot,eplas,ecreep,t,m,nn,layer,ndi,nshear,jpltcd) !******************************************************************** ! This routine sets user defined output variables for Marc diff --git a/trunk/mpie_cpfem_marc2007r1.f90 b/trunk/mpie_cpfem_marc2007r1.f90 index 315391131..1381c01f3 100644 --- a/trunk/mpie_cpfem_marc2007r1.f90 +++ b/trunk/mpie_cpfem_marc2007r1.f90 @@ -33,8 +33,9 @@ include "math.f90" include "IO.f90" include "mesh.f90" - include "crystal.f90" + include "lattice.f90" include "constitutive.f90" + include "crystallite.f90" include "CPFEM.f90" ! @@ -135,7 +136,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) @@ -172,22 +172,32 @@ if (inc == 0) then cycleCounter = 0 else - if (theInc /= inc .or. theCycle /= ncycle .or. theLovl /= lovl) cycleCounter = cycleCounter+1 + if (theCycle > ncycle) cycleCounter = 0 ! reset counter for each cutback + if (theCycle /= ncycle .or. theLovl /= lovl) cycleCounter = cycleCounter+1 ! ping pong + endif + if (cptim > theTime .or. theInc /= inc) then ! reached convergence + lastIncConverged = .true. + outdatedByNewInc = .true. endif - - if (theInc /= inc) outdatedByNewInc = .true. if (mod(cycleCounter,2) /= 0) computationMode = 4 ! recycle if (mod(cycleCounter,4) == 2) computationMode = 3 ! collect if (mod(cycleCounter,4) == 0) computationMode = 2 ! compute + if (computationMode == 4 .and. ncycle == 0 .and. .not. lastIncConverged) & + computationMode = 6 ! recycle but restore known good consistent tangent + if (computationMode == 4 .and. lastIncConverged) then + computationMode = 5 ! recycle and record former consistent tangent + lastIncConverged = .false. + endif if (computationMode == 2 .and. outdatedByNewInc) then + computationMode = 1 ! compute and age former results outdatedByNewInc = .false. - computationMode = 1 ! compute and age former results endif - theInc = inc - theCycle = ncycle - theLovl = lovl + theTime = cptim ! record current starting time + theInc = inc ! record current increment number + theCycle = ncycle ! record current cycle count + theLovl = lovl ! record current lovl call CPFEM_general(computationMode,ffn,ffn1,t(1),timinc,n(1),nn,s,mod(theCycle,2_pInt*ijaco)==0,d,ngens) @@ -201,7 +211,7 @@ END SUBROUTINE ! -! + SUBROUTINE plotv(v,s,sp,etot,eplas,ecreep,t,m,nn,layer,ndi,nshear,jpltcd) !******************************************************************** ! This routine sets user defined output variables for Marc diff --git a/trunk/prec.f90 b/trunk/prec.f90 index 8078a828d..417ab9478 100644 --- a/trunk/prec.f90 +++ b/trunk/prec.f90 @@ -16,11 +16,17 @@ integer(pInt), parameter :: ijaco = 1_pInt ! frequency of FEM Jacobi update integer(pInt), parameter :: nCutback = 10_pInt ! cutbacks in time-step integration 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 + real(pReal), parameter :: pert_Fg = 1.0e-6_pReal ! strain perturbation for FEM Jacobi integer(pInt), parameter :: nOuter = 10_pInt ! outer loop limit integer(pInt), parameter :: nInner = 1000_pInt ! inner loop limit - real(pReal), parameter :: reltol_Outer = 1.0e-4_pReal ! relative tolerance in outer loop (state) + real(pReal), parameter :: reltol_Outer = 1.0e-6_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) +! + real(pReal), parameter :: resToler = 1.0e-6_pReal ! relative tolerance of residual in GIA iteration + real(pReal), parameter :: resAbsol = 1.0e+0_pReal ! absolute tolerance of residual in GIA iteration (corresponds to 1 Pa) + real(pReal), parameter :: resBound = 1.0e+2_pReal ! relative maximum value (upper bound) for GIA residual + integer(pInt), parameter :: NRiterMax = 24_pInt ! maximum number of GIA iteration + real(pReal), parameter :: subStepMin = 1.0e-3_pReal ! minimum (relative) size of sub-step allowed during cutback END MODULE prec