From cac45cff96bb3c598e7b1f1febedfd833010f739 Mon Sep 17 00:00:00 2001 From: Christoph Kords Date: Tue, 24 May 2011 15:57:59 +0000 Subject: [PATCH] ip coordinates are now updated every cycle; this is needed for the nonlocal internal stress fields * Marc: node displacements are added to initial node coordinates (mesh_node0) to get current node positions (mesh_node), then ip coordinates are deduced * Abaqus: ip coordinates are directly updated, no update of node coordinates! * Spectral: for the moment no update of either ip or node coordinates! passing only dummy values with initial ip coordinates --- code/CPFEM.f90 | 122 ++++++++++++++++++++++---------- code/DAMASK_abaqus_exp.f | 4 +- code/DAMASK_abaqus_std.f | 2 +- code/DAMASK_marc.f90 | 2 +- code/DAMASK_spectral.f90 | 13 ++-- code/DAMASK_spectral_single.f90 | 13 ++-- code/constitutive_nonlocal.f90 | 6 +- code/mesh.f90 | 109 ++++++++++++++++------------ 8 files changed, 174 insertions(+), 97 deletions(-) diff --git a/code/CPFEM.f90 b/code/CPFEM.f90 index 17e9b8daf..ec963c3c5 100644 --- a/code/CPFEM.f90 +++ b/code/CPFEM.f90 @@ -224,7 +224,7 @@ endsubroutine !*** perform initialization at first call, update variables and *** !*** call the actual material model *** !*********************************************************************** -subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchyStress,& +subroutine CPFEM_general(mode, coords, ffn, ffn1, Temperature, dt, element, IP, cauchyStress,& & jacobian, pstress, dPdF) ! note: cauchyStress = Cauchy stress cs(6) and jacobian = Consistent tangent dcs/dE @@ -270,7 +270,13 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt mesh_NcpElems, & mesh_maxNips, & mesh_element, & - FE_Nips + mesh_node0, & + mesh_node, & + mesh_ipCenterOfGravity, & + mesh_build_subNodeCoords, & + mesh_build_ipVolumes, & + FE_Nips, & + FE_Nnodes use material, only: homogenization_maxNgrains, & microstructure_elemhomo, & material_phase @@ -284,7 +290,8 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt crystallite_dPdF0, & crystallite_dPdF, & crystallite_Tstar0_v, & - crystallite_Tstar_v + crystallite_Tstar_v, & + crystallite_localConstitution use homogenization, only: homogenization_sizeState, & homogenization_state, & homogenization_state0, & @@ -303,11 +310,16 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt implicit none + !*** input variables ***! + integer(pInt), intent(in) :: element, & ! FE element number IP ! FE integration point number real(pReal), intent(inout) :: Temperature ! temperature real(pReal), intent(in) :: dt ! time increment + real(pReal), dimension (3,*), intent(in) :: coords ! MARC: displacements for each node of the current element + ! ABAQUS: coordinates of the current material point (IP) + ! SPECTRAL: coordinates of the current material point (IP) real(pReal), dimension (3,3), intent(in) :: ffn, & ! deformation gradient for t=t0 ffn1 ! deformation gradient for t=t1 integer(pInt), intent(in) :: mode ! computation mode 1: regular computation plus aging of results @@ -317,13 +329,17 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt ! 5: restore tangent from former converged inc ! 6: recycling of former results (MARC speciality) + !*** output variables ***! + real(pReal), dimension(6), intent(out) :: cauchyStress ! stress vector in Mandel notation real(pReal), dimension(6,6), intent(out) :: jacobian ! jacobian in Mandel notation real(pReal), dimension (3,3), intent(out) :: pstress ! Piola-Kirchhoff stress in Matrix notation real(pReal), dimension (3,3,3,3), intent(out) :: dPdF ! - + + !*** local variables ***! + real(pReal) J_inverse, & ! inverse of Jacobian rnd real(pReal), dimension (3,3) :: Kirchhoff, & ! Piola-Kirchhoff stress in Matrix notation @@ -338,17 +354,11 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt l, & m, & n, & - e + e, & + node, & + FEnodeID logical updateJaco ! flag indicating if JAcobian has to be updated - !*** global variables ***! - ! CPFEM_cs, & - ! CPFEM_dcsdE, & - ! CPFEM_dcsdE_knownGood, & - ! CPFEM_init_done, & - ! CPFEM_calc_done, & - ! CPFEM_odd_stress, & - ! CPFEM_odd_jacobian cp_en = mesh_FEasCP('elem',element) @@ -367,12 +377,18 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt !$OMP END CRITICAL (write2out) endif - ! according to our "mode" we decide what to do + + !*** according to our "mode" we decide what to do + select case (mode) + ! --+>> REGULAR COMPUTATION (WITH AGING OF RESULTS IF MODE == 1) <<+-- + case (1,2,8,9) - ! age results if mode == 1 + + !*** age results + if (mode == 1 .or. mode == 8) then crystallite_F0 = crystallite_partionedF ! crystallite deformation (_subF is perturbed...) crystallite_Fp0 = crystallite_Fp ! crystallite plastic deformation @@ -402,8 +418,8 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt enddo !$OMP END PARALLEL DO - - ! *** dump the last converged values of each essential variable to a binary file + ! * dump the last converged values of each essential variable to a binary file + if (restartWrite) then if (debug_verbosity > 0) then !$OMP CRITICAL (write2out) @@ -459,7 +475,7 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt close (777) endif endif - ! *** end of dumping + ! * end of dumping endif if (mode == 8 .or. mode == 9) then ! Abaqus explicit skips collect @@ -468,7 +484,9 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt materialpoint_F(1:3,1:3,IP,cp_en) = ffn1 endif - ! deformation gradient outdated or any actual deformation gradient differs more than relevantStrain from the stored one + + !*** deformation gradient outdated or any actual deformation gradient differs more than relevantStrain from the stored one + if (terminallyIll .or. outdatedFFN1 .or. any(abs(ffn1 - materialpoint_F(1:3,1:3,IP,cp_en)) > defgradTolerance)) then if (.not. terminallyIll .and. .not. outdatedFFN1) then if (debug_verbosity > 0) then @@ -485,14 +503,15 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt CPFEM_cs(1:6,IP,cp_en) = rnd*CPFEM_odd_stress CPFEM_dcsde(1:6,1:6,IP,cp_en) = CPFEM_odd_jacobian*math_identity2nd(6) - ! deformation gradient is not outdated + + !*** deformation gradient is not outdated + else - ! set flag for Jacobian update updateJaco = mod(cycleCounter,iJacoStiffness) == 0 - ! no parallel computation + !* no parallel computation, so we use just one single element and IP for computation + if (.not. parallelExecution) then - ! we just take one single element and IP FEsolving_execElem(1) = cp_en FEsolving_execElem(2) = cp_en FEsolving_execIP(1,cp_en) = IP @@ -505,13 +524,21 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt call materialpoint_stressAndItsTangent(updateJaco, dt) ! calculate stress and its tangent call materialpoint_postResults(dt) ! post results - ! parallel computation and calulation not yet done + !* parallel computation and calulation not yet done + elseif (.not. CPFEM_calc_done) then if (debug_verbosity > 0) then !$OMP CRITICAL (write2out) write(6,'(a,i8,a,i8)') '<< CPFEM >> Calculation for elements ',FEsolving_execElem(1),' to ',FEsolving_execElem(2) !$OMP END CRITICAL (write2out) endif + if (any(.not. crystallite_localConstitution) .and. FEsolver == 'Marc') then + call mesh_build_subNodeCoords() ! update subnodal coordinates + call mesh_build_ipVolumes() ! update ip center of gravity + endif + !$OMP CRITICAL (write2out) + write(6,'(a,i5,a,i5)') '<< CPFEM >> Start stress and tangent ',FEsolving_execElem(1),' to ',FEsolving_execElem(2) + !$OMP END CRITICAL (write2out) call materialpoint_stressAndItsTangent(updateJaco, dt) ! calculate stress and its tangent (parallel execution inside) call materialpoint_postResults(dt) ! post results !$OMP PARALLEL DO @@ -529,6 +556,9 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt CPFEM_calc_done = .true. endif + + !*** map stress and stiffness (or return odd values if terminally ill) + if ( terminallyIll ) then call random_number(rnd) rnd = 2.0_pReal * rnd - 1.0_pReal @@ -558,7 +588,9 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt endif endif + ! --+>> COLLECTION OF FEM INPUT WITH RETURNING OF RANDOMIZED ODD STRESS AND JACOBIAN <<+-- + case (3,4,5) if (mode == 4) then CPFEM_dcsde_knownGood = CPFEM_dcsde ! --+>> BACKUP JACOBIAN FROM FORMER CONVERGED INC @@ -573,29 +605,41 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt CPFEM_cs(1:6,IP,cp_en) = rnd * CPFEM_odd_stress CPFEM_dcsde(1:6,1:6,IP,cp_en) = CPFEM_odd_jacobian * math_identity2nd(6) CPFEM_calc_done = .false. + select case (FEsolver) + case ('Abaqus','Spectral') + mesh_ipCenterOfGravity(1:3,IP,cp_en) = coords(1:3,1) + case ('Marc') + do node = 1,FE_Nnodes(mesh_element(2,cp_en)) + FEnodeID = mesh_FEasCP('node',mesh_element(4+node,cp_en)) + mesh_node(1:3,FEnodeID) = mesh_node0(1:3,FEnodeID) + coords(1:3,node) + enddo + end select + ! --+>> RECYCLING OF FORMER RESULTS (MARC SPECIALTY) <<+-- + case (6) ! do nothing + + ! --+>> RESTORE CONSISTENT JACOBIAN FROM FORMER CONVERGED INC + case (7) CPFEM_dcsde = CPFEM_dcsde_knownGood end select - ! return the local stress and the jacobian from storage + + !*** fill output variables with computed values + cauchyStress = CPFEM_cs(1:6,IP,cp_en) jacobian = CPFEM_dcsdE(1:6,1:6,IP,cp_en) - - ! copy P and dPdF to the output variables pstress = materialpoint_P(1:3,1:3,IP,cp_en) dPdF = materialpoint_dPdF(1:3,1:3,1:3,1:3,IP,cp_en) - - ! warning for zero stiffness - if (all(abs(jacobian) < 1e-10_pReal)) then - call IO_warning(601,cp_en,IP) + if (theTime > 0.0_pReal) then + Temperature = materialpoint_Temperature(IP,cp_en) ! homogenized result except for potentially non-isothermal starting condition. endif - + if (mode < 6 .and. debug_verbosity > 0 .and. ((debug_e == cp_en .and. debug_i == IP) .or. .not. debug_selectiveDebugger)) then !$OMP CRITICAL (write2out) write(6,'(a,i8,x,i2,/,12(x),6(f10.3,x)/)') '<< CPFEM >> stress/MPa at el ip ', cp_en, IP, cauchyStress/1e6 @@ -604,7 +648,16 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt !$OMP END CRITICAL (write2out) endif - ! remember extreme values of stress and jacobian + + !*** warn if stiffness close to zero + + if (all(abs(jacobian) < 1e-10_pReal)) then + call IO_warning(601,cp_en,IP) + endif + + + !*** remember extreme values of stress and jacobian + if (mode < 3) then cauchyStress33 = math_Mandel6to33(cauchyStress) if (maxval(cauchyStress33) > debug_stressMax) then @@ -626,9 +679,6 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt endif endif - ! return temperature - if (theTime > 0.0_pReal) Temperature = materialpoint_Temperature(IP,cp_en) ! homogenized result except for potentially non-isothermal starting condition. - end subroutine END MODULE CPFEM diff --git a/code/DAMASK_abaqus_exp.f b/code/DAMASK_abaqus_exp.f index 2a17a9bcc..032582465 100644 --- a/code/DAMASK_abaqus_exp.f +++ b/code/DAMASK_abaqus_exp.f @@ -193,7 +193,7 @@ subroutine vumat (jblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal, & stateOld(nblock,nstatev), enerInternOld(nblock), & enerInelasOld(nblock), tempNew(nblock), tempOld(nblock), & stretchNew(nblock,ndir+nshr), defgradNew(nblock,ndir+nshr+nshr), & - stressNew(nblock,ndir+nshr) + stressNew(nblock,ndir+nshr), coordMp(nblock,3) dimension enerInelasNew(nblock),stateNew(nblock,nstatev),enerInternNew(nblock) dimension nElement(nblock),nMatPoint(nblock) @@ -283,7 +283,7 @@ subroutine vumat (jblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal, & defgrd1(3,2) = defgradNew(n,8) endif - call CPFEM_general(computationMode,defgrd0,defgrd1,temp,timeInc,nElement(n),nMatPoint(n),stress,ddsdde, pstress, dPdF) + call CPFEM_general(computationMode,coordMp(n,1:3),defgrd0,defgrd1,temp,timeInc,nElement(n),nMatPoint(n),stress,ddsdde, pstress, dPdF) ! Mandel: 11, 22, 33, SQRT(2)*12, SQRT(2)*23, SQRT(2)*13 ! straight: 11, 22, 33, 12, 23, 13 diff --git a/code/DAMASK_abaqus_std.f b/code/DAMASK_abaqus_std.f index 172a6dc42..c542fe819 100644 --- a/code/DAMASK_abaqus_std.f +++ b/code/DAMASK_abaqus_std.f @@ -261,7 +261,7 @@ subroutine UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,& !$OMP END CRITICAL (write2out) endif - call CPFEM_general(computationMode,dfgrd0,dfgrd1,temp,dtime,noel,npt,stress_h,ddsdde_h, pstress, dPdF) + call CPFEM_general(computationMode,COORDS,dfgrd0,dfgrd1,temp,dtime,noel,npt,stress_h,ddsdde_h, pstress, dPdF) ! Mandel: 11, 22, 33, SQRT(2)*12, SQRT(2)*23, SQRT(2)*13 ! straight: 11, 22, 33, 12, 23, 13 diff --git a/code/DAMASK_marc.f90 b/code/DAMASK_marc.f90 index 7d046d5da..43e74b1af 100644 --- a/code/DAMASK_marc.f90 +++ b/code/DAMASK_marc.f90 @@ -353,7 +353,7 @@ subroutine hypela2(& lastMode = calcMode(nn,cp_en) ! record calculationMode endif - call CPFEM_general(computationMode,ffn,ffn1,t(1),timinc,n(1),nn,stress,ddsdde, pstress, dPdF) + call CPFEM_general(computationMode,dispt,ffn,ffn1,t(1),timinc,n(1),nn,stress,ddsdde, pstress, dPdF) ! Mandel: 11, 22, 33, SQRT(2)*12, SQRT(2)*23, SQRT(2)*13 ! Marc: 11, 22, 33, 12, 23, 13 diff --git a/code/DAMASK_spectral.f90 b/code/DAMASK_spectral.f90 index db62faf24..a36e0b05d 100644 --- a/code/DAMASK_spectral.f90 +++ b/code/DAMASK_spectral.f90 @@ -91,6 +91,7 @@ program DAMASK_spectral real(pReal), dimension(6) :: cstress ! cauchy stress in Mandel notation real(pReal), dimension(6,6) :: dsde, c066, s066 ! Mandel notation of 4th order tensors real(pReal), dimension(:,:,:,:,:), allocatable :: workfft, defgrad, defgradold + real(pReal), dimension(:,:,:,:), allocatable :: coordinates ! variables storing information for spectral method complex(pReal) :: img @@ -261,6 +262,7 @@ program DAMASK_spectral allocate (defgrad (resolution(1),resolution(2),resolution(3),3,3)); defgrad = 0.0_pReal allocate (defgradold(resolution(1),resolution(2),resolution(3),3,3)); defgradold = 0.0_pReal + allocate (coordinates(3,resolution(1),resolution(2),resolution(3))); coordinates = 0.0_pReal wgt = 1.0_pReal/real(resolution(1)*resolution(2)*resolution(3), pReal) defgradAim = math_I3 @@ -275,7 +277,8 @@ program DAMASK_spectral defgradold(i,j,k,:,:) = math_I3 ! no deformation at the beginning defgrad(i,j,k,:,:) = math_I3 ielem = ielem +1 - call CPFEM_general(2,math_I3,math_I3,temperature,0.0_pReal,ielem,1_pInt,cstress,dsde,pstress,dPdF) + coordinates(1:3,i,j,k) = mesh_ipCenterOfGravity(1:3,1,ielem) ! set to initial coordinates (SHOULD BE UPDATED TO CURRENT POSITION IN FUTURE REVISIONS!!!) + call CPFEM_general(2,coordinates(1:3,i,j,k),math_I3,math_I3,temperature,0.0_pReal,ielem,1_pInt,cstress,dsde,pstress,dPdF) c066 = c066 + dsde enddo; enddo; enddo c066 = c066 * wgt @@ -410,7 +413,7 @@ program DAMASK_spectral ielem = 0_pInt do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1) ielem = ielem + 1 - call CPFEM_general(3, defgradold(i,j,k,:,:), defgrad(i,j,k,:,:),& + call CPFEM_general(3, coordinates(1:3,i,j,k), defgradold(i,j,k,:,:), defgrad(i,j,k,:,:),& temperature,timeinc,ielem,1_pInt,& cstress,dsde, pstress, dPdF) enddo; enddo; enddo @@ -419,6 +422,7 @@ program DAMASK_spectral do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1) ielem = ielem + 1_pInt call CPFEM_general(CPFEM_mode,& ! first element in first iteration retains CPFEM_mode 1, + coordinates(1:3,i,j,k), defgradold(i,j,k,:,:), defgrad(i,j,k,:,:),& ! others get 2 (saves winding forward effort) temperature,timeinc,ielem,1_pInt,& cstress,dsde, pstress, dPdF) @@ -469,7 +473,7 @@ program DAMASK_spectral ielem = 0_pInt do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1) ielem = ielem + 1_pInt - call CPFEM_general(3, defgradold(i,j,k,:,:), defgrad(i,j,k,:,:),& + call CPFEM_general(3, coordinates(1:3,i,j,k), defgradold(i,j,k,:,:), defgrad(i,j,k,:,:),& temperature,timeinc,ielem,1_pInt,& cstress,dsde, pstress, dPdF) enddo; enddo; enddo @@ -477,6 +481,7 @@ program DAMASK_spectral do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1) ielem = ielem + 1_pInt call CPFEM_general(CPFEM_mode,& ! first element in first iteration retains CPFEM_mode 1, + coordinates(1:3,i,j,k), defgradold(i,j,k,:,:), defgrad(i,j,k,:,:),& temperature,timeinc,ielem,1_pInt,& cstress,dsde, pstress, dPdF) @@ -594,4 +599,4 @@ subroutine quit(id) integer(pInt) id stop -end subroutine \ No newline at end of file +end subroutine diff --git a/code/DAMASK_spectral_single.f90 b/code/DAMASK_spectral_single.f90 index b4e6be4f8..9adb0914d 100644 --- a/code/DAMASK_spectral_single.f90 +++ b/code/DAMASK_spectral_single.f90 @@ -91,6 +91,7 @@ program DAMASK_spectral real(pReal), dimension(6) :: cstress ! cauchy stress in Mandel notation real(pReal), dimension(6,6) :: dsde, c066, s066 ! Mandel notation of 4th order tensors real(pReal), dimension(:,:,:,:,:), allocatable :: workfft, defgrad, defgradold + real(pReal), dimension(:,:,:,:), allocatable :: coordinates ! variables storing information for spectral method complex(pReal) :: img @@ -261,6 +262,7 @@ program DAMASK_spectral allocate (defgrad (resolution(1),resolution(2),resolution(3),3,3)); defgrad = 0.0_pReal allocate (defgradold(resolution(1),resolution(2),resolution(3),3,3)); defgradold = 0.0_pReal + allocate (coordinates(3,resolution(1),resolution(2),resolution(3))); coordinates = 0.0_pReal wgt = 1.0_pReal/real(resolution(1)*resolution(2)*resolution(3), pReal) defgradAim = math_I3 @@ -275,7 +277,8 @@ program DAMASK_spectral defgradold(i,j,k,:,:) = math_I3 !no deformation at the beginning defgrad(i,j,k,:,:) = math_I3 ielem = ielem +1 - call CPFEM_general(2,math_I3,math_I3,temperature,0.0_pReal,ielem,1_pInt,cstress,dsde,pstress,dPdF) + coordinates(1:3,i,j,k) = mesh_ipCenterOfGravity(1:3,1,ielem) ! set to initial coordinates (SHOULD BE UPDATED TO CURRENT POSITION IN FUTURE REVISIONS!!!) + call CPFEM_general(2,coordinates(1:3,i,j,k),math_I3,math_I3,temperature,0.0_pReal,ielem,1_pInt,cstress,dsde,pstress,dPdF) c066 = c066 + dsde enddo; enddo; enddo c066 = c066 * wgt @@ -411,7 +414,7 @@ program DAMASK_spectral ielem = 0_pInt do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1) ielem = ielem + 1 - call CPFEM_general(3, defgradold(i,j,k,:,:), defgrad(i,j,k,:,:),& + call CPFEM_general(3, coordinates(1:3,i,j,k),defgradold(i,j,k,:,:), defgrad(i,j,k,:,:),& temperature,timeinc,ielem,1_pInt,& cstress,dsde, pstress, dPdF) enddo; enddo; enddo @@ -420,6 +423,7 @@ program DAMASK_spectral do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1) ielem = ielem + 1_pInt call CPFEM_general(CPFEM_mode,& ! first element in first iteration retains CPFEM_mode 1, + coordinates(1:3,i,j,k), defgradold(i,j,k,:,:), defgrad(i,j,k,:,:),& ! others get 2 (saves winding forward effort) temperature,timeinc,ielem,1_pInt,& cstress,dsde, pstress, dPdF) @@ -470,7 +474,7 @@ program DAMASK_spectral ielem = 0_pInt do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1) ielem = ielem + 1_pInt - call CPFEM_general(3, defgradold(i,j,k,:,:), defgrad(i,j,k,:,:),& + call CPFEM_general(3, coordinates(1:3,i,j,k),defgradold(i,j,k,:,:), defgrad(i,j,k,:,:),& temperature,timeinc,ielem,1_pInt,& cstress,dsde, pstress, dPdF) enddo; enddo; enddo @@ -478,6 +482,7 @@ program DAMASK_spectral do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1) ielem = ielem + 1_pInt call CPFEM_general(CPFEM_mode,& ! first element in first iteration retains CPFEM_mode 1, + coordinates(1:3,i,j,k), defgradold(i,j,k,:,:), defgrad(i,j,k,:,:),& temperature,timeinc,ielem,1_pInt,& cstress,dsde, pstress, dPdF) @@ -595,4 +600,4 @@ subroutine quit(id) integer(pInt) id stop -end subroutine \ No newline at end of file +end subroutine diff --git a/code/constitutive_nonlocal.f90 b/code/constitutive_nonlocal.f90 index c3be20d76..48459ce3c 100644 --- a/code/constitutive_nonlocal.f90 +++ b/code/constitutive_nonlocal.f90 @@ -837,7 +837,7 @@ use debug, only: debug_verbosity, & use mesh, only: mesh_NcpElems, & mesh_maxNips, & mesh_element, & - mesh_node, & + mesh_node0, & FE_Nips, & mesh_ipCenterOfGravity, & mesh_ipVolume, & @@ -968,8 +968,8 @@ if (.not. phase_localConstitution(phase)) then !* in case of periodic surfaces we have to find out how many periodic images in each direction we need do dir = 1,3 - maxCoord(dir) = maxval(mesh_node(dir,:)) - minCoord(dir) = minval(mesh_node(dir,:)) + maxCoord(dir) = maxval(mesh_node0(dir,:)) + minCoord(dir) = minval(mesh_node0(dir,:)) enddo meshSize = maxCoord - minCoord ipCoords = mesh_ipCenterOfGravity(1:3,ip,el) diff --git a/code/mesh.f90 b/code/mesh.f90 index ac108a2e8..ff8aa83ac 100644 --- a/code/mesh.f90 +++ b/code/mesh.f90 @@ -35,7 +35,8 @@ ! _maxNsharedElems : max number of CP elements sharing a node ! ! _element : FEid, type(internal representation), material, texture, node indices -! _node : x,y,z coordinates (initially!) +! _node0 : x,y,z coordinates (initially!) +! _node : x,y,z coordinates (after deformation!) ! _sharedElem : entryCount and list of elements containing node ! ! _mapFEtoCPelem : [sorted FEid, corresponding CPid] @@ -56,7 +57,6 @@ ! _nodeOnFace : list of node indices on each face of a specific type of element ! _maxNnodesAtIP : max number of (equivalent) nodes attached to an IP ! _nodesAtIP : map IP index to two node indices in a specific type of element -! _ipNeighborhood : 6 or less neighboring IPs as [element_num, IP_index] ! _NsubNodes : # subnodes required to fully define all IP volumes ! order is +x,-x,+y,-y,+z,-z but meaning strongly depends on Elemtype @@ -70,19 +70,20 @@ mesh_mapMaterial ! name of elementSet for material integer(pInt), dimension(:,:), allocatable :: mesh_mapElemSet ! list of elements in elementSet integer(pInt), dimension(:,:), allocatable, target :: mesh_mapFEtoCPelem, mesh_mapFEtoCPnode - integer(pInt), dimension(:,:), allocatable :: mesh_element, & - mesh_sharedElem, & + integer(pInt), dimension(:,:), allocatable :: mesh_element, & ! FEid, type(internal representation), material, texture, node indices + mesh_sharedElem, & ! entryCount and list of elements containing node mesh_nodeTwins ! node twins are surface nodes that lie exactly on opposite sides of the mesh (surfaces nodes with equal coordinate values in two dimensions) - integer(pInt), dimension(:,:,:,:), allocatable :: mesh_ipNeighborhood + integer(pInt), dimension(:,:,:,:), allocatable :: mesh_ipNeighborhood ! 6 or less neighboring IPs as [element_num, IP_index] real(pReal), dimension(:,:,:), allocatable :: mesh_subNodeCoord ! coordinates of subnodes per element - real(pReal), dimension(:,:), allocatable :: mesh_node, & - mesh_ipVolume ! volume associated with IP - real(pReal), dimension(:,:,:), allocatable :: mesh_ipArea, & ! area of interface to neighboring IP - mesh_ipCenterOfGravity ! center of gravity of IP - real(pReal), dimension(:,:,:,:), allocatable :: mesh_ipAreaNormal ! area normal of interface to neighboring IP + real(pReal), dimension(:,:), allocatable :: mesh_node0, & ! node coordinates (initially!) + mesh_node, & ! node coordinates (after deformation! ONLY FOR MARC!!!) + mesh_ipVolume ! volume associated with IP (initially!) + real(pReal), dimension(:,:,:), allocatable :: mesh_ipArea, & ! area of interface to neighboring IP (initially!) + mesh_ipCenterOfGravity ! center of gravity of IP (after deformation!) + real(pReal), dimension(:,:,:,:), allocatable :: mesh_ipAreaNormal ! area normal of interface to neighboring IP (initially!) - integer(pInt), dimension(:,:,:), allocatable :: FE_nodesAtIP + integer(pInt), dimension(:,:,:), allocatable :: FE_nodesAtIP ! map IP index to two node indices in a specific type of element integer(pInt), dimension(:,:,:), allocatable :: FE_ipNeighbor integer(pInt), dimension(:,:,:), allocatable :: FE_subNodeParent integer(pInt), dimension(:,:,:,:), allocatable :: FE_subNodeOnIPFace @@ -776,7 +777,7 @@ FE_ipNeighbor(:FE_NipNeighbors(8),:FE_Nips(8),8) = & ! element 117 /),(/FE_NipNeighbors(10),FE_Nips(10)/)) ! *** FE_subNodeParent *** - ! lists the group of IPs for which the center of gravity + ! lists the group of nodes for which the center of gravity ! corresponds to the location of a each subnode. ! example: face-centered subnode with faceNodes 1,2,3,4 to be used in, ! e.g., a 8 IP grid, would be encoded: @@ -2417,6 +2418,7 @@ subroutine mesh_marc_count_cpSizes (unit) character(len=64) tag character(len=1024) line + allocate ( mesh_node0 (3,mesh_Nnodes) ); mesh_node0 = 0_pInt allocate ( mesh_node (3,mesh_Nnodes) ); mesh_node = 0_pInt a = 1_pInt @@ -2473,10 +2475,12 @@ subroutine mesh_marc_count_cpSizes (unit) if (x <= 0.0_pReal .or. y <= 0.0_pReal .or. z <= 0.0_pReal) call IO_error(44) forall (n = 0:mesh_Nnodes-1) - mesh_node(1,n+1) = x * dble(mod(n,a) / (a-1.0_pReal)) - mesh_node(2,n+1) = y * dble(mod(n/a,b) / (b-1.0_pReal)) - mesh_node(3,n+1) = z * dble(mod(n/a/b,c) / (c-1.0_pReal)) + mesh_node0(1,n+1) = x * dble(mod(n,a) / (a-1.0_pReal)) + mesh_node0(2,n+1) = y * dble(mod(n/a,b) / (b-1.0_pReal)) + mesh_node0(3,n+1) = z * dble(mod(n/a/b,c) / (c-1.0_pReal)) end forall + + mesh_node = mesh_node0 100 return @@ -2502,6 +2506,7 @@ subroutine mesh_marc_count_cpSizes (unit) integer(pInt) unit,i,j,m + allocate ( mesh_node0 (3,mesh_Nnodes) ); mesh_node0 = 0_pInt allocate ( mesh_node (3,mesh_Nnodes) ); mesh_node = 0_pInt 610 FORMAT(A300) @@ -2515,13 +2520,14 @@ subroutine mesh_marc_count_cpSizes (unit) do i=1,mesh_Nnodes read (unit,610,END=670) line m = mesh_FEasCP('node',IO_fixedIntValue(line,node_ends,1)) - forall (j = 1:3) mesh_node(j,m) = IO_fixedNoEFloatValue (line,node_ends,j+1) + forall (j = 1:3) mesh_node0(j,m) = IO_fixedNoEFloatValue(line,node_ends,j+1) enddo exit endif enddo -670 return +670 mesh_node = mesh_node0 +return endsubroutine @@ -2545,6 +2551,7 @@ subroutine mesh_marc_count_cpSizes (unit) integer(pInt) unit,i,j,m,count logical inPart + allocate ( mesh_node0 (3,mesh_Nnodes) ); mesh_node0 = 0_pInt allocate ( mesh_node (3,mesh_Nnodes) ); mesh_node = 0_pInt 610 FORMAT(A300) @@ -2573,12 +2580,13 @@ subroutine mesh_marc_count_cpSizes (unit) read (unit,610,END=670) line pos = IO_stringPos(line,maxNchunks) m = mesh_FEasCP('node',IO_intValue(line,pos,1)) - forall (j=1:3) mesh_node(j,m) = IO_floatValue(line,pos,j+1) + forall (j=1:3) mesh_node0(j,m) = IO_floatValue(line,pos,j+1) enddo endif enddo -670 if (size(mesh_node,2) /= mesh_Nnodes) call IO_error(909) +670 if (size(mesh_node0,2) /= mesh_Nnodes) call IO_error(909) + mesh_node = mesh_node0 return endsubroutine @@ -3098,25 +3106,26 @@ endsubroutine integer(pInt) e,t,n,p - allocate(mesh_subNodeCoord(3,mesh_maxNnodes+mesh_maxNsubNodes,mesh_NcpElems)) ; mesh_subNodeCoord = 0.0_pReal + if (.not. allocated(mesh_subNodeCoord)) then + allocate(mesh_subNodeCoord(3,mesh_maxNnodes+mesh_maxNsubNodes,mesh_NcpElems)) + endif + mesh_subNodeCoord = 0.0_pReal do e = 1,mesh_NcpElems ! loop over cpElems t = mesh_element(2,e) ! get elemType do n = 1,FE_Nnodes(t) - mesh_subNodeCoord(:,n,e) = mesh_node(:,mesh_FEasCP('node',mesh_element(4+n,e))) ! loop over nodes of this element type + mesh_subNodeCoord(1:3,n,e) = mesh_node(1:3,mesh_FEasCP('node',mesh_element(4+n,e))) ! loop over nodes of this element type enddo do n = 1,FE_NsubNodes(t) ! now for the true subnodes do p = 1,FE_Nips(t) ! loop through possible parent nodes if (FE_subNodeParent(p,n,t) > 0) & ! valid parent node - mesh_subNodeCoord(:,n+FE_Nnodes(t),e) = & - mesh_subNodeCoord(:,n+FE_Nnodes(t),e) + & - mesh_node(:,mesh_FEasCP('node',mesh_element(4+FE_subNodeParent(p,n,t),e))) ! add up parents + mesh_subNodeCoord(1:3,FE_Nnodes(t)+n,e) = mesh_subNodeCoord(1:3,FE_Nnodes(t)+n,e) & + + mesh_node(1:3,mesh_FEasCP('node',mesh_element(4+FE_subNodeParent(p,n,t),e))) ! add up parents enddo - mesh_subNodeCoord(:,n+FE_Nnodes(t),e) = mesh_subNodeCoord(:,n+FE_Nnodes(t),e) / count(FE_subNodeParent(:,n,t) > 0) + mesh_subNodeCoord(1:3,n+FE_Nnodes(t),e) = mesh_subNodeCoord(1:3,n+FE_Nnodes(t),e) & + / count(FE_subNodeParent(:,n,t) > 0) enddo enddo - - return endsubroutine @@ -3140,9 +3149,15 @@ endsubroutine real(pReal), dimension(3,FE_NipFaceNodes) :: nPos ! coordinates of nodes on IP face real(pReal), dimension(Ntriangles,FE_NipFaceNodes) :: volume ! volumes of possible tetrahedra real(pReal), dimension(3) :: centerOfGravity + logical :: calcIPvolume = .false. - allocate(mesh_ipVolume(mesh_maxNips,mesh_NcpElems)) ; mesh_ipVolume = 0.0_pReal - allocate(mesh_ipCenterOfGravity(3,mesh_maxNips,mesh_NcpElems)) ; mesh_ipCenterOfGravity = 0.0_pReal + if (.not. allocated(mesh_ipVolume)) then + allocate(mesh_ipVolume(mesh_maxNips,mesh_NcpElems)) + allocate(mesh_ipCenterOfGravity(3,mesh_maxNips,mesh_NcpElems)) + mesh_ipVolume = 0.0_pReal + mesh_ipCenterOfGravity = 0.0_pReal + calcIPvolume = .true. + endif do e = 1,mesh_NcpElems ! loop over cpElems t = mesh_element(2,e) ! get elemType @@ -3168,21 +3183,23 @@ endsubroutine endif enddo centerOfGravity = sum(gravityNodePos,2)/count(gravityNode) - - do f = 1,FE_NipNeighbors(t) ! loop over interfaces of IP and add tetrahedra which connect to CoG - forall (n = 1:FE_NipFaceNodes) nPos(:,n) = mesh_subNodeCoord(:,FE_subNodeOnIPFace(n,f,i,t),e) - forall (n = 1:FE_NipFaceNodes, j = 1:Ntriangles) & ! start at each interface node and build valid triangles to cover interface - volume(j,n) = math_volTetrahedron(nPos(:,n), & ! calc volume of respective tetrahedron to CoG - nPos(:,1+mod(n-1 +j ,FE_NipFaceNodes)), & ! start at offset j - nPos(:,1+mod(n-1 +j+1,FE_NipFaceNodes)), & ! and take j's neighbor - centerOfGravity) - mesh_ipVolume(i,e) = mesh_ipVolume(i,e) + sum(volume) ! add contribution from this interface - enddo - mesh_ipVolume(i,e) = mesh_ipVolume(i,e) / FE_NipFaceNodes ! renormalize with interfaceNodeNum due to loop over them mesh_ipCenterOfGravity(:,i,e) = centerOfGravity + + if (calcIPvolume) then + do f = 1,FE_NipNeighbors(t) ! loop over interfaces of IP and add tetrahedra which connect to CoG + forall (n = 1:FE_NipFaceNodes) nPos(:,n) = mesh_subNodeCoord(:,FE_subNodeOnIPFace(n,f,i,t),e) + forall (n = 1:FE_NipFaceNodes, j = 1:Ntriangles) & ! start at each interface node and build valid triangles to cover interface + volume(j,n) = math_volTetrahedron(nPos(:,n), & ! calc volume of respective tetrahedron to CoG + nPos(:,1+mod(n-1 +j ,FE_NipFaceNodes)), & ! start at offset j + nPos(:,1+mod(n-1 +j+1,FE_NipFaceNodes)), & ! and take j's neighbor + centerOfGravity) + mesh_ipVolume(i,e) = mesh_ipVolume(i,e) + sum(volume) ! add contribution from this interface + enddo + mesh_ipVolume(i,e) = mesh_ipVolume(i,e) / FE_NipFaceNodes ! renormalize with interfaceNodeNum due to loop over them + endif + enddo enddo - return endsubroutine @@ -3268,13 +3285,13 @@ do dir = 1,3 ! check periodicity in direction minimumNodes = 0_pInt maximumNodes = 0_pInt - minCoord = minval(mesh_node(dir,:)) - maxCoord = maxval(mesh_node(dir,:)) + minCoord = minval(mesh_node0(dir,:)) + maxCoord = maxval(mesh_node0(dir,:)) do node = 1,mesh_Nnodes ! loop through all nodes and find surface nodes - if (abs(mesh_node(dir,node) - minCoord) <= tolerance) then + if (abs(mesh_node0(dir,node) - minCoord) <= tolerance) then minimumNodes(1) = minimumNodes(1) + 1_pInt minimumNodes(minimumNodes(1)+1) = node - elseif (abs(mesh_node(dir,node) - maxCoord) <= tolerance) then + elseif (abs(mesh_node0(dir,node) - maxCoord) <= tolerance) then maximumNodes(1) = maximumNodes(1) + 1_pInt maximumNodes(maximumNodes(1)+1) = node endif @@ -3290,7 +3307,7 @@ do dir = 1,3 ! check periodicity in direction cycle do n2 = 1,maximumNodes(1) maximumNode = maximumNodes(n2+1) - distance = abs(mesh_node(:,minimumNode) - mesh_node(:,maximumNode)) + distance = abs(mesh_node0(:,minimumNode) - mesh_node0(:,maximumNode)) if (sum(distance) - distance(dir) <= tolerance) then ! minimum possible distance (within tolerance) mesh_nodeTwins(dir,minimumNode) = maximumNode mesh_nodeTwins(dir,maximumNode) = minimumNode