diff --git a/src/CPFEM.f90 b/src/CPFEM.f90 index 3a1340764..06d1ca6a8 100644 --- a/src/CPFEM.f90 +++ b/src/CPFEM.f90 @@ -191,7 +191,8 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS else validCalculation if (debugCPFEM%extensive) print'(a,i8,1x,i2)', '<< CPFEM >> calculation for elFE ip ',elFE,ip call materialpoint_stressAndItsTangent(dt,[ip,ip],[elCP,elCP]) - call materialpoint_stressAndItsTangent2(dt,[ip,ip],[elCP,elCP]) + if (.not. terminallyIll) & + call materialpoint_stressAndItsTangent2(dt,[ip,ip],[elCP,elCP]) terminalIllness: if (terminallyIll) then diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index d164c33a2..607928c1b 100644 --- a/src/grid/spectral_utilities.f90 +++ b/src/grid/spectral_utilities.f90 @@ -816,11 +816,13 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,& homogenization_F = reshape(F,[3,3,product(grid(1:2))*grid3]) ! set materialpoint target F to estimated field call materialpoint_stressAndItsTangent(timeinc,[1,1],[1,product(grid(1:2))*grid3]) ! calculate P field - call materialpoint_stressAndItsTangent3(timeinc,[1,1],[1,product(grid(1:2))*grid3]) - call materialpoint_stressAndItsTangent2(timeinc,[1,1],[1,product(grid(1:2))*grid3]) + if (.not. terminallyIll) & + call materialpoint_stressAndItsTangent3(timeinc,[1,1],[1,product(grid(1:2))*grid3]) + if (.not. terminallyIll) & + call materialpoint_stressAndItsTangent2(timeinc,[1,1],[1,product(grid(1:2))*grid3]) P = reshape(homogenization_P, [3,3,grid(1),grid(2),grid3]) - P_av = sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt ! average of P + P_av = sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt call MPI_Allreduce(MPI_IN_PLACE,P_av,9,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,ierr) if (debugRotation) print'(/,a,/,2(3(2x,f12.4,1x)/),3(2x,f12.4,1x))', & ' Piola--Kirchhoff stress (lab) / MPa =', transpose(P_av)*1.e-6_pReal diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 9e80c6dda..f4f701d7c 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -309,30 +309,28 @@ subroutine materialpoint_stressAndItsTangent3(dt,FEsolving_execIP,FEsolving_exec logical, dimension(2) :: & doneAndHappy - if (.not. terminallyIll) then - !$OMP PARALLEL DO PRIVATE(ho,ph,ce) - do el = FEsolving_execElem(1),FEsolving_execElem(2) - if (terminallyIll) continue - do ip = FEsolving_execIP(1),FEsolving_execIP(2) - ce = (el-1)*discretization_nIPs + ip - ho = material_homogenizationID(ce) - call thermal_partition(ce) - do co = 1, homogenization_Nconstituents(ho) - ph = material_phaseID(co,ce) - if (.not. thermal_stress(dt,ph,material_phaseMemberAt(co,ip,el))) then - if (.not. terminallyIll) & ! so first signals terminally ill... - print*, ' Integration point ', ip,' at element ', el, ' terminally ill' - terminallyIll = .true. ! ...and kills all others - endif - enddo + !$OMP PARALLEL DO PRIVATE(ho,ph,ce) + do el = FEsolving_execElem(1),FEsolving_execElem(2) + if (terminallyIll) continue + do ip = FEsolving_execIP(1),FEsolving_execIP(2) + ce = (el-1)*discretization_nIPs + ip + ho = material_homogenizationID(ce) + call thermal_partition(ce) + do co = 1, homogenization_Nconstituents(ho) + ph = material_phaseID(co,ce) + if (.not. thermal_stress(dt,ph,material_phaseMemberAt(co,ip,el))) then + if (.not. terminallyIll) & ! so first signals terminally ill... + print*, ' Integration point ', ip,' at element ', el, ' terminally ill' + terminallyIll = .true. ! ...and kills all others + endif enddo enddo - !$OMP END PARALLEL DO - else - print'(/,a,/)', ' << HOMOG >> Material Point terminally ill' - endif + enddo + !$OMP END PARALLEL DO + end subroutine materialpoint_stressAndItsTangent3 + !-------------------------------------------------------------------------------------------------- !> @brief !-------------------------------------------------------------------------------------------------- @@ -351,22 +349,19 @@ subroutine materialpoint_stressAndItsTangent2(dt,FEsolving_execIP,FEsolving_exec doneAndHappy - if (.not. terminallyIll) then - !$OMP PARALLEL DO PRIVATE(ho,ce) - elementLooping3: do el = FEsolving_execElem(1),FEsolving_execElem(2) - IpLooping3: do ip = FEsolving_execIP(1),FEsolving_execIP(2) - ce = (el-1)*discretization_nIPs + ip - ho = material_homogenizationID(ce) - do co = 1, homogenization_Nconstituents(ho) - call crystallite_orientations(co,ip,el) - enddo - call mechanical_homogenize(dt,ce) - enddo IpLooping3 - enddo elementLooping3 - !$OMP END PARALLEL DO - else - print'(/,a,/)', ' << HOMOG >> Material Point terminally ill' - endif + !$OMP PARALLEL DO PRIVATE(ho,ce) + elementLooping3: do el = FEsolving_execElem(1),FEsolving_execElem(2) + IpLooping3: do ip = FEsolving_execIP(1),FEsolving_execIP(2) + ce = (el-1)*discretization_nIPs + ip + ho = material_homogenizationID(ce) + do co = 1, homogenization_Nconstituents(ho) + call crystallite_orientations(co,ip,el) + enddo + call mechanical_homogenize(dt,ce) + enddo IpLooping3 + enddo elementLooping3 + !$OMP END PARALLEL DO + end subroutine materialpoint_stressAndItsTangent2 diff --git a/src/mesh/FEM_utilities.f90 b/src/mesh/FEM_utilities.f90 index b388bc95b..402c6ddf7 100644 --- a/src/mesh/FEM_utilities.f90 +++ b/src/mesh/FEM_utilities.f90 @@ -163,8 +163,9 @@ subroutine utilities_constitutiveResponse(timeinc,P_av,forwardData) print'(/,a)', ' ... evaluating constitutive response ......................................' call materialpoint_stressAndItsTangent(timeinc,[1,mesh_maxNips],[1,mesh_NcpElems]) ! calculate P field - call materialpoint_stressAndItsTangent2(timeinc,[1,mesh_maxNips],[1,mesh_NcpElems]) - cutBack = .false. ! reset cutBack status + if (.not. terminallyIll) & + call materialpoint_stressAndItsTangent2(timeinc,[1,mesh_maxNips],[1,mesh_NcpElems]) + cutBack = .false. P_av = sum(homogenization_P,dim=3) * wgt call MPI_Allreduce(MPI_IN_PLACE,P_av,9,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,ierr)