From 5cedba0721417c08e51894f9d7a86e314b4e8d5f Mon Sep 17 00:00:00 2001 From: Fengbo Han Date: Thu, 27 Jul 2017 16:21:02 +0200 Subject: [PATCH 01/27] implemented yield stop criteria --- src/DAMASK_spectral.f90 | 85 +++++++++++++++++++++- src/crystallite.f90 | 4 +- src/math.f90 | 65 +++++++++++++++++ src/spectral_utilities.f90 | 142 ++++++++++++++++++++++++++++++++++++- 4 files changed, 290 insertions(+), 6 deletions(-) mode change 100644 => 100755 src/crystallite.f90 mode change 100644 => 100755 src/math.f90 mode change 100644 => 100755 src/spectral_utilities.f90 diff --git a/src/DAMASK_spectral.f90 b/src/DAMASK_spectral.f90 index dfa1746b2..bffc89790 100755 --- a/src/DAMASK_spectral.f90 +++ b/src/DAMASK_spectral.f90 @@ -75,7 +75,8 @@ program DAMASK_spectral FIELD_UNDEFINED_ID, & FIELD_MECH_ID, & FIELD_THERMAL_ID, & - FIELD_DAMAGE_ID + FIELD_DAMAGE_ID, & + utilities_calcPlasticity use spectral_mech_Basic use spectral_mech_AL use spectral_mech_Polarisation @@ -153,7 +154,21 @@ program DAMASK_spectral MPI_finalize, & MPI_allreduce, & PETScFinalize - +!-------------------------------------------------------------------------------------------------- +! variables related to stop criterion for yielding + real(pReal) :: plasticWorkOld, plasticWorkNew, & ! plastic work + eqTotalStrainOld, eqTotalStrainNew, & ! total equivalent strain + eqPlasticStrainOld, eqPlasticStrainNew, & ! total equivalent plastic strain + eqStressOld, eqStressNew , & ! equivalent stress + yieldStopValue + real(pReal), dimension(3,3) :: yieldStress,yieldStressOld,yieldStressNew, plasticStrainOld, plasticStrainNew, plasticStrainRate + integer(pInt) :: yieldResUnit = 0_pInt + character(len=13) :: stopFlag + logical :: yieldStop, yieldStopSatisfied + ! logical :: & + ! stop_totalStrain, & !< stop criterion + ! stop_plasticStrain, & + ! stop_plasticWork !-------------------------------------------------------------------------------------------------- ! init DAMASK (all modules) call CPFEM_initAll(el = 1_pInt, ip = 1_pInt) @@ -210,6 +225,8 @@ program DAMASK_spectral !-------------------------------------------------------------------------------------------------- ! reading the load case and assign values to the allocated data structure + yieldStop = .False. + yieldStopSatisfied = .False. rewind(FILEUNIT) do line = IO_read(FILEUNIT) @@ -284,6 +301,18 @@ program DAMASK_spectral temp_valueVector(j) = IO_floatValue(line,chunkPos,i+j) enddo loadCases(currentLoadCase)%rotation = math_plain9to33(temp_valueVector) + case('totalStrain', 'totalstrain') + yieldStop = .True. + stopFlag = 'totalStrain' + yieldStopValue = IO_floatValue(line,chunkPos,i+1_pInt) + case('plasticStrain', 'plasticstrain') + yieldStop = .True. + stopFlag = 'plasticStrain' + yieldStopValue = IO_floatValue(line,chunkPos,i+1_pInt) + case('plasticWork', 'plasticwork') + yieldStop = .True. + stopFlag = 'plasticWork' + yieldStopValue = IO_floatValue(line,chunkPos,i+1_pInt) end select enddo; enddo close(FILEUNIT) @@ -658,10 +687,60 @@ program DAMASK_spectral time = time + timeinc guess = .true. endif forwarding + + yieldCheck: if(yieldStop) then ! check if it yields or satisfies the certain stop condition + yieldStressOld = yieldStressNew + plasticStrainOld = plasticStrainNew + eqStressOld = eqStressNew + eqTotalStrainOld = eqTotalStrainNew + eqPlasticStrainOld = eqPlasticStrainNew + plasticWorkOld = plasticWorkNew + + call utilities_calcPlasticity(yieldStressNew, plasticStrainNew, eqStressNew, eqTotalStrainNew, eqPlasticStrainNew, plasticWorkNew) + + if(stopFlag .eq. 'totalStrain') then + if(eqTotalStrainNew > yieldStopValue) then + yieldStress = yieldStressOld * (eqTotalStrainNew - yieldStopValue)/(eqTotalStrainNew - eqTotalStrainOld) & ! linear interpolation of stress values + + yieldStressNew * (yieldStopValue - eqTotalStrainOld)/(eqTotalStrainNew - eqTotalStrainOld) + plasticStrainRate = (plasticStrainNew - plasticStrainOld)/(time - time0) ! calculate plastic strain rate + yieldStopSatisfied = .True. + endif + elseif(stopFlag .eq. 'plasticStrain') then + if(eqPlasticStrainNew > yieldStopValue) then + yieldStress = yieldStressOld * (eqPlasticStrainNew - yieldStopValue)/(eqPlasticStrainNew - eqPlasticStrainOld) & + + yieldStressNew * (yieldStopValue - eqPlasticStrainOld)/(eqPlasticStrainNew - eqPlasticStrainOld) + plasticStrainRate = (plasticStrainNew - plasticStrainOld)/(time - time0) + yieldStopSatisfied = .True. + endif + elseif(stopFlag .eq. 'plasticWork') then + if(plasticWorkNew > yieldStopValue) then + yieldStress = yieldStressOld * (plasticWorkNew - yieldStopValue)/(plasticWorkNew - plasticWorkOld) & + + yieldStressNew * (yieldStopValue - plasticWorkOld)/(plasticWorkNew - plasticWorkOld) + plasticStrainRate = (plasticStrainNew - plasticStrainOld)/(time - time0) + yieldStopSatisfied = .True. + endif + endif + endif yieldCheck + + if (yieldStopSatisfied) then ! when yield, write the yield stress and strain rate to file and quit the job + if (worldrank == 0) then + open(newunit=yieldResUnit,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//& + '.yield',form='FORMATTED',status='REPLACE') + do i = 1_pInt,3_pInt + write(yieldResUnit,*) (yieldStress(i,j), j=1,3) + enddo + do i = 1_pInt,3_pInt + write(yieldResUnit,*) (plasticStrainRate(i,j), j=1,3) + enddo + close(yieldResUnit) + call quit(0_pInt) + endif + endif enddo incLooping enddo loadCaseLooping - + + !-------------------------------------------------------------------------------------------------- ! report summary of whole calculation write(6,'(/,a)') ' ###########################################################################' diff --git a/src/crystallite.f90 b/src/crystallite.f90 old mode 100644 new mode 100755 index c5bd4d979..d9f01e3d1 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -55,14 +55,14 @@ module crystallite crystallite_Li0, & !< intermediate velocitiy grad at start of FE inc crystallite_partionedLi0,& !< intermediate velocity grad at start of homog inc crystallite_Fe, & !< current "elastic" def grad (end of converged time step) - crystallite_P !< 1st Piola-Kirchhoff stress per grain + crystallite_P, & !< 1st Piola-Kirchhoff stress per grain + crystallite_subF !< def grad to be reached at end of crystallite inc real(pReal), dimension(:,:,:,:,:), allocatable, private :: & crystallite_subFe0,& !< "elastic" def grad at start of crystallite inc crystallite_invFp, & !< inverse of current plastic def grad (end of converged time step) crystallite_subFp0,& !< plastic def grad at start of crystallite inc crystallite_invFi, & !< inverse of current intermediate def grad (end of converged time step) crystallite_subFi0,& !< intermediate def grad at start of crystallite inc - crystallite_subF, & !< def grad to be reached at end of crystallite inc crystallite_subF0, & !< def grad at start of crystallite inc crystallite_subLp0,& !< plastic velocity grad at start of crystallite inc crystallite_subLi0,& !< intermediate velocity grad at start of crystallite inc diff --git a/src/math.f90 b/src/math.f90 old mode 100644 new mode 100755 index 3a0533ddf..937f00d7a --- a/src/math.f90 +++ b/src/math.f90 @@ -145,6 +145,7 @@ module math math_sampleGaussVar, & math_symmetricEulers, & math_eigenvectorBasisSym33, & + math_eigenvectorBasisSym33_log, & math_eigenvectorBasisSym, & math_eigenValuesVectorsSym33, & math_eigenValuesVectorsSym, & @@ -2064,6 +2065,70 @@ function math_eigenvectorBasisSym33(m) end function math_eigenvectorBasisSym33 +!-------------------------------------------------------------------------------------------------- +!> @brief logarithm eigenvector basis of symmetric 33 matrix m +!-------------------------------------------------------------------------------------------------- +function math_eigenvectorBasisSym33_log(m) + + implicit none + real(pReal), dimension(3,3) :: math_eigenvectorBasisSym33_log + real(pReal), dimension(3) :: invariants, values + real(pReal), dimension(3,3), intent(in) :: m + real(pReal) :: P, Q, rho, phi + real(pReal), parameter :: TOL=1.e-14_pReal + real(pReal), dimension(3,3,3) :: N, EB + + invariants = math_invariantsSym33(m) + EB = 0.0_pReal + + P = invariants(2)-invariants(1)**2.0_pReal/3.0_pReal + Q = -2.0_pReal/27.0_pReal*invariants(1)**3.0_pReal+product(invariants(1:2))/3.0_pReal-invariants(3) + + threeSimilarEigenvalues: if(all(abs([P,Q]) < TOL)) then + values = invariants(1)/3.0_pReal +! this is not really correct, but at least the basis is correct + EB(1,1,1)=1.0_pReal + EB(2,2,2)=1.0_pReal + EB(3,3,3)=1.0_pReal + else threeSimilarEigenvalues + rho=sqrt(-3.0_pReal*P**3.0_pReal)/9.0_pReal + phi=acos(math_limit(-Q/rho*0.5_pReal,-1.0_pReal,1.0_pReal)) + values = 2.0_pReal*rho**(1.0_pReal/3.0_pReal)* & + [cos(phi/3.0_pReal), & + cos((phi+2.0_pReal*PI)/3.0_pReal), & + cos((phi+4.0_pReal*PI)/3.0_pReal) & + ] + invariants(1)/3.0_pReal + N(1:3,1:3,1) = m-values(1)*math_I3 + N(1:3,1:3,2) = m-values(2)*math_I3 + N(1:3,1:3,3) = m-values(3)*math_I3 + twoSimilarEigenvalues: if(abs(values(1)-values(2)) < TOL) then + EB(1:3,1:3,3)=math_mul33x33(N(1:3,1:3,1),N(1:3,1:3,2))/ & + ((values(3)-values(1))*(values(3)-values(2))) + EB(1:3,1:3,1)=math_I3-EB(1:3,1:3,3) + elseif(abs(values(2)-values(3)) < TOL) then twoSimilarEigenvalues + EB(1:3,1:3,1)=math_mul33x33(N(1:3,1:3,2),N(1:3,1:3,3))/ & + ((values(1)-values(2))*(values(1)-values(3))) + EB(1:3,1:3,2)=math_I3-EB(1:3,1:3,1) + elseif(abs(values(3)-values(1)) < TOL) then twoSimilarEigenvalues + EB(1:3,1:3,2)=math_mul33x33(N(1:3,1:3,1),N(1:3,1:3,3))/ & + ((values(2)-values(1))*(values(2)-values(3))) + EB(1:3,1:3,1)=math_I3-EB(1:3,1:3,2) + else twoSimilarEigenvalues + EB(1:3,1:3,1)=math_mul33x33(N(1:3,1:3,2),N(1:3,1:3,3))/ & + ((values(1)-values(2))*(values(1)-values(3))) + EB(1:3,1:3,2)=math_mul33x33(N(1:3,1:3,1),N(1:3,1:3,3))/ & + ((values(2)-values(1))*(values(2)-values(3))) + EB(1:3,1:3,3)=math_mul33x33(N(1:3,1:3,1),N(1:3,1:3,2))/ & + ((values(3)-values(1))*(values(3)-values(2))) + endif twoSimilarEigenvalues + endif threeSimilarEigenvalues + + math_eigenvectorBasisSym33_log = log(sqrt(values(1))) * EB(1:3,1:3,1) & + + log(sqrt(values(2))) * EB(1:3,1:3,2) & + + log(sqrt(values(3))) * EB(1:3,1:3,3) + +end function math_eigenvectorBasisSym33_log + !-------------------------------------------------------------------------------------------------- !> @brief rotational part from polar decomposition of 33 tensor m !-------------------------------------------------------------------------------------------------- diff --git a/src/spectral_utilities.f90 b/src/spectral_utilities.f90 old mode 100644 new mode 100755 index 0773a9065..d0d720155 --- a/src/spectral_utilities.f90 +++ b/src/spectral_utilities.f90 @@ -145,7 +145,8 @@ module spectral_utilities FIELD_UNDEFINED_ID, & FIELD_MECH_ID, & FIELD_THERMAL_ID, & - FIELD_DAMAGE_ID + FIELD_DAMAGE_ID, & + utilities_calcPlasticity private :: & utilities_getFreqDerivative @@ -1040,6 +1041,145 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,timeinc, & end subroutine utilities_constitutiveResponse +!-------------------------------------------------------------------------------------------------- +!> @brief calculates yield stress, plastic strain, total strain and their equivalent values +!-------------------------------------------------------------------------------------------------- +subroutine utilities_calcPlasticity(yieldStress, plasticStrain, eqStress, eqTotalStrain, eqPlasticStrain, plasticWork) + use crystallite, only: & + crystallite_Fp, & + crystallite_P, & + crystallite_subF + use material, only: & + homogenization_maxNgrains + use mesh, only: & + mesh_maxNips,& + mesh_NcpElems + use math, only: & + math_det33, & + math_inv33, & + math_mul33x33, & + math_trace33, & + math_transpose33, & + math_identity2nd, & + math_crossproduct, & + math_eigenvectorBasisSym, & + math_eigenvectorBasisSym33, & + math_eigenvectorBasisSym33_log, & + math_eigenValuesVectorsSym33 + + implicit none + + real(pReal), intent(out) :: eqStress, eqTotalStrain, eqPlasticStrain, plasticWork + real(pReal), dimension(3,3),intent(out) :: yieldStress, plasticStrain + real(pReal), dimension(3,3) :: cauchy, cauchy_av, P, Vp_av, V_total_av !< average + real(pReal), dimension(3) :: Values, S + real(pReal), dimension(3,3) :: Vectors, diag + real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & + Bp, Vp, F, F_temp, Fp, U, VT, R, V, V_total + real(pReal), dimension(15) :: WORK !< previous deformation gradient + integer(pInt) :: INFO, i, j, k, l, ierr + real(pReal) :: stress, stress_av, strain_total, strain_total_av, strain_plastic, strain_plastic_av, wgtm + + external :: dgesvd + + wgtm = 1.0/real(mesh_NcpElems*mesh_maxNips*homogenization_maxNgrains,pReal) + Vp_av = 0.0_pReal + V_total_av = 0.0_pReal + strain_plastic_av = 0.0_pReal + strain_total_av = 0.0_pReal + cauchy_av = 0.0_pReal + diag = 0.0_pReal + do k = 1_pInt, mesh_NcpElems; do j = 1_pInt, mesh_maxNips; do i = 1_pInt,homogenization_maxNgrains + F(1:3,1:3,i,j,k) = crystallite_subF(1:3,1:3,i,j,k) + Fp(1:3,1:3,i,j,k) = crystallite_Fp(1:3,1:3,i,j,k) + + P = crystallite_P(1:3,1:3,i,j,k) + cauchy = 1.0/math_det33(F(1:3,1:3,i,j,k))*math_mul33x33(P,transpose(F(1:3,1:3,i,j,k))) + cauchy_av = cauchy_av + cauchy + stress = Mises(cauchy, 'stress') + stress_av = stress_av + stress + + Bp(1:3,1:3,i,j,k) = math_mul33x33(Fp(1:3,1:3,i,j,k),math_transpose33(Fp(1:3,1:3,i,j,k))) ! plastic part of left Cauchy–Green deformation tensor + Vp(1:3,1:3,i,j,k) = math_eigenvectorBasisSym33_log(Bp(1:3,1:3,i,j,k)) + strain_plastic = Mises(Vp(1:3,1:3,i,j,k), 'strain') + strain_plastic_av = strain_plastic_av + strain_plastic + Vp_av = Vp_av + Vp(1:3,1:3,i,j,k) + + F_temp(1:3,1:3,i,j,k) = F(1:3,1:3,i,j,k) + call dgesvd ('A', 'A', 3, 3, F_temp(1:3,1:3,i,j,k), 3, & + S, U(1:3,1:3,i,j,k), 3, VT(1:3,1:3,i,j,k), 3, WORK, 15, INFO) ! singular value decomposition + R(1:3,1:3,i,j,k) = math_mul33x33(U(1:3,1:3,i,j,k), VT(1:3,1:3,i,j,k)) ! rotation of polar decomposition + V(1:3,1:3,i,j,k) = math_mul33x33(F(1:3,1:3,i,j,k),math_inv33(R(1:3,1:3,i,j,k))) + + call math_eigenValuesVectorsSym33(V(1:3,1:3,i,j,k),Values,Vectors) + do l = 1_pInt, 3_pInt + if (Values(l) < 0.0_pReal) then + Values(l) = -Values(l) + Vectors(1:3, l) = -Vectors(1:3, l) + endif + Values(l) = log(Values(l)) + diag(l,l) = Values(l) + enddo + if (dot_product(Vectors(1:3,1),Vectors(1:3,2)) /= 0) then + Vectors(1:3,2) = math_crossproduct(Vectors(1:3,3), Vectors(1:3,1)) + Vectors(1:3,2) = Vectors(1:3,2)/sqrt(dot_product(Vectors(1:3,2),Vectors(1:3,2))) + endif + if (dot_product(Vectors(1:3,2),Vectors(1:3,3)) /= 0) then + Vectors(1:3,3) = math_crossproduct(Vectors(1:3,1), Vectors(1:3,2)) + Vectors(1:3,3) = Vectors(1:3,3)/sqrt(dot_product(Vectors(1:3,3),Vectors(1:3,3))) + endif + if (dot_product(Vectors(1:3,3),Vectors(1:3,1)) /= 0) then + Vectors(1:3,1) = math_crossproduct(Vectors(1:3,2), Vectors(1:3,3)) + Vectors(1:3,1) = Vectors(1:3,1)/sqrt(dot_product(Vectors(1:3,1),Vectors(1:3,1))) + endif + + V_total(1:3,1:3,i,j,k) = REAL(math_mul33x33(Vectors, math_mul33x33(diag, transpose(Vectors)))) + strain_total = Mises(V_total(1:3,1:3,i,j,k), 'strain') + strain_total_av = strain_total_av + strain_total + + enddo; enddo; enddo + + + yieldStress = cauchy_av * wgtm + call MPI_Allreduce(MPI_IN_PLACE,yieldStress,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) + plasticStrain = Vp_av * wgtm + call MPI_Allreduce(MPI_IN_PLACE,plasticStrain,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) + + plasticWork = plasticWork + 0.5*(eqStress + stress_av * wgtm) * (strain_total_av * wgtm - eqTotalStrain) + call MPI_Allreduce(MPI_IN_PLACE,plasticWork,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) + eqStress = stress_av * wgtm + call MPI_Allreduce(MPI_IN_PLACE,eqStress,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) + eqTotalStrain = strain_total_av * wgtm + call MPI_Allreduce(MPI_IN_PLACE,eqTotalStrain,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) + eqPlasticStrain = strain_plastic_av * wgtm + call MPI_Allreduce(MPI_IN_PLACE,eqPlasticStrain,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) + +end subroutine utilities_calcPlasticity + +!-------------------------------------------------------------------------------------------------- +!> @brief vonMises equivalent values for symmetric part of requested strains and/or stresses. +!-------------------------------------------------------------------------------------------------- +function Mises(m, mType) + use math, only: & + math_trace33, & + math_identity2nd + + implicit none + real(pReal), intent(in), dimension(3,3) :: m + character(len=*), intent(in) :: mType + real(pReal) :: trace, Mises + real(pReal), dimension(3,3) :: eye, dev, symdev + + eye = math_identity2nd(3) + trace = math_trace33(m) + dev = m - trace/3.0*eye + symdev = 0.5*(dev+transpose(dev)) + if (mType .eq. 'stress') then + Mises = sqrt(3.0/2.0*sum(symdev*transpose(symdev))) + else + Mises = sqrt(2.0/3.0*sum(symdev*transpose(symdev))) + endif +end function Mises !-------------------------------------------------------------------------------------------------- !> @brief calculates forward rate, either guessing or just add delta/timeinc From 3f0284496687471360ade4866da052d84805972e Mon Sep 17 00:00:00 2001 From: Fengbo Han Date: Thu, 27 Jul 2017 16:24:56 +0200 Subject: [PATCH 02/27] implemented yield stop criteria --- src/DAMASK_spectral.f90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/DAMASK_spectral.f90 b/src/DAMASK_spectral.f90 index bffc89790..29d9b672a 100755 --- a/src/DAMASK_spectral.f90 +++ b/src/DAMASK_spectral.f90 @@ -696,7 +696,8 @@ program DAMASK_spectral eqPlasticStrainOld = eqPlasticStrainNew plasticWorkOld = plasticWorkNew - call utilities_calcPlasticity(yieldStressNew, plasticStrainNew, eqStressNew, eqTotalStrainNew, eqPlasticStrainNew, plasticWorkNew) + call utilities_calcPlasticity(yieldStressNew, plasticStrainNew, eqStressNew, eqTotalStrainNew, & + eqPlasticStrainNew, plasticWorkNew) if(stopFlag .eq. 'totalStrain') then if(eqTotalStrainNew > yieldStopValue) then From b33d7e05857bfd7fd8d4a195761436f06e6c63df Mon Sep 17 00:00:00 2001 From: Fengbo Han Date: Thu, 27 Jul 2017 16:28:33 +0200 Subject: [PATCH 03/27] implemented yield stop criteria --- src/DAMASK_spectral.f90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/DAMASK_spectral.f90 b/src/DAMASK_spectral.f90 index 29d9b672a..877970e74 100755 --- a/src/DAMASK_spectral.f90 +++ b/src/DAMASK_spectral.f90 @@ -699,21 +699,21 @@ program DAMASK_spectral call utilities_calcPlasticity(yieldStressNew, plasticStrainNew, eqStressNew, eqTotalStrainNew, & eqPlasticStrainNew, plasticWorkNew) - if(stopFlag .eq. 'totalStrain') then + if(stopFlag == 'totalStrain') then if(eqTotalStrainNew > yieldStopValue) then yieldStress = yieldStressOld * (eqTotalStrainNew - yieldStopValue)/(eqTotalStrainNew - eqTotalStrainOld) & ! linear interpolation of stress values + yieldStressNew * (yieldStopValue - eqTotalStrainOld)/(eqTotalStrainNew - eqTotalStrainOld) plasticStrainRate = (plasticStrainNew - plasticStrainOld)/(time - time0) ! calculate plastic strain rate yieldStopSatisfied = .True. endif - elseif(stopFlag .eq. 'plasticStrain') then + elseif(stopFlag == 'plasticStrain') then if(eqPlasticStrainNew > yieldStopValue) then yieldStress = yieldStressOld * (eqPlasticStrainNew - yieldStopValue)/(eqPlasticStrainNew - eqPlasticStrainOld) & + yieldStressNew * (yieldStopValue - eqPlasticStrainOld)/(eqPlasticStrainNew - eqPlasticStrainOld) plasticStrainRate = (plasticStrainNew - plasticStrainOld)/(time - time0) yieldStopSatisfied = .True. endif - elseif(stopFlag .eq. 'plasticWork') then + elseif(stopFlag == 'plasticWork') then if(plasticWorkNew > yieldStopValue) then yieldStress = yieldStressOld * (plasticWorkNew - yieldStopValue)/(plasticWorkNew - plasticWorkOld) & + yieldStressNew * (yieldStopValue - plasticWorkOld)/(plasticWorkNew - plasticWorkOld) From 414faa53d3a620eb217a8b340c97bc05bbc90c84 Mon Sep 17 00:00:00 2001 From: Fengbo Han Date: Thu, 27 Jul 2017 16:31:16 +0200 Subject: [PATCH 04/27] implemented yield stop criteria --- src/DAMASK_spectral.f90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/DAMASK_spectral.f90 b/src/DAMASK_spectral.f90 index 877970e74..fe7bdf2b3 100755 --- a/src/DAMASK_spectral.f90 +++ b/src/DAMASK_spectral.f90 @@ -721,12 +721,12 @@ program DAMASK_spectral yieldStopSatisfied = .True. endif endif - endif yieldCheck + endif yieldCheck - if (yieldStopSatisfied) then ! when yield, write the yield stress and strain rate to file and quit the job + if (yieldStopSatisfied) then ! when yield, write the yield stress and strain rate to file and quit the job if (worldrank == 0) then open(newunit=yieldResUnit,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//& - '.yield',form='FORMATTED',status='REPLACE') + '.yield',form='FORMATTED',status='REPLACE') do i = 1_pInt,3_pInt write(yieldResUnit,*) (yieldStress(i,j), j=1,3) enddo From d51fa10ae55453338b58fabec62932dddd80ccc6 Mon Sep 17 00:00:00 2001 From: Fengbo Han Date: Thu, 27 Jul 2017 16:33:05 +0200 Subject: [PATCH 05/27] implemented yield stop criteria --- src/spectral_utilities.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/spectral_utilities.f90 b/src/spectral_utilities.f90 index d0d720155..a1795c485 100755 --- a/src/spectral_utilities.f90 +++ b/src/spectral_utilities.f90 @@ -1174,7 +1174,7 @@ function Mises(m, mType) trace = math_trace33(m) dev = m - trace/3.0*eye symdev = 0.5*(dev+transpose(dev)) - if (mType .eq. 'stress') then + if (mType == 'stress') then Mises = sqrt(3.0/2.0*sum(symdev*transpose(symdev))) else Mises = sqrt(2.0/3.0*sum(symdev*transpose(symdev))) From afda166fd897aeea08b54a822792493c90e7a889 Mon Sep 17 00:00:00 2001 From: Fengbo Han Date: Tue, 1 Aug 2017 18:02:53 +0200 Subject: [PATCH 06/27] calculate platic strain by subtracting elastic strain from total strain --- src/spectral_utilities.f90 | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/src/spectral_utilities.f90 b/src/spectral_utilities.f90 index a1795c485..ecd47b594 100755 --- a/src/spectral_utilities.f90 +++ b/src/spectral_utilities.f90 @@ -1047,6 +1047,7 @@ end subroutine utilities_constitutiveResponse subroutine utilities_calcPlasticity(yieldStress, plasticStrain, eqStress, eqTotalStrain, eqPlasticStrain, plasticWork) use crystallite, only: & crystallite_Fp, & + crystallite_Fe, & crystallite_P, & crystallite_subF use material, only: & @@ -1075,7 +1076,7 @@ subroutine utilities_calcPlasticity(yieldStress, plasticStrain, eqStress, eqTota real(pReal), dimension(3) :: Values, S real(pReal), dimension(3,3) :: Vectors, diag real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & - Bp, Vp, F, F_temp, Fp, U, VT, R, V, V_total + Be, Ve, Vp, F, F_temp, Fe, Fp, U, VT, R, V, V_total real(pReal), dimension(15) :: WORK !< previous deformation gradient integer(pInt) :: INFO, i, j, k, l, ierr real(pReal) :: stress, stress_av, strain_total, strain_total_av, strain_plastic, strain_plastic_av, wgtm @@ -1092,6 +1093,7 @@ subroutine utilities_calcPlasticity(yieldStress, plasticStrain, eqStress, eqTota do k = 1_pInt, mesh_NcpElems; do j = 1_pInt, mesh_maxNips; do i = 1_pInt,homogenization_maxNgrains F(1:3,1:3,i,j,k) = crystallite_subF(1:3,1:3,i,j,k) Fp(1:3,1:3,i,j,k) = crystallite_Fp(1:3,1:3,i,j,k) + Fe(1:3,1:3,i,j,k) = crystallite_Fe(1:3,1:3,i,j,k) P = crystallite_P(1:3,1:3,i,j,k) cauchy = 1.0/math_det33(F(1:3,1:3,i,j,k))*math_mul33x33(P,transpose(F(1:3,1:3,i,j,k))) @@ -1099,12 +1101,6 @@ subroutine utilities_calcPlasticity(yieldStress, plasticStrain, eqStress, eqTota stress = Mises(cauchy, 'stress') stress_av = stress_av + stress - Bp(1:3,1:3,i,j,k) = math_mul33x33(Fp(1:3,1:3,i,j,k),math_transpose33(Fp(1:3,1:3,i,j,k))) ! plastic part of left Cauchy–Green deformation tensor - Vp(1:3,1:3,i,j,k) = math_eigenvectorBasisSym33_log(Bp(1:3,1:3,i,j,k)) - strain_plastic = Mises(Vp(1:3,1:3,i,j,k), 'strain') - strain_plastic_av = strain_plastic_av + strain_plastic - Vp_av = Vp_av + Vp(1:3,1:3,i,j,k) - F_temp(1:3,1:3,i,j,k) = F(1:3,1:3,i,j,k) call dgesvd ('A', 'A', 3, 3, F_temp(1:3,1:3,i,j,k), 3, & S, U(1:3,1:3,i,j,k), 3, VT(1:3,1:3,i,j,k), 3, WORK, 15, INFO) ! singular value decomposition @@ -1137,6 +1133,14 @@ subroutine utilities_calcPlasticity(yieldStress, plasticStrain, eqStress, eqTota strain_total = Mises(V_total(1:3,1:3,i,j,k), 'strain') strain_total_av = strain_total_av + strain_total + Be(1:3,1:3,i,j,k) = math_mul33x33(Fe(1:3,1:3,i,j,k),math_transpose33(Fe(1:3,1:3,i,j,k))) ! plastic part of left Cauchy–Green deformation tensor + Ve(1:3,1:3,i,j,k) = math_eigenvectorBasisSym33_log(Be(1:3,1:3,i,j,k)) + Vp(1:3,1:3,i,j,k) = V_total(1:3,1:3,i,j,k) - Ve(1:3,1:3,i,j,k) + strain_plastic = Mises(Vp(1:3,1:3,i,j,k), 'strain') + strain_plastic_av = strain_plastic_av + strain_plastic + Vp_av = Vp_av + Vp(1:3,1:3,i,j,k) + + flush(6) enddo; enddo; enddo From 36c370e6682199c1ba2f02af36799f54230921a2 Mon Sep 17 00:00:00 2001 From: Fengbo Han Date: Tue, 8 Aug 2017 17:25:38 +0200 Subject: [PATCH 07/27] implemented fast yield surface fitting with yield stop criteria --- lib/damask/util.py | 3 +- processing/misc/yieldSurfaceFast.py | 1390 +++++++++++++++++++++++++++ src/spectral_utilities.f90 | 1 - 3 files changed, 1392 insertions(+), 2 deletions(-) mode change 100644 => 100755 lib/damask/util.py create mode 100755 processing/misc/yieldSurfaceFast.py diff --git a/lib/damask/util.py b/lib/damask/util.py old mode 100644 new mode 100755 index 54de9acf9..dfe50c06e --- a/lib/damask/util.py +++ b/lib/damask/util.py @@ -39,7 +39,7 @@ def srepr(arg,glue = '\n'): if (not hasattr(arg, "strip") and hasattr(arg, "__getitem__") or hasattr(arg, "__iter__")): - return glue.join(srepr(x) for x in arg) + return glue.join(str(x) for x in arg) return arg if isinstance(arg,str) else repr(arg) # ----------------------------- @@ -221,6 +221,7 @@ def leastsqBound(func, x0, args=(), bounds=None, Dfun=None, full_output=0, def _check_func(checker, argname, thefunc, x0, args, numinputs, output_shape=None): + from numpy import shape """The same as that of minpack.py""" res = np.atleast_1d(thefunc(*((x0[:numinputs],) + args))) if (output_shape is not None) and (shape(res) != output_shape): diff --git a/processing/misc/yieldSurfaceFast.py b/processing/misc/yieldSurfaceFast.py new file mode 100755 index 000000000..7b00124e7 --- /dev/null +++ b/processing/misc/yieldSurfaceFast.py @@ -0,0 +1,1390 @@ +#!/usr/bin/env python2.7 +# -*- coding: UTF-8 no BOM -*- + +import threading,time,os +import numpy as np +from optparse import OptionParser +import damask +from damask.util import leastsqBound + +scriptName = os.path.splitext(os.path.basename(__file__))[0] +scriptID = ' '.join([scriptName,damask.version]) + +def runFit(exponent, eqStress, dimension, criterion): + global threads, myFit, myLoad + global fitResidual + global Guess, dDim + + dDim = dimension - 3 + nParas = len(fitCriteria[criterion]['bound'][dDim]) + nExpo = fitCriteria[criterion]['nExpo'] + + if exponent > 0.0: # User defined exponents + nParas = nParas-nExpo + fitCriteria[criterion]['bound'][dDim] = fitCriteria[criterion]['bound'][dDim][:nParas] + + for i in range(nParas): + temp = fitCriteria[criterion]['bound'][dDim][i] + if fitCriteria[criterion]['bound'][dDim][i] == (None,None): + Guess.append(1.0) + else: + g = (temp[0]+temp[1])/2.0 + if g == 0: g = temp[1]*0.5 + Guess.append(g) + + myLoad = Loadcase(options.load[0],options.load[1],options.load[2],options.flag,options.yieldValue, + nSet = 10, dimension = dimension, vegter = options.criterion=='vegter') + + + myFit = Criterion(exponent,eqStress, dimension, criterion) + for t in range(options.threads): + threads.append(myThread(t)) + threads[t].start() + + for t in range(options.threads): + threads[t].join() + damask.util.croak('Residuals') + damask.util.croak(fitResidual) + +def principalStresses(sigmas): + """ + Computes principal stresses (i.e. eigenvalues) for a set of Cauchy stresses. + + sorted in descending order. + """ + lambdas=np.zeros(0,'d') + for i in range(np.shape(sigmas)[1]): + eigenvalues = np.linalg.eigvalsh(sym6toT33(sigmas[:,i])) + lambdas = np.append(lambdas,np.sort(eigenvalues)[::-1]) #append eigenvalues in descending order + lambdas = np.transpose(lambdas.reshape(np.shape(sigmas)[1],3)) + return lambdas + +def principalStress(p): + I = invariant(p) + + I1s3I2= (I[0]**2 - 3.0*I[1])**0.5 + numer = 2.0*I[0]**3 - 9.0*I[0]*I[1] + 27.0*I[2] + denom = 2.0*I1s3I2**3 + cs = numer/denom + + phi = np.arccos(cs)/3.0 + t1 = I[0]/3.0; t2 = 2.0/3.0*I1s3I2 + return np.array( [t1 + t2*np.cos(phi), + t1 + t2*np.cos(phi+np.pi*2.0/3.0), + t1 + t2*np.cos(phi+np.pi*4.0/3.0)]) + +def principalStrs_Der(p, s, dim, Karafillis=False): + """Derivative of principal stress with respect to stress""" + third = 1.0/3.0 + third2 = 2.0*third + + I = invariant(p) + I1s3I2= np.sqrt(I[0]**2 - 3.0*I[1]) + numer = 2.0*I[0]**3 - 9.0*I[0]*I[1] + 27.0*I[2] + denom = 2.0*I1s3I2**3 + cs = numer/denom + phi = np.arccos(cs)/3.0 + + dphidcs = -third/np.sqrt(1.0 - cs**2) + dcsddenom = 0.5*numer*(-1.5)*I1s3I2**(-5.0) + dcsdI1 = (6.0*I[0]**2 - 9.0*I[1])*denom + dcsddenom*(2.0*I[0]) + dcsdI2 = ( - 9.0*I[0])*denom + dcsddenom*(-3.0) + dcsdI3 = 27.0*denom + dphidI1, dphidI2, dphidI3 = dphidcs*dcsdI1, dphidcs*dcsdI2, dphidcs*dcsdI3 + + dI1s3I2dI1 = I[0]/I1s3I2 + dI1s3I2dI2 = -1.5/I1s3I2 + tcoeff = third2*I1s3I2 + + dSidIj = lambda theta : ( tcoeff*(-np.sin(theta))*dphidI1 + third2*dI1s3I2dI1*np.cos(theta) + third, + tcoeff*(-np.sin(theta))*dphidI2 + third2*dI1s3I2dI2*np.cos(theta), + tcoeff*(-np.sin(theta))*dphidI3) + dSdI = np.array([dSidIj(phi),dSidIj(phi+np.pi*2.0/3.0),dSidIj(phi+np.pi*4.0/3.0)]) # i=1,2,3; j=1,2,3 + +# calculate the derivation of principal stress with regards to the anisotropic coefficients + one = np.ones_like(s); zero = np.zeros_like(s); num = len(s) + dIdp = np.array([[one, one, one, zero, zero, zero], + [p[1]+p[2], p[2]+p[0], p[0]+p[1], -2.0*p[3], -2.0*p[4], -2.0*p[5]], + [p[1]*p[2]-p[4]**2, p[2]*p[0]-p[5]**2, p[0]*p[1]-p[3]**2, + -2.0*p[3]*p[2]+2.0*p[4]*p[5], -2.0*p[4]*p[0]+2.0*p[5]*p[3], -2.0*p[5]*p[1]+2.0*p[3]*p[4]] ]) + if Karafillis: + dpdc = np.array([[zero,s[0]-s[2],s[0]-s[1]], [s[1]-s[2],zero,s[1]-s[0]], [s[2]-s[1],s[2]-s[0],zero]])/3.0 + dSdp = np.array([np.dot(dSdI[:,:,i],dIdp[:,:,i]).T for i in range(num)]).T + if dim == 2: + temp = np.vstack([dSdp[:,3]*s[3]]).T.reshape(num,1,3).T + else: + temp = np.vstack([dSdp[:,3]*s[3],dSdp[:,4]*s[4],dSdp[:,5]*s[5]]).T.reshape(num,3,3).T + + return np.concatenate((np.array([np.dot(dSdp[:,0:3,i], dpdc[:,:,i]).T for i in range(num)]).T, + temp), axis=1) + else: + if dim == 2: + dIdc=np.array([[-dIdp[i,0]*s[1], -dIdp[i,1]*s[0], -dIdp[i,1]*s[2], + -dIdp[i,2]*s[1], -dIdp[i,2]*s[0], -dIdp[i,0]*s[2], + dIdp[i,3]*s[3] ] for i in range(3)]) + else: + dIdc=np.array([[-dIdp[i,0]*s[1], -dIdp[i,1]*s[0], -dIdp[i,1]*s[2], + -dIdp[i,2]*s[1], -dIdp[i,2]*s[0], -dIdp[i,0]*s[2], + dIdp[i,3]*s[3], dIdp[i,4]*s[4], dIdp[i,5]*s[5] ] for i in range(3)]) + return np.array([np.dot(dSdI[:,:,i],dIdc[:,:,i]).T for i in range(num)]).T + +def invariant(sigmas): + I = np.zeros(3) + s11,s22,s33,s12,s23,s31 = sigmas + I[0] = s11 + s22 + s33 + I[1] = s11*s22 + s22*s33 + s33*s11 - s12**2 - s23**2 - s31**2 + I[2] = s11*s22*s33 + 2.0*s12*s23*s31 - s12**2*s33 - s23**2*s11 - s31**2*s22 + return I + +def math_ln(x): + return np.log(x + 1.0e-32) + +def sym6toT33(sym6): + """Shape the symmetric stress tensor(6) into (3,3)""" + return np.array([[sym6[0],sym6[3],sym6[5]], + [sym6[3],sym6[1],sym6[4]], + [sym6[5],sym6[4],sym6[2]]]) + +def t33toSym6(t33): + """Shape the stress tensor(3,3) into symmetric (6)""" + return np.array([ t33[0,0], + t33[1,1], + t33[2,2], + (t33[0,1] + t33[1,0])/2.0, # 0 3 5 + (t33[1,2] + t33[2,1])/2.0, # * 1 4 + (t33[2,0] + t33[0,2])/2.0,]) # * * 2 + +class Criteria(object): + def __init__(self, criterion, uniaxialStress,exponent, dimension): + self.stress0 = uniaxialStress + if exponent < 0.0: # Fitting exponent m + self.mFix = [False, exponent] + else: # fixed exponent m + self.mFix = [True, exponent] + self.func = fitCriteria[criterion]['func'] + self.criteria = criterion + self.dim = dimension + def fun(self, paras, ydata, sigmas): + return self.func(self.stress0, paras, sigmas,self.mFix,self.criteria,self.dim) + def jac(self, paras, ydata, sigmas): + return self.func(self.stress0, paras, sigmas,self.mFix,self.criteria,self.dim,Jac=True) + +class Vegter(object): + """Vegter yield criterion""" + + def __init__(self, refPts, refNormals,nspace=11): + self.refPts, self.refNormals = self._getRefPointsNormals(refPts, refNormals) + self.hingePts = self._getHingePoints() + self.nspace = nspace + def _getRefPointsNormals(self,refPtsQtr,refNormalsQtr): + if len(refPtsQtr) == 12: + refPts = refPtsQtr + refNormals = refNormalsQtr + else: + refPts = np.empty([13,2]) + refNormals = np.empty([13,2]) + refPts[12] = refPtsQtr[0] + refNormals[12] = refNormalsQtr[0] + for i in range(3): + refPts[i] = refPtsQtr[i] + refPts[i+3] = refPtsQtr[3-i][::-1] + refPts[i+6] =-refPtsQtr[i] + refPts[i+9] =-refPtsQtr[3-i][::-1] + refNormals[i] = refNormalsQtr[i] + refNormals[i+3] = refNormalsQtr[3-i][::-1] + refNormals[i+6] =-refNormalsQtr[i] + refNormals[i+9] =-refNormalsQtr[3-i][::-1] + return refPts,refNormals + + def _getHingePoints(self): + """ + Calculate the hinge point B according to the reference points A,C and the normals n,m + + refPoints = np.array([[p1_x, p1_y], [p2_x, p2_y]]); + refNormals = np.array([[n1_x, n1_y], [n2_x, n2_y]]) + """ + def hingPoint(points, normals): + A1 = points[0][0]; A2 = points[0][1] + C1 = points[1][0]; C2 = points[1][1] + n1 = normals[0][0]; n2 = normals[0][1] + m1 = normals[1][0]; m2 = normals[1][1] + B1 = (m2*(n1*A1 + n2*A2) - n2*(m1*C1 + m2*C2))/(n1*m2-m1*n2) + B2 = (n1*(m1*C1 + m2*C2) - m1*(n1*A1 + n2*A2))/(n1*m2-m1*n2) + return np.array([B1,B2]) + return np.array([hingPoint(self.refPts[i:i+2],self.refNormals[i:i+2]) for i in range(len(self.refPts)-1)]) + + def getBezier(self): + def bezier(R,H): + b = [] + for mu in np.linspace(0.0,1.0,self.nspace): + b.append(np.array(R[0]*np.ones_like(mu) + 2.0*mu*(H - R[0]) + mu**2*(R[0]+R[1] - 2.0*H))) + return b + return np.array([bezier(self.refPts[i:i+2],self.hingePts[i]) for i in range(len(self.refPts)-1)]) + +def VetgerCriterion(stress,lankford, rhoBi0, theta=0.0): + """0-pure shear; 1-uniaxial; 2-plane strain; 3-equi-biaxial""" + def getFourierParas(r): + # get the value after Fourier transformation + nset = len(r) + lmatrix = np.empty([nset,nset]) + theta = np.linspace(0.0,np.pi/2,nset) + for i,th in enumerate(theta): + lmatrix[i] = np.array([np.cos(2*j*th) for j in range(nset)]) + return np.linalg.solve(lmatrix, r) + + nps = len(stress) + if nps%4 != 0: + damask.util.croak('Warning: the number of stress points is uncorrect, stress points of %s are missing in set %i'%( + ['eq-biaxial, plane strain & uniaxial', 'eq-biaxial & plane strain','eq-biaxial'][nps%4-1],nps/4+1)) + else: + nset = nps/4 + strsSet = stress.reshape(nset,4,2) + refPts = np.empty([4,2]) + + fouriercoeffs = np.array([np.cos(2.0*i*theta) for i in range(nset)]) + for i in range(2): + refPts[3,i] = sum(strsSet[:,3,i])/nset + for j in range(3): + refPts[j,i] = np.dot(getFourierParas(strsSet[:,j,i]), fouriercoeffs) + + +def Tresca(eqStress=None, #not needed/supported + paras=None, + sigmas=None, + mFix=None, #not needed/supported + criteria=None, #not needed/supported + dim=3, + Jac=False): + """ + Tresca yield criterion + + the fitted parameter is paras(sigma0) + """ + if not Jac: + lambdas = principalStresses(sigmas) + r = np.amax(np.array([abs(lambdas[2,:]-lambdas[1,:]),\ + abs(lambdas[1,:]-lambdas[0,:]),\ + abs(lambdas[0,:]-lambdas[2,:])]),0) - paras + return r.ravel() + else: + return -np.ones(len(sigmas)) + +def Cazacu_Barlat(eqStress=None, + paras=None, + sigmas=None, + mFix=None,#not needed/supported + criteria=None, + dim=3, #2D also possible + Jac=False): + """ + Cazacu-Barlat (CB) yield criterion + + the fitted parameters are: + a1,a2,a3,a6; b1,b2,b3,b4,b5,b10; c for plane stress + a1,a2,a3,a4,a5,a6; b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11; c: for general case + mFix is ignored + """ + s11,s22,s33,s12,s23,s31 = sigmas + if dim == 2: + (a1,a2,a3,a4), (b1,b2,b3,b4,b5,b10), c = paras[0:4],paras[4:10],paras[10] + a5 = a6 = b6 = b7 = b8 = b9 = b11 = 0.0 + s33 = s23 = s31 = np.zeros_like(s11) + else: + (a1,a2,a3,a4,a5,a6), (b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11), c = paras[0:6],paras[6:17],paras[17] + + s1_2, s2_2, s3_2, s12_2, s23_2, s31_2 = np.array([s11,s22,s33,s12,s23,s31])**2 + s1_3, s2_3, s3_3, s123, s321 = s11*s1_2, s22*s2_2, s33*s3_2,s11*s22*s33, s12*s23*s31 + d12_2,d23_2,d31_2 = (s11-s22)**2, (s22-s33)**2, (s33-s11)**2 + + J20 = ( a1*d12_2 + a2*d23_2 + a3*d31_2 )/6.0 + a4*s12_2 + a5*s23_2 + a6*s31_2 + J30 = ( (b1 +b2 )*s1_3 + (b3 +b4 )*s2_3 + ( b1+b4-b2 + b1+b4-b3 )*s3_3 )/27.0- \ + ( (b1*s22+b2*s33)*s1_2 + (b3*s33+b4*s11)*s2_2 + ((b1+b4-b2)*s11 + (b1+b4-b3)*s22)*s3_2 )/9.0 + \ + ( (b1+b4)*s123/9.0 + b11*s321 )*2.0 - \ + ( ( 2.0*b9 *s22 - b8*s33 - (2.0*b9 -b8)*s11 )*s31_2 + + ( 2.0*b10*s33 - b5*s22 - (2.0*b10-b5)*s11 )*s12_2 + + ( (b6+b7)*s11 - b6*s22 - b7*s33 )*s23_2 + )/3.0 + f0 = J20**3 - c*J30**2 + r = f0**(1.0/6.0)*np.sqrt(3.0)/eqStress + + if not Jac: + return (r - 1.0).ravel() + else: + drdf = r/f0/6.0 + dj2, dj3 = drdf*3.0*J20**2, -drdf*2.0*J30*c + jc = -drdf*J30**2 + + ja1,ja2,ja3 = dj2*d12_2/6.0, dj2*d23_2/6.0, dj2*d31_2/6.0 + ja4,ja5,ja6 = dj2*s12_2, dj2*s23_2, dj2*s31_2 + jb1 = dj3*( (s1_3 + 2.0*s3_3)/27.0 - s22*s1_2/9.0 - (s11+s22)*s3_2/9.0 + s123/4.5 ) + jb2 = dj3*( (s1_3 - s3_3)/27.0 - s33*s1_2/9.0 + s11 *s3_2/9.0 ) + jb3 = dj3*( (s2_3 - s3_3)/27.0 - s33*s2_2/9.0 + s22 *s3_2/9.0 ) + jb4 = dj3*( (s2_3 + 2.0*s3_3)/27.0 - s11*s2_2/9.0 - (s11+s22)*s3_2/9.0 + s123/4.5 ) + + jb5, jb10 = dj3*(s22 - s11)*s12_2/3.0, dj3*(s11 - s33)*s12_2/1.5 + jb6, jb7 = dj3*(s22 - s11)*s23_2/3.0, dj3*(s33 - s11)*s23_2/3.0 + jb8, jb9 = dj3*(s33 - s11)*s31_2/3.0, dj3*(s11 - s22)*s31_2/1.5 + jb11 = dj3*s321*2.0 + if dim == 2: + return np.vstack((ja1,ja2,ja3,ja4,jb1,jb2,jb3,jb4,jb5,jb10,jc)).T + else: + return np.vstack((ja1,ja2,ja3,ja4,ja5,ja6,jb1,jb2,jb3,jb4,jb5,jb6,jb7,jb8,jb9,jb10,jb11,jc)).T + +def Drucker(eqStress=None,#not needed/supported + paras=None, + sigmas=None, + mFix=None, #not needed/supported + criteria=None, + dim=3, + Jac=False): + """ + Drucker yield criterion + + the fitted parameters are + sigma0, C_D for Drucker(p=1); + sigma0, C_D, p for general Drucker + eqStress, mFix are invalid inputs + """ + if criteria == 'drucker': + sigma0, C_D= paras + p = 1.0 + else: + sigma0, C_D = paras[0:2] + if mFix[0]: p = mFix[1] + else: p = paras[-1] + I = invariant(sigmas) + J = np.zeros([3]) + J[1] = I[0]**2/3.0 - I[1] + J[2] = I[0]**3/13.5 - I[0]*I[1]/3.0 + I[2] + J2_3p = J[1]**(3.0*p) + J3_2p = J[2]**(2.0*p) + left = J2_3p - C_D*J3_2p + r = left**(1.0/(6.0*p))*3.0**0.5/sigma0 + + if not Jac: + return (r - 1.0).ravel() + else: + drdl = r/left/(6.0*p) + if criteria == 'drucker': + return np.vstack((-r/sigma0, -drdl*J3_2p)).T + else: + dldp = 3.0*J2_3p*math_ln(J[1]) - 2.0*C_D*J3_2p*math_ln(J[2]) + jp = drdl*dldp + r*math_ln(left)/(-6.0*p*p) + + if mFix[0]: return np.vstack((-r/sigma0, -drdl*J3_2p)).T + else: return np.vstack((-r/sigma0, -drdl*J3_2p, jp)).T + +def Hill1948(eqStress=None,#not needed/supported + paras=None, + sigmas=None, + mFix=None, #not needed/supported + criteria=None,#not needed/supported + dim=3, + Jac=False): + """ + Hill 1948 yield criterion + + the fitted parameters are: + F, G, H, L, M, N for 3D + F, G, H, N for 2D + """ + s11,s22,s33,s12,s23,s31 = sigmas + if dim == 2: # plane stress + jac = np.array([ s22**2, s11**2, (s11-s22)**2, 2.0*s12**2]) + else: # general case + jac = np.array([(s22-s33)**2,(s33-s11)**2,(s11-s22)**2, 2.0*s23**2,2.0*s31**2,2.0*s12**2]) + + if not Jac: + return (np.dot(paras,jac)/2.0-0.5).ravel() + else: + return jac.T + +def Hill1979(eqStress=None,#not needed/supported + paras=None, + sigmas=None, + mFix=None, + criteria=None,#not needed/supported + dim=3, + Jac=False): + """ + Hill 1979 yield criterion + + the fitted parameters are: f,g,h,a,b,c,m + """ + if mFix[0]: + m = mFix[1] + else: + m = paras[-1] + + coeff = paras[0:6] + s = principalStresses(sigmas) + diffs = np.array([s[1]-s[2], s[2]-s[0], s[0]-s[1],\ + 2.0*s[0]-s[1]-s[2], 2.0*s[1]-s[2]-s[0], 2.0*s[2]-s[0]-s[1]])**2 + + diffsm = diffs**(m/2.0) + left = np.dot(coeff,diffsm) + r = (0.5*left)**(1.0/m)/eqStress #left = base**mi + + if not Jac: + return (r-1.0).ravel() + else: + drdl, dldm = r/left/m, np.dot(coeff,diffsm*math_ln(diffs))*0.5 + jm = drdl*dldm + r*math_ln(0.5*left)*(-1.0/m/m) #/(-m**2) + + if mFix[0]: return np.vstack((drdl*diffsm)).T + else: return np.vstack((drdl*diffsm, jm)).T + +def Hosford(eqStress=None, + paras=None, + sigmas=None, + mFix=None, + criteria=None, + dim=3, + Jac=False): + """ + Hosford family criteria + + the fitted parameters are: + von Mises: sigma0 + Hershey: (1) sigma0, a, when a is not fixed; (2) sigma0, when a is fixed + general Hosford: (1) F,G,H, a, when a is not fixed; (2) F,G,H, when a is fixed + """ + if criteria == 'vonmises': + sigma0 = paras + coeff = np.ones(3) + a = 2.0 + elif criteria == 'hershey': + sigma0 = paras[0] + coeff = np.ones(3) + if mFix[0]: a = mFix[1] + else: a = paras[1] + else: + sigma0 = eqStress + coeff = paras[0:3] + if mFix[0]: a = mFix[1] + else: a = paras[3] + + s = principalStresses(sigmas) + diffs = np.array([s[1]-s[2], s[2]-s[0], s[0]-s[1]])**2 + diffsm = diffs**(a/2.0) + left = np.dot(coeff,diffsm) + r = (0.5*left)**(1.0/a)/sigma0 + + if not Jac: + return (r-1.0).ravel() + else: + if criteria == 'vonmises': # von Mises + return -r/sigma0 + else: + drdl, dlda = r/left/a, np.dot(coeff,diffsm*math_ln(diffs))*0.5 + ja = drdl*dlda + r*math_ln(0.5*left)*(-1.0/a/a) + if criteria == 'hershey': # Hershey + if mFix[0]: return -r/sigma0 + else: return np.vstack((-r/sigma0, ja)).T + else: # Anisotropic Hosford + if mFix[0]: return np.vstack((drdl*diffsm)).T + else: return np.vstack((drdl*diffsm, ja)).T + +def Barlat1989(eqStress=None, + paras=None, + sigmas=None, + mFix=None, + criteria=None, + dim=3, + Jac=False): + """ + Barlat-Lian 1989 yield criteria + + the fitted parameters are: + Anisotropic: a, h, p, m; m is optional + """ + a, h, p = paras[0:3] + if mFix[0]: m = mFix[1] + else: m = paras[-1] + + c = 2.0-a + s11,s22,s12 = sigmas[0], sigmas[1], sigmas[3] + k1,k2 = 0.5*(s11 + h*s22), (0.25*(s11 - h*s22)**2 + (p*s12)**2)**0.5 + fs = np.array([ (k1+k2)**2, (k1-k2)**2, 4.0*k2**2 ]); fm = fs**(m/2.0) + left = np.dot(np.array([a,a,c]),fm) + r = (0.5*left)**(1.0/m)/eqStress + + if not Jac: + return (r-1.0).ravel() + else: + dk1dh = 0.5*s22 + dk2dh, dk2dp = 0.25*(s11-h*s22)*(-s22)/k2, p*s12**2/k2 + dlda, dldc = fm[0]+fm[1], fm[2] + fm1 = fs**(m/2.0-1.0)*m + dldk1, dldk2 = a*fm1[0]*(k1+k2)+a*fm1[1]*(k1-k2), a*fm1[0]*(k1+k2)-a*fm1[1]*(k1-k2)+c*fm1[2]*k2*4.0 + drdl, drdm = r/m/left, r*math_ln(0.5*left)*(-1.0/m/m) + dldm = np.dot(np.array([a,a,c]),fm*math_ln(fs))*0.5 + + ja,jc = drdl*dlda, drdl*dldc + jh,jp = drdl*(dldk1*dk1dh + dldk2*dk2dh), drdl*dldk2*dk2dp + jm = drdl*dldm + drdm + + if mFix[0]: return np.vstack((ja,jc,jh,jp)).T + else: return np.vstack((ja,jc,jh,jp,jm)).T + +def Barlat1991(eqStress, paras, sigmas, mFix, criteria, dim, Jac=False): + """ + Barlat 1991 criteria + + the fitted parameters are: + Anisotropic: a, b, c, f, g, h, m for 3D + a, b, c, h, m for plane stress + m is optional + """ + if dim == 2: coeff = paras[0:4] # plane stress + else: coeff = paras[0:6] # general case + if mFix[0]: m = mFix[1] + else: m = paras[-1] + + s11,s22,s33,s12,s23,s31 = sigmas + if dim == 2: + dXdx = np.array([s22,-s11,s11-s22,s12]) + A,B,C,H = np.array(coeff)[:,None]*dXdx; F=G=0.0 + else: + dXdx = np.array([s22-s33,s33-s11,s11-s22,s23,s31,s12]) + A,B,C,F,G,H = np.array(coeff)[:,None]*dXdx + + I2 = (F*F + G*G + H*H)/3.0+ ((A-C)**2+(C-B)**2+(B-A)**2)/54.0 + I3 = (C-B)*(A-C)*(B-A)/54.0 + F*G*H - ((C-B)*F*F + (A-C)*G*G + (B-A)*H*H)/6.0 + phi1 = np.arccos(I3/I2**1.5)/3.0 + np.pi/6.0; absc1 = 2.0*np.abs(np.cos(phi1)) + phi2 = phi1 + np.pi/3.0; absc2 = 2.0*np.abs(np.cos(phi2)) + phi3 = phi2 + np.pi/3.0; absc3 = 2.0*np.abs(np.cos(phi3)) + left = ( absc1**m + absc2**m + absc3**m ) + r = (0.5*left)**(1.0/m)*np.sqrt(3.0*I2)/eqStress + + if not Jac: + return (r - 1.0).ravel() + else: + dfdl = r/left/m + jm = r*math_ln(0.5*left)*(-1.0/m/m) + dfdl*0.5*( + absc1**m*math_ln(absc1) + absc2**m*math_ln(absc2) + absc3**m*math_ln(absc3) ) + + da,db,dc = (2.0*A-B-C)/18.0, (2.0*B-C-A)/18.0, (2.0*C-A-B)/18.0 + if dim == 2: + dI2dx = np.array([da, db, dc, H])/1.5*dXdx + dI3dx = np.array([ da*(B-C) + (H**2-G**2)/2.0, + db*(C-A) + (F**2-H**2)/2.0, + dc*(A-B) + (G**2-F**2)/2.0, + (G*F + (A-B))*H ])/3.0*dXdx + else: + dI2dx = np.array([da, db, dc, F,G,H])/1.5*dXdx + dI3dx = np.array([ da*(B-C) + (H**2-G**2)/2.0, + db*(C-A) + (F**2-H**2)/2.0, + dc*(A-B) + (G**2-F**2)/2.0, + (H*G*3.0 + (B-C))*F, + (F*H*3.0 + (C-A))*G, + (G*F*3.0 + (A-B))*H ])/3.0*dXdx + darccos = -1.0/np.sqrt(1.0 - I3**2/I2**3) + + dfdcos = lambda phi : dfdl*m*(2.0*abs(np.cos(phi)))**(m-1.0)*np.sign(np.cos(phi))*(-np.sin(phi)/1.5) + + dfdthe= (dfdcos(phi1) + dfdcos(phi2) + dfdcos(phi3)) + dfdI2, dfdI3 = dfdthe*darccos*I3*(-1.5)*I2**(-2.5)+r/2.0/I2, dfdthe*darccos*I2**(-1.5) + + if mFix[0]: return np.vstack((dfdI2*dI2dx + dfdI3*dI3dx)).T + else: return np.vstack((dfdI2*dI2dx + dfdI3*dI3dx, jm)).T + +def BBC2000(eqStress, paras, sigmas, mFix, criteria, dim, Jac=False): + """ + BBC2000 yield criterion + + the fitted parameters are + d,e,f,g, b,c,a, k; k is optional + criteria are invalid input + """ + d,e,f,g, b,c,a= paras[0:7] + if mFix[0]: k = mFix[1] + else: k = paras[-1] + + s11,s22,s12 = sigmas[0], sigmas[1], sigmas[3] + k2 = 2.0*k; k1 = k - 1.0 + M,N,P,Q,R = d+e, e+f, (d-e)/2.0, (e-f)/2.0, g**2 + Gamma = M*s11 + N*s22 + Psi = ( (P*s11 + Q*s22)**2 + s12**2*R )**0.5 + + l1, l2, l3 = b*Gamma + c*Psi, b*Gamma - c*Psi, 2.0*c*Psi + l1s,l2s,l3s = l1**2, l2**2, l3**2 + + left = a*l1s**k + a*l2s**k + (1-a)*l3s**k + r = left**(1.0/k2)/eqStress + if not Jac: + return (r - 1.0).ravel() + else: + drdl,drdk = r/left/k2, r*math_ln(left)*(-1.0/k2/k) + dldl1,dldl2,dldl3 = a*k2*(l1s**k1)*l1, a*k2*(l2s**k1)*l2, (1-a)*k2*(l3s**k1)*l3 + dldGama, dldPsi = (dldl1 + dldl2)*b, (dldl1 - dldl2 + 2.0*dldl3)*c + temp = (P*s11 + Q*s22)/Psi + dPsidP, dPsidQ, dPsidR = temp*s11, temp*s22, 0.5*s12**2/Psi + dlda = l1s**k + l2s**k - l3s**k + dldb = dldl1*Gamma + dldl2*Gamma + dldc = dldl1*Psi - dldl2*Psi + dldl3*2.0*Psi + dldk = a*math_ln(l1s)*l1s**k + a*math_ln(l2s)*l2s**k + (1-a)*math_ln(l3s)*l3s**k + + J = drdl*np.array([dldGama*s11+dldPsi*dPsidP*0.5, dldGama*(s11+s22)+dldPsi*(-dPsidP+dPsidQ)*0.5, #jd,je + dldGama*s22-dldPsi*dPsidQ*0.5, dldPsi*dPsidR*2.0*g, #jf,jg + dldb, dldc, dlda]) #jb,jc,ja + if mFix[0]: return np.vstack(J).T + else: return np.vstack((J, drdl*dldk + drdk)).T + + +def BBC2003(eqStress, paras, sigmas, mFix, criteria, dim, Jac=False): + """ + BBC2003 yield criterion + + the fitted parameters are + M,N,P,Q,R,S,T,a, k; k is optional + criteria are invalid input + """ + M,N,P,Q,R,S,T,a = paras[0:8] + if mFix[0]: k = mFix[1] + else: k = paras[-1] + + s11,s22,s12 = sigmas[0], sigmas[1], sigmas[3] + k2 = 2.0*k; k1 = k - 1.0 + Gamma = 0.5 * (s11 + M*s22) + Psi = ( 0.25*(N*s11 - P*s22)**2 + Q*Q*s12**2 )**0.5 + Lambda = ( 0.25*(R*s11 - S*s22)**2 + T*T*s12**2 )**0.5 + + l1, l2, l3 = Gamma + Psi, Gamma - Psi, 2.0*Lambda + l1s,l2s,l3s = l1**2, l2**2, l3**2 + left = a*l1s**k + a*l2s**k + (1-a)*l3s**k + r = left**(1.0/k2)/eqStress + if not Jac: + return (r - 1.0).ravel() + else: + drdl,drdk = r/left/k2, r*math_ln(left)*(-1.0/k2/k) + dldl1,dldl2,dldl3 = a*k2*(l1s**k1)*l1, a*k2*(l2s**k1)*l2, (1-a)*k2*(l3s**k1)*l3 + + dldGamma, dldPsi, dldLambda = dldl1+dldl2, dldl1-dldl2, 2.0*dldl3 + temp = 0.25/Psi*(N*s11 - P*s22) + dPsidN, dPsidP, dPsidQ = s11*temp, -s22*temp, Q*s12**2/Psi + temp = 0.25/Lambda*(R*s11 - S*s22) + dLambdadR, dLambdadS, dLambdadT = s11*temp, -s22*temp, T*s12**2/Psi + dldk = a*math_ln(l1s)*l1s**k + a*math_ln(l2s)*l2s**k + (1-a)*math_ln(l3s)*l3s**k + + J = drdl * np.array([dldGamma*s22*0.5, #jM + dldPsi*dPsidN, dldPsi*dPsidP, dldPsi*dPsidQ, #jN, jP, jQ + dldLambda*dLambdadR, dldLambda*dLambdadS, dldLambda*dLambdadT, #jR, jS, jT + l1s**k + l2s**k - l3s**k ]) #ja + + if mFix[0]: return np.vstack(J).T + else : return np.vstack((J, drdl*dldk+drdk)).T + +def BBC2005(eqStress, paras, sigmas, mFix, criteria, dim, Jac=False): + """ + BBC2005 yield criterion + + the fitted parameters are + a, b, L ,M, N, P, Q, R, k k are optional + criteria is invalid input + """ + a,b,L, M, N, P, Q, R = paras[0:8] + if mFix[0]: k = mFix[1] + else: k = paras[-1] + + s11 = sigmas[0]; s22 = sigmas[1]; s12 = sigmas[3] + k2 = 2.0*k + Gamma = L*s11 + M*s22 + Lambda = ( (N*s11 - P*s22)**2 + s12**2 )**0.5 + Psi = ( (Q*s11 - R*s22)**2 + s12**2 )**0.5 + + l1 = Lambda + Gamma; l2 = Lambda - Gamma; l3 = Lambda + Psi; l4 = Lambda - Psi + l1s = l1**2; l2s = l2**2; l3s = l3**2; l4s = l4**2 + left = a*l1s**k + a*l2s**k + b*l3s**k + b*l4s**k + sBar = left**(1.0/k2); r = sBar/eqStress - 1.0 + if not Jac: + return r.ravel() + else: + ln = lambda x : np.log(x + 1.0e-32) + expo = 0.5/k; k1 = k-1.0 + + dsBardl = expo*sBar/left/eqStress + dsBarde = sBar*ln(left); dedk = expo/(-k) + dldl1 = a*k*(l1s**k1)*(2.0*l1) + dldl2 = a*k*(l2s**k1)*(2.0*l2) + dldl3 = b*k*(l3s**k1)*(2.0*l3) + dldl4 = b*k*(l4s**k1)*(2.0*l4) + + dldLambda = dldl1 + dldl2 + dldl3 + dldl4 + dldGama = dldl1 - dldl2 + dldPsi = dldl3 - dldl4 + temp = (N*s11 - P*s22)/Lambda + dLambdadN = s11*temp; dLambdadP = -s22*temp + temp = (Q*s11 - R*s22)/Psi + dPsidQ = s11*temp; dPsidR = -s22*temp + dldk = a*ln(l1s)*l1s**k + a*ln(l2s)*l2s**k + b*ln(l3s)*l3s**k + b*ln(l4s)*l4s**k + + J = dsBardl * np.array( [ + l1s**k+l2s**k, l3s**k+l4s**k,dldGama*s11,dldGama*s22,dldLambda*dLambdadN, + dldLambda*dLambdadP, dldPsi*dPsidQ, dldPsi*dPsidR]) + + if mFix[0]: return np.vstack(J).T + else : return np.vstack(J, dldk+dsBarde*dedk).T + +def Yld2000(eqStress, paras, sigmas, mFix, criteria, dim, Jac=False): + """ + Yld2000 yield criterion + + C: c11,c22,c66 c12=c21=1.0 JAC NOT PASS + D: d11,d12,d21,d22,d66 + """ + C,D = paras[0:3], paras[3:8] + if mFix[0]: m = mFix[1] + else: m = paras[-1] + + s11, s22, s12 = sigmas[0],sigmas[1],sigmas[3] + X = np.array([ 2.0*C[0]*s11-C[0]*s22, 2.0*C[1]*s22-C[1]*s11, 3.0*C[2]*s12 ])/3.0 # a1,a2,a7 + Y = np.array([ (8.0*D[2]-2.0*D[0]-2.0*D[3]+2.0*D[1])*s11 + (4.0*D[3]-4.0*D[1]-4.0*D[2]+ D[0])*s22, + (4.0*D[0]-4.0*D[2]-4.0*D[1]+ D[3])*s11 + (8.0*D[1]-2.0*D[3]-2.0*D[0]+2.0*D[2])*s22, + 9.0*D[4]*s12 ])/9.0 + + def priStrs(s): + temp = np.sqrt( (s[0]-s[1])**2 + 4.0*s[2]**2 ) + return 0.5*(s[0]+s[1] + temp), 0.5*(s[0]+s[1] - temp) + m2 = m/2.0; m21 = m2 - 1.0 + (X1,X2), (Y1,Y2) = priStrs(X), priStrs(Y) # Principal values of X, Y + phi1s, phi21s, phi22s = (X1-X2)**2, (2.0*Y2+Y1)**2, (2.0*Y1+Y2)**2 + phi1, phi21, phi22 = phi1s**m2, phi21s**m2, phi22s**m2 + left = phi1 + phi21 + phi22 + r = (0.5*left)**(1.0/m)/eqStress + + if not Jac: + return (r-1.0).ravel() + else: + drdl, drdm = r/m/left, r*math_ln(0.5*left)*(-1.0/m/m) #/(-m*m) + dldm = ( phi1*math_ln(phi1s) + phi21*math_ln(phi21s) + phi22*math_ln(phi22s) )*0.5 + zero = np.zeros_like(s11); num = len(s11) + def dPrincipalds(X): + """Derivative of principla with respect to stress""" + temp = 1.0/np.sqrt( (X[0]-X[1])**2 + 4.0*X[2]**2 ) + dP1dsi = 0.5*np.array([ 1.0+temp*(X[0]-X[1]), 1.0-temp*(X[0]-X[1]), temp*4.0*X[2]]) + dP2dsi = 0.5*np.array([ 1.0-temp*(X[0]-X[1]), 1.0+temp*(X[0]-X[1]), -temp*4.0*X[2]]) + return np.array([dP1dsi, dP2dsi]) + + dXdXi, dYdYi = dPrincipalds(X), dPrincipalds(Y) + dXidC = np.array([ [ 2.0*s11-s22, zero, zero ], #dX11dC + [ zero, 2.0*s22-s11, zero ], #dX22dC + [ zero, zero, 3.0*s12 ] ])/3.0 #dX12dC + dYidD = np.array([ [ -2.0*s11+ s22, 2.0*s11-4.0*s22, 8.0*s11-4.0*s22, -2.0*s11+4.0*s22, zero ], #dY11dD + [ 4.0*s11-2.0*s22, -4.0*s11+8.0*s22, -4.0*s11+2.0*s22, s11-2.0*s22, zero ], #dY22dD + [ zero, zero, zero, zero, 9.0*s12 ] ])/9.0 #dY12dD + + dXdC=np.array([np.dot(dXdXi[:,:,i], dXidC[:,:,i]).T for i in range(num)]).T + dYdD=np.array([np.dot(dYdYi[:,:,i], dYidD[:,:,i]).T for i in range(num)]).T + + dldX = m*np.array([ phi1s**m21*(X1-X2), phi1s**m21*(X2-X1)]) + dldY = m*np.array([phi21s**m21*(2.0*Y2+Y1) + 2.0*phi22s**m21*(2.0*Y1+Y2), \ + phi22s**m21*(2.0*Y1+Y2) + 2.0*phi21s**m21*(2.0*Y2+Y1) ]) + jC = drdl*np.array([np.dot(dldX[:,i], dXdC[:,:,i]) for i in range(num)]).T + jD = drdl*np.array([np.dot(dldY[:,i], dYdD[:,:,i]) for i in range(num)]).T + + jm = drdl*dldm + drdm + if mFix[0]: return np.vstack((jC,jD)).T + else: return np.vstack((jC,jD,jm)).T + +def Yld200418p(eqStress, paras, sigmas, mFix, criteria, dim, Jac=False): + """ + Yld2004-18p yield criterion + + the fitted parameters are + C: c12,c21,c23,c32,c31,c13,c44,c55,c66; D: d12,d21,d23,d32,d31,d13,d44,d55,d66 for 3D + C: c12,c21,c23,c32,c31,c13,c44; D: d12,d21,d23,d32,d31,d13,d44 for 2D + and m, m are optional + criteria is ignored + """ + if dim == 2: C,D = np.append(paras[0:7],[0.0,0.0]), np.append(paras[7:14],[0.0,0.0]) + else: C,D = paras[0:9], paras[9:18] + if mFix[0]: m = mFix[1] + else: m = paras[-1] + + sv = (sigmas[0] + sigmas[1] + sigmas[2])/3.0 + sdev = np.vstack((sigmas[0:3]-sv,sigmas[3:6])) + ys = lambda sdev, C: np.array([-C[0]*sdev[1]-C[5]*sdev[2], -C[1]*sdev[0]-C[2]*sdev[2], + -C[4]*sdev[0]-C[3]*sdev[1], C[6]*sdev[3], C[7]*sdev[4], C[8]*sdev[5]]) + p,q = ys(sdev, C), ys(sdev, D) + pLambdas, qLambdas = principalStress(p), principalStress(q) # no sort + + m2 = m/2.0; x3 = range(3); num = len(sv) + PiQj = np.array([(pLambdas[i,:]-qLambdas[j,:]) for i in x3 for j in x3]) + QiPj = np.array([(qLambdas[i,:]-pLambdas[j,:]) for i in x3 for j in x3]).reshape(3,3,num) + PiQjs = PiQj**2 + left = np.sum(PiQjs**m2,axis=0) + r = (0.25*left)**(1.0/m)/eqStress + + if not Jac: + return (r - 1.0).ravel() + else: + drdl, drdm = r/m/left, r*math_ln(0.25*left)*(-1.0/m/m) + dldm = np.sum(PiQjs**m2*math_ln(PiQjs),axis=0)*0.5 + dPdc, dQdd = principalStrs_Der(p, sdev, dim), principalStrs_Der(q, sdev, dim) + PiQjs3d = ( PiQjs**(m2-1.0) ).reshape(3,3,num) + dldP = -m*np.array([np.diag(np.dot(PiQjs3d[:,:,i], QiPj [:,:,i])) for i in range(num)]).T + dldQ = m*np.array([np.diag(np.dot(QiPj [:,:,i], PiQjs3d[:,:,i])) for i in range(num)]).T + + jm = drdl*dldm + drdm + jc = drdl*np.sum([dldP[i]*dPdc[i] for i in x3],axis=0) + jd = drdl*np.sum([dldQ[i]*dQdd[i] for i in x3],axis=0) + + if mFix[0]: return np.vstack((jc,jd)).T + else: return np.vstack((jc,jd,jm)).T + +def KarafillisBoyce(eqStress, paras, sigmas, mFix, criteria, dim, Jac=False): + """ + Karafillis-Boyce + + the fitted parameters are + c11,c12,c13,c14,c15,c16,c,m for 3D + c11,c12,c13,c14,c,m for plane stress + 0 1 and self.dimen == 2: + return fitCriteria[self.name]['labels'][1] + else: + return fitCriteria[self.name]['labels'][0] + + def report_name(self): + return fitCriteria[self.name]['name'] + + def fit(self,stress): + global fitResults; fitErrors; fitResidual + if options.exponent > 0.0: nExponent = options.exponent + else: nExponent = 0 + nameCriterion = self.name.lower() + criteria = Criteria(nameCriterion,self.uniaxial,self.expo, self.dimen) + bounds = fitCriteria[nameCriterion]['bound'][dDim] # Default bounds, no bound + guess0 = Guess # Default initial guess, depends on bounds + + if fitResults == []: + initialguess = guess0 + else: + initialguess = np.array(fitResults[-1]) + + ydata = np.zeros(np.shape(stress)[1]) + try: + popt, pcov, infodict, errmsg, ierr = \ + leastsqBound (criteria.fun, initialguess, args=(ydata,stress), + bounds=bounds, Dfun=criteria.jac, full_output=True) + if ierr not in [1, 2, 3, 4]: + raise RuntimeError("Optimal parameters not found: "+errmsg) + else: + residual = criteria.fun(popt, ydata, stress) + fitResidual.append(np.linalg.norm(residual)/np.sqrt(len(residual))) + if (len(ydata) > len(initialguess)) and pcov is not None: + s_sq = (criteria.fun(popt, *(ydata,stress))**2).sum()/(len(ydata)-len(initialguess)) + pcov = pcov * s_sq + perr = np.sqrt(np.diag(pcov)) + fitResults.append(popt.tolist()) + fitErrors .append(perr.tolist()) + + popt = np.concatenate((np.array(popt), np.repeat(options.exponent,nExponent))) + perr = np.concatenate((np.array(perr), np.repeat(0.0,nExponent))) + + damask.util.croak('Needed {} function calls for fitting'.format(infodict['nfev'])) + except Exception as detail: + damask.util.croak(detail) + pass + return popt + +#--------------------------------------------------------------------------------------------------- +class myThread (threading.Thread): + """Runner""" + + def __init__(self, threadID): + threading.Thread.__init__(self) + self.threadID = threadID + def run(self): + semaphore.acquire() + conv=converged() + semaphore.release() + while not conv: + doSim(self.name) + semaphore.acquire() + conv=converged() + semaphore.release() + +def doSim(thread): + semaphore.acquire() + global myLoad + loadNo=loadcaseNo() + if not os.path.isfile('%s.load'%loadNo): + damask.util.croak('Generating load case for simulation %s (%s)'%(loadNo,thread)) + f=open('%s.load'%loadNo,'w') + f.write(myLoad.getLoadcase(loadNo)) + f.close() + semaphore.release() + else: semaphore.release() + +# if spectralOut does not exist, run simulation + semaphore.acquire() + if not os.path.isfile('%s_%i.spectralOut'%(options.geometry,loadNo)): + damask.util.croak('Starting simulation %i (%s)'%(loadNo,thread)) + semaphore.release() + damask.util.execute('DAMASK_spectral -g %s -l %i'%(options.geometry,loadNo)) + else: semaphore.release() + +# reading values from ASCII file + semaphore.acquire() + damask.util.croak('Reading values from simulation %i (%s)'%(loadNo,thread)) + semaphore.release() + refFile = '%s_%i.yield'%(options.geometry,loadNo) + yieldStress = np.empty((6),'d') + if not os.path.isfile(refFile): + validity = False + else: + validity = True + yieldData = np.loadtxt(refFile) + stress = yieldData[:3] + yieldStress = t33toSym6(stress) +# do the actual fitting procedure and write results to file + semaphore.acquire() + global stressAll + f=open(options.geometry+'_'+options.criterion+'_'+str(time.time())+'.txt','w') + f.write(' '.join([options.fitting]+myFit.report_labels())+'\n') + try: + if validity: + stressAll=np.append(stressAll, yieldStress/stressUnit) + f.write(' '.join(map(str,myFit.fit(stressAll.reshape(len(stressAll)//6,6).transpose())))+'\n') + except Exception: + damask.util.croak('Could not fit results of simulation (%s)'%thread) + semaphore.release() + return + damask.util.croak('\n') + semaphore.release() + +def loadcaseNo(): + global N_simulations + N_simulations+=1 + return N_simulations + +def converged(): + global N_simulations; fitResidual + + if N_simulations < options.max: + if len(fitResidual) > 5 and N_simulations >= options.min: + residualList = np.array(fitResidual[len(fitResidual)-5:]) + if np.std(residualList)/np.max(residualList) < 0.05: + return True + return False + else: + return True + +# -------------------------------------------------------------------- +# MAIN +# -------------------------------------------------------------------- + +parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """ +Performs calculations with various loads on given geometry file and fits yield surface. + +""", version = scriptID) + +# maybe make an option to specifiy if 2D/3D fitting should be done? + +parser.add_option('-l','--load' , dest='load', type='float', nargs=3, + help='load: final strain; increments; time %default', metavar='float int float') +parser.add_option('-g','--geometry', dest='geometry', type='string', + help='name of the geometry file [%default]', metavar='string') +parser.add_option('-c','--criterion', dest='criterion', choices=fitCriteria.keys(), + help='criterion for stopping simulations [%default]', metavar='string') +parser.add_option('-f','--fitting', dest='fitting', choices=thresholdParameter, + help='yield criterion [%default]', metavar='string') +parser.add_option('-y','--yieldvalue', dest='yieldValue', type='float', + help='yield points %default', metavar='float') +parser.add_option('--min', dest='min', type='int', + help='minimum number of simulations [%default]', metavar='int') +parser.add_option('--max', dest='max', type='int', + help='maximum number of iterations [%default]', metavar='int') +parser.add_option('-t','--threads', dest='threads', type='int', + help='number of parallel executions [%default]', metavar='int') +parser.add_option('-b','--bound', dest='bounds', type='float', nargs=2, + help='yield points: start; end; count %default', metavar='float float') +parser.add_option('-d','--dimension', dest='dimension', type='choice', choices=['2','3'], + help='dimension of the virtual test [%default]', metavar='int') +parser.add_option('-e', '--exponent', dest='exponent', type='float', + help='exponent of non-quadratic criteria', metavar='int') +parser.add_option('-u', '--uniaxial', dest='eqStress', type='float', + help='Equivalent stress', metavar='float') +parser.add_option('--flag', dest='flag', type='string', + help='the variable based which the yield will be judged, totalStrain, plasticStrain or plasticWork', metavar='string') + +parser.set_defaults(min = 12, + max = 20, + threads = 4, + yieldValue = 0.002, + load = (0.010,100,100.0), + criterion = 'vonmises', + fitting = 'totalshear', + geometry = '20grains16x16x16', + bounds = None, + dimension = '3', + exponent = -1.0, + flag = 'totalStrain', + ) + +options = parser.parse_args()[0] + +if options.threads < 1: + parser.error('invalid number of threads {}'.format(options.threads)) +if options.min < 0: + parser.error('invalid minimum number of simulations {}'.format(options.min)) +if options.max < options.min: + parser.error('invalid maximum number of simulations (below minimum)') + +for check in [options.geometry+'.geom','numerics.config','material.config']: + if not os.path.isfile(check): + damask.util.croak('"{}" file not found'.format(check)) + +options.dimension = int(options.dimension) + +stressUnit = 1.0e9 if options.criterion == 'hill1948' else 1.0e6 + + +if options.dimension not in fitCriteria[options.criterion]['dimen']: + parser.error('invalid dimension for selected criterion') + +if options.criterion not in ['vonmises','tresca','drucker','hill1948'] and options.eqStress is None: + parser.error('please specify an equivalent stress (e.g. fitting to von Mises)') + +# global variables +fitResults = [] +fitErrors = [] +fitResidual = [] +stressAll= np.zeros(0,'d').reshape(0,0) +N_simulations=0 +Guess = [] +threads=[] +semaphore=threading.Semaphore(1) +dDim = None +myLoad = None +myFit = None + +run = runFit(options.exponent, options.eqStress, options.dimension, options.criterion) + +damask.util.croak('Finished fitting to yield criteria') diff --git a/src/spectral_utilities.f90 b/src/spectral_utilities.f90 index ecd47b594..c5560752e 100755 --- a/src/spectral_utilities.f90 +++ b/src/spectral_utilities.f90 @@ -1140,7 +1140,6 @@ subroutine utilities_calcPlasticity(yieldStress, plasticStrain, eqStress, eqTota strain_plastic_av = strain_plastic_av + strain_plastic Vp_av = Vp_av + Vp(1:3,1:3,i,j,k) - flush(6) enddo; enddo; enddo From 9bbc0d480340fdf99f10ec36bceb18f381440a04 Mon Sep 17 00:00:00 2001 From: Fengbo Han Date: Tue, 8 Aug 2017 17:29:14 +0200 Subject: [PATCH 08/27] implemented fast yield surface fitting with yield stop criteria --- processing/misc/yieldSurfaceFast.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/processing/misc/yieldSurfaceFast.py b/processing/misc/yieldSurfaceFast.py index 7b00124e7..654e90145 100755 --- a/processing/misc/yieldSurfaceFast.py +++ b/processing/misc/yieldSurfaceFast.py @@ -1332,7 +1332,7 @@ parser.add_option('-e', '--exponent', dest='exponent', type='float', parser.add_option('-u', '--uniaxial', dest='eqStress', type='float', help='Equivalent stress', metavar='float') parser.add_option('--flag', dest='flag', type='string', - help='the variable based which the yield will be judged, totalStrain, plasticStrain or plasticWork', metavar='string') + help='yield stop flag, totalStrain, plasticStrain or plasticWork', metavar='string') parser.set_defaults(min = 12, max = 20, From 0750f7fd01ff3697e6f24d6bd18b4bdf0e247647 Mon Sep 17 00:00:00 2001 From: Fengbo Han Date: Thu, 10 Aug 2017 15:40:18 +0200 Subject: [PATCH 09/27] fixed plastic work calculation --- src/DAMASK_spectral.f90 | 10 +++++++++- src/spectral_utilities.f90 | 7 ++++--- 2 files changed, 13 insertions(+), 4 deletions(-) diff --git a/src/DAMASK_spectral.f90 b/src/DAMASK_spectral.f90 index fe7bdf2b3..38684f450 100755 --- a/src/DAMASK_spectral.f90 +++ b/src/DAMASK_spectral.f90 @@ -316,7 +316,15 @@ program DAMASK_spectral end select enddo; enddo close(FILEUNIT) - + + if(yieldStop) then ! initialize variables related to yield stop + yieldStressNew = 0.0_pReal + plasticStrainNew = 0.0_pReal + eqStressNew = 0.0_pReal + eqTotalStrainNew = 0.0_pReal + eqPlasticStrainNew = 0.0_pReal + plasticWorkNew = 0.0_pReal + endif !-------------------------------------------------------------------------------------------------- ! consistency checks and output of load case loadCases(1)%followFormerTrajectory = .false. ! cannot guess along trajectory for first inc of first currentLoadCase diff --git a/src/spectral_utilities.f90 b/src/spectral_utilities.f90 index c5560752e..ecbce38bf 100755 --- a/src/spectral_utilities.f90 +++ b/src/spectral_utilities.f90 @@ -1080,9 +1080,12 @@ subroutine utilities_calcPlasticity(yieldStress, plasticStrain, eqStress, eqTota real(pReal), dimension(15) :: WORK !< previous deformation gradient integer(pInt) :: INFO, i, j, k, l, ierr real(pReal) :: stress, stress_av, strain_total, strain_total_av, strain_plastic, strain_plastic_av, wgtm + real(pReal) :: eqStressOld, eqPlasticStrainOld external :: dgesvd + eqStressOld = eqStress + eqPlasticStrainOld = eqPlasticStrain wgtm = 1.0/real(mesh_NcpElems*mesh_maxNips*homogenization_maxNgrains,pReal) Vp_av = 0.0_pReal V_total_av = 0.0_pReal @@ -1142,20 +1145,18 @@ subroutine utilities_calcPlasticity(yieldStress, plasticStrain, eqStress, eqTota enddo; enddo; enddo - yieldStress = cauchy_av * wgtm call MPI_Allreduce(MPI_IN_PLACE,yieldStress,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) plasticStrain = Vp_av * wgtm call MPI_Allreduce(MPI_IN_PLACE,plasticStrain,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) - plasticWork = plasticWork + 0.5*(eqStress + stress_av * wgtm) * (strain_total_av * wgtm - eqTotalStrain) - call MPI_Allreduce(MPI_IN_PLACE,plasticWork,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) eqStress = stress_av * wgtm call MPI_Allreduce(MPI_IN_PLACE,eqStress,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) eqTotalStrain = strain_total_av * wgtm call MPI_Allreduce(MPI_IN_PLACE,eqTotalStrain,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) eqPlasticStrain = strain_plastic_av * wgtm call MPI_Allreduce(MPI_IN_PLACE,eqPlasticStrain,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) + plasticWork = plasticWork + 0.5*(eqStressOld + eqStress) * (eqPlasticStrain - eqPlasticStrainOld) end subroutine utilities_calcPlasticity From 82758bd90f550c3ad45fa6a3fec762b7feb89e17 Mon Sep 17 00:00:00 2001 From: Fengbo Han Date: Tue, 24 Oct 2017 11:15:34 +0200 Subject: [PATCH 10/27] added yield criterion of facet potential --- processing/misc/yieldSurfaceFast.py | 182 ++++++++++++++++++++++++---- src/math.f90 | 4 +- 2 files changed, 158 insertions(+), 28 deletions(-) diff --git a/processing/misc/yieldSurfaceFast.py b/processing/misc/yieldSurfaceFast.py index 654e90145..df8dcbd28 100755 --- a/processing/misc/yieldSurfaceFast.py +++ b/processing/misc/yieldSurfaceFast.py @@ -6,6 +6,8 @@ import numpy as np from optparse import OptionParser import damask from damask.util import leastsqBound +from scipy.optimize import nnls +import math scriptName = os.path.splitext(os.path.basename(__file__))[0] scriptID = ' '.join([scriptName,damask.version]) @@ -14,23 +16,24 @@ def runFit(exponent, eqStress, dimension, criterion): global threads, myFit, myLoad global fitResidual global Guess, dDim + + if options.criterion!='facet': + dDim = dimension - 3 + nParas = len(fitCriteria[criterion]['bound'][dDim]) + nExpo = fitCriteria[criterion]['nExpo'] - dDim = dimension - 3 - nParas = len(fitCriteria[criterion]['bound'][dDim]) - nExpo = fitCriteria[criterion]['nExpo'] + if exponent > 0.0: # User defined exponents + nParas = nParas-nExpo + fitCriteria[criterion]['bound'][dDim] = fitCriteria[criterion]['bound'][dDim][:nParas] - if exponent > 0.0: # User defined exponents - nParas = nParas-nExpo - fitCriteria[criterion]['bound'][dDim] = fitCriteria[criterion]['bound'][dDim][:nParas] - - for i in range(nParas): - temp = fitCriteria[criterion]['bound'][dDim][i] - if fitCriteria[criterion]['bound'][dDim][i] == (None,None): - Guess.append(1.0) - else: - g = (temp[0]+temp[1])/2.0 - if g == 0: g = temp[1]*0.5 - Guess.append(g) + for i in range(nParas): + temp = fitCriteria[criterion]['bound'][dDim][i] + if fitCriteria[criterion]['bound'][dDim][i] == (None,None): + Guess.append(1.0) + else: + g = (temp[0]+temp[1])/2.0 + if g == 0: g = temp[1]*0.5 + Guess.append(g) myLoad = Loadcase(options.load[0],options.load[1],options.load[2],options.flag,options.yieldValue, nSet = 10, dimension = dimension, vegter = options.criterion=='vegter') @@ -43,9 +46,92 @@ def runFit(exponent, eqStress, dimension, criterion): for t in range(options.threads): threads[t].join() + + if options.criterion=='facet': + doFacetFit() + damask.util.croak('Residuals') damask.util.croak(fitResidual) +def doFacetFit(): + n = options.order + Data = np.zeros((options.numpoints, 10)) + for i in range(options.numpoints): + fileName = options.geometry + '_' + str(i+1) + '.yield' + while os.path.exists(fileName): + break + data_i = np.loadtxt(fileName) + + sv = (data_i[0,0] + data_i[1,1] + data_i[2,2])/3.0 + + #convert stress and strain form the 6D to 5D space + S1 = math.sqrt(2.0)*(data_i[0,0] - data_i[1,1])/2.0 + S2 = math.sqrt(6.0)*(data_i[0,0] + data_i[1,1] - 2.0*sv)/2.0 + S3 = math.sqrt(2.0)*data_i[1,2] + S4 = math.sqrt(2.0)*data_i[2,0] + S5 = math.sqrt(2.0)*data_i[0,1] + + E1 = math.sqrt(2.0)*(data_i[3,0]-data_i[4,1])/2.0 + E2 = math.sqrt(6.0)*(data_i[3,0]+data_i[4,1])/2.0 + E3 = math.sqrt(2.0)*data_i[4,2] + E4 = math.sqrt(2.0)*data_i[5,0] + E5 = math.sqrt(2.0)*data_i[3,1] + + Data[i,:] = [E1,E2,E3,E4,E5,S1,S2,S3,S4,S5] + + Data[:,5:] = Data[:,5:] / 100000000.0 + + path=os.path.join(os.getcwd(),'final.mmm') + np.savetxt(path, Data, header='', comments='', fmt='% 15.10f') + + if options.dimension == 2: + reducedIndices = [0,1,4,5,6,9] + strainRateIndices = [0,1,4] + stressIndices = [5,6,9] + elif options.dimension == 3: + reducedIndices = [i for i in range(10)] + strainRateIndices = [0,1,2,3,4] + stressIndices = [5,6,7,8,9] + + numDirections = Data.shape[0] + Indices = np.arange(numDirections) + sdPairs = Data[:,reducedIndices][Indices,:] + numPairs = sdPairs.shape[0] + dimensionality = sdPairs.shape[1] / 2 + ds = sdPairs[:,0:dimensionality] + s = sdPairs[:,dimensionality::] + + A = np.zeros((numPairs, numPairs)) + B = np.ones((numPairs,)) + for i in range(numPairs): + for j in range(numPairs): + lamb = 1.0 + s_i = s[i,:] + ds_j = ds[j,:] + A[i,j] = lamb * (np.dot(s_i.ravel(), ds_j.ravel()) ** n) + + lambdas, residuals = nnls(A, B) + nonZeroTerms = np.logical_not(np.isclose(lambdas, 0.)) + numNonZeroTerms = np.sum(nonZeroTerms) + dataOut = np.zeros((numNonZeroTerms, 6)) + + if options.dimension == 2: + dataOut[:,0] = lambdas[nonZeroTerms] + dataOut[:,1] = ds[nonZeroTerms,:][:,0] + dataOut[:,2] = ds[nonZeroTerms,:][:,1] + dataOut[:,5] = ds[nonZeroTerms,:][:,2] + elif options.dimension == 3: + dataOut[:,0] = lambdas[nonZeroTerms] + dataOut[:,1] = ds[nonZeroTerms,:][:,0] + dataOut[:,2] = ds[nonZeroTerms,:][:,1] + dataOut[:,3] = ds[nonZeroTerms,:][:,2] + dataOut[:,4] = ds[nonZeroTerms,:][:,3] + dataOut[:,5] = ds[nonZeroTerms,:][:,4] + + headerText = 'facet\n 1 \n F \n {0:<3d} \n {1:<3d} '.format(n, numNonZeroTerms) + path=os.path.join(os.getcwd(),'facet_o{0}.fac'.format(n)) + np.savetxt(path, dataOut, header=headerText, comments='', fmt='% 15.10f') + def principalStresses(sigmas): """ Computes principal stresses (i.e. eigenvalues) for a set of Cauchy stresses. @@ -1014,7 +1100,13 @@ fitCriteria = { 'vegter' :{'name': 'Vegter', 'labels': 'a,b,c,d,e,f,g,h', 'dimen': [2], - }, + }, + 'facet' :{'name': 'Facet', + 'nExpo': None, + 'bound': [(None,None)], + 'labels': 'lambdas', + 'dimen': [2,3], + }, } thresholdParameter = ['totalshear','equivalentStrain'] @@ -1225,7 +1317,10 @@ class myThread (threading.Thread): conv=converged() semaphore.release() while not conv: - doSim(self.name) + if options.criterion=='facet': + doSimForFacet(self.name) + else: + doSim(self.name) semaphore.acquire() conv=converged() semaphore.release() @@ -1278,6 +1373,26 @@ def doSim(thread): return damask.util.croak('\n') semaphore.release() + +def doSimForFacet(thread): + semaphore.acquire() + global myLoad + loadNo=loadcaseNo() + if not os.path.isfile('%s.load'%loadNo): + damask.util.croak('Generating load case for simulation %s (%s)'%(loadNo,thread)) + f=open('%s.load'%loadNo,'w') + f.write(myLoad.getLoadcase(loadNo)) + f.close() + semaphore.release() + else: semaphore.release() + +# if spectralOut does not exist, run simulation + semaphore.acquire() + if not os.path.isfile('%s_%i.spectralOut'%(options.geometry,loadNo)): + damask.util.croak('Starting simulation %i (%s)'%(loadNo,thread)) + semaphore.release() + damask.util.execute('DAMASK_spectral -g %s -l %i'%(options.geometry,loadNo)) + else: semaphore.release() def loadcaseNo(): global N_simulations @@ -1286,15 +1401,21 @@ def loadcaseNo(): def converged(): global N_simulations; fitResidual - - if N_simulations < options.max: - if len(fitResidual) > 5 and N_simulations >= options.min: - residualList = np.array(fitResidual[len(fitResidual)-5:]) - if np.std(residualList)/np.max(residualList) < 0.05: - return True - return False + + if options.criterion=='facet': + if N_simulations == options.numpoints: + return True + else: + return False else: - return True + if N_simulations < options.max: + if len(fitResidual) > 5 and N_simulations >= options.min: + residualList = np.array(fitResidual[len(fitResidual)-5:]) + if np.std(residualList)/np.max(residualList) < 0.05: + return True + return False + else: + return True # -------------------------------------------------------------------- # MAIN @@ -1333,6 +1454,10 @@ parser.add_option('-u', '--uniaxial', dest='eqStress', type='float', help='Equivalent stress', metavar='float') parser.add_option('--flag', dest='flag', type='string', help='yield stop flag, totalStrain, plasticStrain or plasticWork', metavar='string') +parser.add_option('--numpoints', dest='numpoints', type='int', + help='number of yield points to fit facet potential [%default]', metavar='int') +parser.add_option('--order', dest='order', type='int', + help='order of facet potential [%default]', metavar='int') parser.set_defaults(min = 12, max = 20, @@ -1346,6 +1471,8 @@ parser.set_defaults(min = 12, dimension = '3', exponent = -1.0, flag = 'totalStrain', + numpoints = 100, + order = 8 ) options = parser.parse_args()[0] @@ -1385,6 +1512,9 @@ dDim = None myLoad = None myFit = None -run = runFit(options.exponent, options.eqStress, options.dimension, options.criterion) +if options.criterion == 'facet': + run = runFit(options.exponent, options.eqStress, options.dimension, options.criterion) +else: + run = runFit(options.exponent, options.eqStress, options.dimension, options.criterion) damask.util.croak('Finished fitting to yield criteria') diff --git a/src/math.f90 b/src/math.f90 index 937f00d7a..11ecd6e0d 100755 --- a/src/math.f90 +++ b/src/math.f90 @@ -982,7 +982,7 @@ real(pReal) pure function math_detSym33(m) real(pReal), dimension(3,3), intent(in) :: m math_detSym33 = -(m(1,1)*m(2,3)**2_pInt + m(2,2)*m(1,3)**2_pInt + m(3,3)*m(1,2)**2_pInt) & - + m(1,1)*m(2,2)*m(3,3) - 2.0_pReal * m(1,2)*m(1,3)*m(2,3) + + m(1,1)*m(2,2)*m(3,3) + 2.0_pReal * m(1,2)*m(1,3)*m(2,3) end function math_detSym33 @@ -1994,7 +1994,7 @@ function math_eigenvectorBasisSym(m) do i=1_pInt, size(m,1) math_eigenvectorBasisSym = math_eigenvectorBasisSym & - + sqrt(values(i)) * math_tensorproduct(vectors(:,i),vectors(:,i)) + + sqrt(values(i)) * math_tensorproduct(vectors(:,i),vectors(:,i)) enddo end function math_eigenvectorBasisSym From d77789b89f3e425a7c1abcfdc0716deaffc3d992 Mon Sep 17 00:00:00 2001 From: Fengbo Han Date: Tue, 24 Oct 2017 12:01:10 +0200 Subject: [PATCH 11/27] minor change --- processing/misc/yieldSurfaceFast.py | 23 ++++++++++------------- 1 file changed, 10 insertions(+), 13 deletions(-) diff --git a/processing/misc/yieldSurfaceFast.py b/processing/misc/yieldSurfaceFast.py index df8dcbd28..7be33041f 100755 --- a/processing/misc/yieldSurfaceFast.py +++ b/processing/misc/yieldSurfaceFast.py @@ -7,7 +7,6 @@ from optparse import OptionParser import damask from damask.util import leastsqBound from scipy.optimize import nnls -import math scriptName = os.path.splitext(os.path.basename(__file__))[0] scriptID = ' '.join([scriptName,damask.version]) @@ -58,24 +57,22 @@ def doFacetFit(): Data = np.zeros((options.numpoints, 10)) for i in range(options.numpoints): fileName = options.geometry + '_' + str(i+1) + '.yield' - while os.path.exists(fileName): - break data_i = np.loadtxt(fileName) sv = (data_i[0,0] + data_i[1,1] + data_i[2,2])/3.0 #convert stress and strain form the 6D to 5D space - S1 = math.sqrt(2.0)*(data_i[0,0] - data_i[1,1])/2.0 - S2 = math.sqrt(6.0)*(data_i[0,0] + data_i[1,1] - 2.0*sv)/2.0 - S3 = math.sqrt(2.0)*data_i[1,2] - S4 = math.sqrt(2.0)*data_i[2,0] - S5 = math.sqrt(2.0)*data_i[0,1] + S1 = np.sqrt(2.0)*(data_i[0,0] - data_i[1,1])/2.0 + S2 = np.sqrt(6.0)*(data_i[0,0] + data_i[1,1] - 2.0*sv)/2.0 + S3 = np.sqrt(2.0)*data_i[1,2] + S4 = np.sqrt(2.0)*data_i[2,0] + S5 = np.sqrt(2.0)*data_i[0,1] - E1 = math.sqrt(2.0)*(data_i[3,0]-data_i[4,1])/2.0 - E2 = math.sqrt(6.0)*(data_i[3,0]+data_i[4,1])/2.0 - E3 = math.sqrt(2.0)*data_i[4,2] - E4 = math.sqrt(2.0)*data_i[5,0] - E5 = math.sqrt(2.0)*data_i[3,1] + E1 = np.sqrt(2.0)*(data_i[3,0]-data_i[4,1])/2.0 + E2 = np.sqrt(6.0)*(data_i[3,0]+data_i[4,1])/2.0 + E3 = np.sqrt(2.0)*data_i[4,2] + E4 = np.sqrt(2.0)*data_i[5,0] + E5 = np.sqrt(2.0)*data_i[3,1] Data[i,:] = [E1,E2,E3,E4,E5,S1,S2,S3,S4,S5] From 2b69cb8ea947d1cc9a70fe0a52bc47d5dccda681 Mon Sep 17 00:00:00 2001 From: Fengbo Han Date: Tue, 24 Oct 2017 13:23:44 +0200 Subject: [PATCH 12/27] removed unnecessary liens --- processing/misc/yieldSurfaceFast.py | 4 ---- 1 file changed, 4 deletions(-) diff --git a/processing/misc/yieldSurfaceFast.py b/processing/misc/yieldSurfaceFast.py index 7be33041f..c58dca733 100755 --- a/processing/misc/yieldSurfaceFast.py +++ b/processing/misc/yieldSurfaceFast.py @@ -83,12 +83,8 @@ def doFacetFit(): if options.dimension == 2: reducedIndices = [0,1,4,5,6,9] - strainRateIndices = [0,1,4] - stressIndices = [5,6,9] elif options.dimension == 3: reducedIndices = [i for i in range(10)] - strainRateIndices = [0,1,2,3,4] - stressIndices = [5,6,7,8,9] numDirections = Data.shape[0] Indices = np.arange(numDirections) From 2b5a536458ca88c81c7020fb826943604609875c Mon Sep 17 00:00:00 2001 From: Fengbo Han Date: Wed, 22 Nov 2017 08:52:48 +0100 Subject: [PATCH 13/27] calculate stress and strain from the average PK stress and average deformation gradient of the whole RVE --- src/DAMASK_spectral.f90 | 2 +- src/spectral_utilities.f90 | 136 ++++++++++++++++++------------------- 2 files changed, 67 insertions(+), 71 deletions(-) mode change 100644 => 100755 src/DAMASK_spectral.f90 diff --git a/src/DAMASK_spectral.f90 b/src/DAMASK_spectral.f90 old mode 100644 new mode 100755 index 05327a985..e51f7636a --- a/src/DAMASK_spectral.f90 +++ b/src/DAMASK_spectral.f90 @@ -708,7 +708,7 @@ program DAMASK_spectral plasticWorkOld = plasticWorkNew call utilities_calcPlasticity(yieldStressNew, plasticStrainNew, eqStressNew, eqTotalStrainNew, & - eqPlasticStrainNew, plasticWorkNew) + eqPlasticStrainNew, plasticWorkNew, loadCases(currentLoadCase)%rotation) if(stopFlag == 'totalStrain') then if(eqTotalStrainNew > yieldStopValue) then diff --git a/src/spectral_utilities.f90 b/src/spectral_utilities.f90 index ee7e27e4c..742068c9a 100755 --- a/src/spectral_utilities.f90 +++ b/src/spectral_utilities.f90 @@ -1048,9 +1048,9 @@ end subroutine utilities_constitutiveResponse !-------------------------------------------------------------------------------------------------- !> @brief calculates yield stress, plastic strain, total strain and their equivalent values !-------------------------------------------------------------------------------------------------- -subroutine utilities_calcPlasticity(yieldStress, plasticStrain, eqStress, eqTotalStrain, eqPlasticStrain, plasticWork) +subroutine utilities_calcPlasticity(yieldStress, plasticStrain, eqStress, eqTotalStrain, & + eqPlasticStrain, plasticWork, rotation_BC) use crystallite, only: & - crystallite_Fp, & crystallite_Fe, & crystallite_P, & crystallite_subF @@ -1065,6 +1065,7 @@ subroutine utilities_calcPlasticity(yieldStress, plasticStrain, eqStress, eqTota math_mul33x33, & math_trace33, & math_transpose33, & + math_rotate_forward33, & math_identity2nd, & math_crossproduct, & math_eigenvectorBasisSym, & @@ -1074,95 +1075,90 @@ subroutine utilities_calcPlasticity(yieldStress, plasticStrain, eqStress, eqTota implicit none - real(pReal), intent(out) :: eqStress, eqTotalStrain, eqPlasticStrain, plasticWork + real(pReal), intent(inout) :: eqStress, eqPlasticStrain, plasticWork + real(pReal), intent(out) :: eqTotalStrain real(pReal), dimension(3,3),intent(out) :: yieldStress, plasticStrain - real(pReal), dimension(3,3) :: cauchy, cauchy_av, P, Vp_av, V_total_av !< average + real(pReal), intent(in), dimension(3,3) :: rotation_BC !< rotation of load frame + real(pReal), dimension(3,3) :: cauchy, P_av, F_av, Ve_av !< average real(pReal), dimension(3) :: Values, S real(pReal), dimension(3,3) :: Vectors, diag + real(pReal), dimension(3,3) :: & + Vp, F_temp, U, VT, R, V, V_total real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & - Be, Ve, Vp, F, F_temp, Fe, Fp, U, VT, R, V, V_total + Be, Ve, Fe real(pReal), dimension(15) :: WORK !< previous deformation gradient integer(pInt) :: INFO, i, j, k, l, ierr - real(pReal) :: stress, stress_av, strain_total, strain_total_av, strain_plastic, strain_plastic_av, wgtm - real(pReal) :: eqStressOld, eqPlasticStrainOld + real(pReal) :: wgtm + real(pReal) :: eqStressOld, eqPlasticStrainOld, plasticWorkOld external :: dgesvd eqStressOld = eqStress eqPlasticStrainOld = eqPlasticStrain + plasticWorkOld = plasticWork wgtm = 1.0/real(mesh_NcpElems*mesh_maxNips*homogenization_maxNgrains,pReal) - Vp_av = 0.0_pReal - V_total_av = 0.0_pReal - strain_plastic_av = 0.0_pReal - strain_total_av = 0.0_pReal - cauchy_av = 0.0_pReal diag = 0.0_pReal + + P_av = sum(sum(sum(crystallite_P,dim=5),dim=4),dim=3) * wgtm + call MPI_Allreduce(MPI_IN_PLACE,P_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) + P_av = math_rotate_forward33(P_av,rotation_BC) + + F_av = sum(sum(sum(crystallite_subF,dim=5),dim=4),dim=3) * wgtm + call MPI_Allreduce(MPI_IN_PLACE,F_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) + + cauchy = 1.0/math_det33(F_av)*math_mul33x33(P_av,transpose(F_av)) + yieldStress = cauchy + eqStress = Mises(cauchy, 'stress') + + F_temp = F_av + call dgesvd ('A', 'A', 3, 3, F_temp, 3, S, U, 3, VT, 3, WORK, 15, INFO) ! singular value decomposition + + R = math_mul33x33(U, VT) ! rotation of polar decomposition + V = math_mul33x33(F_av,math_inv33(R)) + + call math_eigenValuesVectorsSym33(V,Values,Vectors) + do l = 1_pInt, 3_pInt + if (Values(l) < 0.0_pReal) then + Values(l) = -Values(l) + Vectors(1:3, l) = -Vectors(1:3, l) + endif + Values(l) = log(Values(l)) + diag(l,l) = Values(l) + enddo + if (dot_product(Vectors(1:3,1),Vectors(1:3,2)) /= 0) then + Vectors(1:3,2) = math_crossproduct(Vectors(1:3,3), Vectors(1:3,1)) + Vectors(1:3,2) = Vectors(1:3,2)/sqrt(dot_product(Vectors(1:3,2),Vectors(1:3,2))) + endif + if (dot_product(Vectors(1:3,2),Vectors(1:3,3)) /= 0) then + Vectors(1:3,3) = math_crossproduct(Vectors(1:3,1), Vectors(1:3,2)) + Vectors(1:3,3) = Vectors(1:3,3)/sqrt(dot_product(Vectors(1:3,3),Vectors(1:3,3))) + endif + if (dot_product(Vectors(1:3,3),Vectors(1:3,1)) /= 0) then + Vectors(1:3,1) = math_crossproduct(Vectors(1:3,2), Vectors(1:3,3)) + Vectors(1:3,1) = Vectors(1:3,1)/sqrt(dot_product(Vectors(1:3,1),Vectors(1:3,1))) + endif + + V_total = REAL(math_mul33x33(Vectors, math_mul33x33(diag, transpose(Vectors)))) + eqTotalStrain = Mises(V_total, 'strain') + do k = 1_pInt, mesh_NcpElems; do j = 1_pInt, mesh_maxNips; do i = 1_pInt,homogenization_maxNgrains - F(1:3,1:3,i,j,k) = crystallite_subF(1:3,1:3,i,j,k) - Fp(1:3,1:3,i,j,k) = crystallite_Fp(1:3,1:3,i,j,k) Fe(1:3,1:3,i,j,k) = crystallite_Fe(1:3,1:3,i,j,k) - - P = crystallite_P(1:3,1:3,i,j,k) - cauchy = 1.0/math_det33(F(1:3,1:3,i,j,k))*math_mul33x33(P,transpose(F(1:3,1:3,i,j,k))) - cauchy_av = cauchy_av + cauchy - stress = Mises(cauchy, 'stress') - stress_av = stress_av + stress - - F_temp(1:3,1:3,i,j,k) = F(1:3,1:3,i,j,k) - call dgesvd ('A', 'A', 3, 3, F_temp(1:3,1:3,i,j,k), 3, & - S, U(1:3,1:3,i,j,k), 3, VT(1:3,1:3,i,j,k), 3, WORK, 15, INFO) ! singular value decomposition - R(1:3,1:3,i,j,k) = math_mul33x33(U(1:3,1:3,i,j,k), VT(1:3,1:3,i,j,k)) ! rotation of polar decomposition - V(1:3,1:3,i,j,k) = math_mul33x33(F(1:3,1:3,i,j,k),math_inv33(R(1:3,1:3,i,j,k))) - - call math_eigenValuesVectorsSym33(V(1:3,1:3,i,j,k),Values,Vectors) - do l = 1_pInt, 3_pInt - if (Values(l) < 0.0_pReal) then - Values(l) = -Values(l) - Vectors(1:3, l) = -Vectors(1:3, l) - endif - Values(l) = log(Values(l)) - diag(l,l) = Values(l) - enddo - if (dot_product(Vectors(1:3,1),Vectors(1:3,2)) /= 0) then - Vectors(1:3,2) = math_crossproduct(Vectors(1:3,3), Vectors(1:3,1)) - Vectors(1:3,2) = Vectors(1:3,2)/sqrt(dot_product(Vectors(1:3,2),Vectors(1:3,2))) - endif - if (dot_product(Vectors(1:3,2),Vectors(1:3,3)) /= 0) then - Vectors(1:3,3) = math_crossproduct(Vectors(1:3,1), Vectors(1:3,2)) - Vectors(1:3,3) = Vectors(1:3,3)/sqrt(dot_product(Vectors(1:3,3),Vectors(1:3,3))) - endif - if (dot_product(Vectors(1:3,3),Vectors(1:3,1)) /= 0) then - Vectors(1:3,1) = math_crossproduct(Vectors(1:3,2), Vectors(1:3,3)) - Vectors(1:3,1) = Vectors(1:3,1)/sqrt(dot_product(Vectors(1:3,1),Vectors(1:3,1))) - endif - - V_total(1:3,1:3,i,j,k) = REAL(math_mul33x33(Vectors, math_mul33x33(diag, transpose(Vectors)))) - strain_total = Mises(V_total(1:3,1:3,i,j,k), 'strain') - strain_total_av = strain_total_av + strain_total - Be(1:3,1:3,i,j,k) = math_mul33x33(Fe(1:3,1:3,i,j,k),math_transpose33(Fe(1:3,1:3,i,j,k))) ! plastic part of left Cauchy–Green deformation tensor Ve(1:3,1:3,i,j,k) = math_eigenvectorBasisSym33_log(Be(1:3,1:3,i,j,k)) - Vp(1:3,1:3,i,j,k) = V_total(1:3,1:3,i,j,k) - Ve(1:3,1:3,i,j,k) - strain_plastic = Mises(Vp(1:3,1:3,i,j,k), 'strain') - strain_plastic_av = strain_plastic_av + strain_plastic - Vp_av = Vp_av + Vp(1:3,1:3,i,j,k) - enddo; enddo; enddo - yieldStress = cauchy_av * wgtm - call MPI_Allreduce(MPI_IN_PLACE,yieldStress,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) - plasticStrain = Vp_av * wgtm - call MPI_Allreduce(MPI_IN_PLACE,plasticStrain,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) + Ve_av = sum(sum(sum(Ve,dim=5),dim=4),dim=3) * wgtm + call MPI_Allreduce(MPI_IN_PLACE,Ve_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) - eqStress = stress_av * wgtm - call MPI_Allreduce(MPI_IN_PLACE,eqStress,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) - eqTotalStrain = strain_total_av * wgtm - call MPI_Allreduce(MPI_IN_PLACE,eqTotalStrain,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) - eqPlasticStrain = strain_plastic_av * wgtm - call MPI_Allreduce(MPI_IN_PLACE,eqPlasticStrain,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) - plasticWork = plasticWork + 0.5*(eqStressOld + eqStress) * (eqPlasticStrain - eqPlasticStrainOld) + Vp = V_total - Ve_av -end subroutine utilities_calcPlasticity + eqPlasticStrain = Mises(Vp, 'strain') + + plasticStrain = Vp + + plasticWork = plasticWorkOld + 0.5*(eqStressOld + eqStress) * (eqPlasticStrain - eqPlasticStrainOld) + +end subroutine utilities_calcPlasticity !-------------------------------------------------------------------------------------------------- !> @brief vonMises equivalent values for symmetric part of requested strains and/or stresses. From d81870dc57e5e04cfc8251f9544af81d04341e4f Mon Sep 17 00:00:00 2001 From: Fengbo Han Date: Wed, 22 Nov 2017 09:02:35 +0100 Subject: [PATCH 14/27] output the stress-strain curve to file if yield stop criterion is used --- src/DAMASK_spectral.f90 | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/src/DAMASK_spectral.f90 b/src/DAMASK_spectral.f90 index e51f7636a..3f6709d97 100755 --- a/src/DAMASK_spectral.f90 +++ b/src/DAMASK_spectral.f90 @@ -166,6 +166,7 @@ program DAMASK_spectral yieldStopValue real(pReal), dimension(3,3) :: yieldStress,yieldStressOld,yieldStressNew, plasticStrainOld, plasticStrainNew, plasticStrainRate integer(pInt) :: yieldResUnit = 0_pInt + integer(pInt) :: stressstrainUnit = 0_pInt character(len=13) :: stopFlag logical :: yieldStop, yieldStopSatisfied ! logical :: & @@ -710,6 +711,21 @@ program DAMASK_spectral call utilities_calcPlasticity(yieldStressNew, plasticStrainNew, eqStressNew, eqTotalStrainNew, & eqPlasticStrainNew, plasticWorkNew, loadCases(currentLoadCase)%rotation) + if (worldrank == 0) then ! output the stress-strain curve to file if yield stop criterion is used + if ((currentLoadCase == 1_pInt) .and. (inc == 1_pInt)) then + open(newunit=stressstrainUnit,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//& + '.stressstrain',form='FORMATTED',status='REPLACE') + write(stressstrainUnit,*) 0.0_pReal, 0.0_pReal + write(stressstrainUnit,*) eqTotalStrainNew, eqStressNew + close(stressstrainUnit) + else + open(newunit=stressstrainUnit,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//& + '.stressstrain',form='FORMATTED', position='APPEND', status='OLD') + write(stressstrainUnit,*) eqTotalStrainNew, eqStressNew + close(stressstrainUnit) + endif + endif + if(stopFlag == 'totalStrain') then if(eqTotalStrainNew > yieldStopValue) then yieldStress = yieldStressOld * (eqTotalStrainNew - yieldStopValue)/(eqTotalStrainNew - eqTotalStrainOld) & ! linear interpolation of stress values From 710970d728443310dde5da1c0c914cfc1da4c9e1 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Mon, 29 Jan 2018 21:27:05 -0500 Subject: [PATCH 15/27] updated option --label for addCurl/Div/Grad --- PRIVATE | 2 +- lib/damask/test/test.py | 2 +- processing/post/addCurl.py | 2 +- processing/post/addDivergence.py | 145 ++++++++++++++++--------------- processing/post/addGradient.py | 2 +- 5 files changed, 81 insertions(+), 72 deletions(-) diff --git a/PRIVATE b/PRIVATE index 1c1e80084..bdbf2da71 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 1c1e8008489d81773a13a247604144a8d7ee3723 +Subproject commit bdbf2da71cd9e0825d17f673ec2fbabc2c8027f8 diff --git a/lib/damask/test/test.py b/lib/damask/test/test.py index cbf330044..341a85c97 100644 --- a/lib/damask/test/test.py +++ b/lib/damask/test/test.py @@ -109,7 +109,7 @@ class Test(): except Exception as e: logging.critical('exception during variant execution: "{}"'.format(str(e))) - return variant+1 # return culprit + return variant+1 # return culprit return 0 def feasible(self): diff --git a/processing/post/addCurl.py b/processing/post/addCurl.py index 506b14282..01499b672 100755 --- a/processing/post/addCurl.py +++ b/processing/post/addCurl.py @@ -66,7 +66,7 @@ parser.add_option('-p','--pos','--periodiccellcenter', dest = 'pos', type = 'string', metavar = 'string', help = 'label of coordinates [%default]') -parser.add_option('-d','--data', +parser.add_option('-l','--label', dest = 'data', action = 'extend', metavar = '', help = 'label(s) of field values') diff --git a/processing/post/addDivergence.py b/processing/post/addDivergence.py index 232b5bc21..d7905630e 100755 --- a/processing/post/addDivergence.py +++ b/processing/post/addDivergence.py @@ -9,36 +9,43 @@ import damask scriptName = os.path.splitext(os.path.basename(__file__))[0] scriptID = ' '.join([scriptName,damask.version]) +def merge_dicts(*dict_args): + """Given any number of dicts, shallow copy and merge into a new dict, with precedence going to key value pairs in latter dicts.""" + result = {} + for dictionary in dict_args: + result.update(dictionary) + return result + def divFFT(geomdim,field): - shapeFFT = np.array(np.shape(field))[0:3] - grid = np.array(np.shape(field)[2::-1]) - N = grid.prod() # field size - n = np.array(np.shape(field)[3:]).prod() # data size + """Calculate divergence of a vector or tensor field by transforming into Fourier space.""" + shapeFFT = np.array(np.shape(field))[0:3] + grid = np.array(np.shape(field)[2::-1]) + N = grid.prod() # field size + n = np.array(np.shape(field)[3:]).prod() # data size - if n == 3: dataType = 'vector' - elif n == 9: dataType = 'tensor' + field_fourier = np.fft.rfftn(field,axes=(0,1,2),s=shapeFFT) + div_fourier = np.empty(field_fourier.shape[0:len(np.shape(field))-1],'c16') - field_fourier = np.fft.rfftn(field,axes=(0,1,2),s=shapeFFT) - div_fourier = np.empty(field_fourier.shape[0:len(np.shape(field))-1],'c16') + # differentiation in Fourier space + TWOPIIMG = 2.0j*math.pi + einsums = { + 3:'ijkl,ijkl->ijk', # vector, 3 -> 1 + 9:'ijkm,ijklm->ijkl', # tensor, 3x3 -> 3 + } + k_sk = np.where(np.arange(grid[2])>grid[2]//2,np.arange(grid[2])-grid[2],np.arange(grid[2]))/geomdim[0] + if grid[2]%2 == 0: k_sk[grid[2]//2] = 0 # Nyquist freq=0 for even grid (Johnson, MIT, 2011) -# differentiation in Fourier space - TWOPIIMG = 2.0j*math.pi - k_sk = np.where(np.arange(grid[2])>grid[2]//2,np.arange(grid[2])-grid[2],np.arange(grid[2]))/geomdim[0] - if grid[2]%2 == 0: k_sk[grid[2]//2] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) - - k_sj = np.where(np.arange(grid[1])>grid[1]//2,np.arange(grid[1])-grid[1],np.arange(grid[1]))/geomdim[1] - if grid[1]%2 == 0: k_sj[grid[1]//2] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) + k_sj = np.where(np.arange(grid[1])>grid[1]//2,np.arange(grid[1])-grid[1],np.arange(grid[1]))/geomdim[1] + if grid[1]%2 == 0: k_sj[grid[1]//2] = 0 # Nyquist freq=0 for even grid (Johnson, MIT, 2011) - k_si = np.arange(grid[0]//2+1)/geomdim[2] - - kk, kj, ki = np.meshgrid(k_sk,k_sj,k_si,indexing = 'ij') - k_s = np.concatenate((ki[:,:,:,None],kj[:,:,:,None],kk[:,:,:,None]),axis = 3).astype('c16') - if dataType == 'tensor': # tensor, 3x3 -> 3 - div_fourier = np.einsum('ijklm,ijkm->ijkl',field_fourier,k_s)*TWOPIIMG - elif dataType == 'vector': # vector, 3 -> 1 - div_fourier = np.einsum('ijkl,ijkl->ijk',field_fourier,k_s)*TWOPIIMG - - return np.fft.irfftn(div_fourier,axes=(0,1,2),s=shapeFFT).reshape([N,n/3]) + k_si = np.arange(grid[0]//2+1)/geomdim[2] + + kk, kj, ki = np.meshgrid(k_sk,k_sj,k_si,indexing = 'ij') + k_s = np.concatenate((ki[:,:,:,None],kj[:,:,:,None],kk[:,:,:,None]),axis = 3).astype('c16') + + div_fourier = np.einsum(einsums[n],k_s,field_fourier)*TWOPIIMG + + return np.fft.irfftn(div_fourier,axes=(0,1,2),s=shapeFFT).reshape([N,n/3]) # -------------------------------------------------------------------- @@ -46,32 +53,38 @@ def divFFT(geomdim,field): # -------------------------------------------------------------------- parser = OptionParser(option_class=damask.extendableOption, usage='%prog option(s) [ASCIItable(s)]', description = """ -Add column(s) containing divergence of requested column(s). -Operates on periodic ordered three-dimensional data sets. -Deals with both vector- and tensor-valued fields. - +Add column(s) containing curl of requested column(s). +Operates on periodic ordered three-dimensional data sets +of vector and tensor fields. """, version = scriptID) parser.add_option('-p','--pos','--periodiccellcenter', dest = 'pos', type = 'string', metavar = 'string', help = 'label of coordinates [%default]') -parser.add_option('-v','--vector', - dest = 'vector', +parser.add_option('-l','--label', + dest = 'data', action = 'extend', metavar = '', - help = 'label(s) of vector field values') -parser.add_option('-t','--tensor', - dest = 'tensor', - action = 'extend', metavar = '', - help = 'label(s) of tensor field values') + help = 'label(s) of field values') parser.set_defaults(pos = 'pos', ) + (options,filenames) = parser.parse_args() -if options.vector is None and options.tensor is None: - parser.error('no data column specified.') +if options.data is None: parser.error('no data column specified.') + +# --- define possible data types ------------------------------------------------------------------- + +datatypes = { + 3: {'name': 'vector', + 'shape': [3], + }, + 9: {'name': 'tensor', + 'shape': [3,3], + }, + } # --- loop over input files ------------------------------------------------------------------------ @@ -82,30 +95,27 @@ for name in filenames: except: continue damask.util.report(scriptName,name) -# ------------------------------------------ read header ------------------------------------------ +# --- interpret header ---------------------------------------------------------------------------- table.head_read() -# ------------------------------------------ sanity checks ---------------------------------------- - - items = { - 'tensor': {'dim': 9, 'shape': [3,3], 'labels':options.tensor, 'active':[], 'column': []}, - 'vector': {'dim': 3, 'shape': [3], 'labels':options.vector, 'active':[], 'column': []}, - } - errors = [] remarks = [] - column = {} - - if table.label_dimension(options.pos) != 3: errors.append('coordinates {} are not a vector.'.format(options.pos)) - else: colCoord = table.label_index(options.pos) + errors = [] + active = [] - for type, data in items.iteritems(): - for what in (data['labels'] if data['labels'] is not None else []): - dim = table.label_dimension(what) - if dim != data['dim']: remarks.append('column {} is not a {}.'.format(what,type)) - else: - items[type]['active'].append(what) - items[type]['column'].append(table.label_index(what)) + coordDim = table.label_dimension(options.pos) + if coordDim != 3: + errors.append('coordinates "{}" must be three-dimensional.'.format(options.pos)) + else: coordCol = table.label_index(options.pos) + + for me in options.data: + dim = table.label_dimension(me) + if dim in datatypes: + active.append(merge_dicts({'label':me},datatypes[dim])) + remarks.append('differentiating {} "{}"...'.format(datatypes[dim]['name'],me)) + else: + remarks.append('skipping "{}" of dimension {}...'.format(me,dim) if dim != -1 else \ + '"{}" not found...'.format(me) ) if remarks != []: damask.util.croak(remarks) if errors != []: @@ -116,17 +126,17 @@ for name in filenames: # ------------------------------------------ assemble header -------------------------------------- table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:])) - for type, data in items.iteritems(): - for label in data['active']: - table.labels_append(['divFFT({})'.format(label) if type == 'vector' else - '{}_divFFT({})'.format(i+1,label) for i in range(data['dim']//3)]) # extend ASCII header with new labels + for data in active: + table.labels_append(['divFFT({})'.format(data['label']) if data['shape'] == [3] \ + else '{}_divFFT({})'.format(i+1,data['label']) + for i in range(np.prod(np.array(data['shape']))//3)]) # extend ASCII header with new labels table.head_write() # --------------- figure out size and grid --------------------------------------------------------- table.data_readArray() - coords = [np.unique(table.data[:,colCoord+i]) for i in range(3)] + coords = [np.unique(table.data[:,coordCol+i]) for i in range(3)] mincorner = np.array(map(min,coords)) maxcorner = np.array(map(max,coords)) grid = np.array(map(len,coords),'i') @@ -136,12 +146,11 @@ for name in filenames: # ------------------------------------------ process value field ----------------------------------- stack = [table.data] - for type, data in items.iteritems(): - for i,label in enumerate(data['active']): - # we need to reverse order here, because x is fastest,ie rightmost, but leftmost in our x,y,z notation - stack.append(divFFT(size[::-1], - table.data[:,data['column'][i]:data['column'][i]+data['dim']]. - reshape(grid[::-1].tolist()+data['shape']))) + for data in active: + # we need to reverse order here, because x is fastest,ie rightmost, but leftmost in our x,y,z notation + stack.append(divFFT(size[::-1], + table.data[:,table.label_indexrange(data['label'])]. + reshape(grid[::-1].tolist()+data['shape']))) # ------------------------------------------ output result ----------------------------------------- diff --git a/processing/post/addGradient.py b/processing/post/addGradient.py index 4b9ef5a22..5a73af3a3 100755 --- a/processing/post/addGradient.py +++ b/processing/post/addGradient.py @@ -63,7 +63,7 @@ parser.add_option('-p','--pos','--periodiccellcenter', dest = 'pos', type = 'string', metavar = 'string', help = 'label of coordinates [%default]') -parser.add_option('-d','--data', +parser.add_option('-l','--label', dest = 'data', action = 'extend', metavar = '', help = 'label(s) of field values') From dc77e996cc02a3ad219eb93b17bf872344bfc963 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Tue, 30 Jan 2018 09:18:59 -0500 Subject: [PATCH 16/27] added utility function coordGridAndSize --- lib/damask/util.py | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) diff --git a/lib/damask/util.py b/lib/damask/util.py index 54de9acf9..d35071b7f 100644 --- a/lib/damask/util.py +++ b/lib/damask/util.py @@ -100,6 +100,18 @@ def execute(cmd, if process.returncode != 0: raise RuntimeError('{} failed with returncode {}'.format(cmd,process.returncode)) return out,error + +def coordGridAndSize(coordinates): + """Determines grid count and overall physical size along each dimension of an ordered array of coordinates""" + dim = coordinates.shape[1] + coords = [np.unique(coordinates[:,i]) for i in range(dim)] + mincorner = np.array(map(min,coords)) + maxcorner = np.array(map(max,coords)) + grid = np.array(map(len,coords),'i') + size = grid/np.maximum(np.ones(dim,'d'), grid-1.0) * (maxcorner-mincorner) # size from edge to edge = dim * n/(n-1) + size = np.where(grid > 1, size, min(size[grid > 1]/grid[grid > 1])) # spacing for grid==1 equal to smallest among other ones + return grid,size + # ----------------------------- class extendableOption(Option): """ @@ -130,7 +142,7 @@ class backgroundMessage(threading.Thread): 'hexagon': ['⬢', '⬣'], 'square': ['▖', '▘', '▝', '▗'], 'triangle': ['ᐊ', 'ᐊ', 'ᐃ', 'ᐅ', 'ᐅ', 'ᐃ'], - 'amoeba': ['▖', '▏', '▘', '▔', '▝', '▕', '▗', '▂'], + 'amoeba': ['▖', '▏', '▘', '▔', '▝', '▕', '▗', '▁'], 'beat': ['▁', '▂', '▃', '▅', '▆', '▇', '▇', '▆', '▅', '▃', '▂'], 'prison': ['ᚋ', 'ᚌ', 'ᚍ', 'ᚏ', 'ᚎ', 'ᚍ', 'ᚌ', 'ᚋ'], 'breath': ['ᚐ', 'ᚑ', 'ᚒ', 'ᚓ', 'ᚔ', 'ᚓ', 'ᚒ', 'ᚑ', 'ᚐ'], From 355d576b4dbe292f3c94f73c7b2900ca73cd626a Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Tue, 30 Jan 2018 09:20:47 -0500 Subject: [PATCH 17/27] shortened code with utility function coordGridAndSize --- processing/post/addCurl.py | 7 +------ processing/post/addDivergence.py | 7 +------ processing/post/addGaussian.py | 14 ++++---------- processing/post/addGradient.py | 7 +------ 4 files changed, 7 insertions(+), 28 deletions(-) diff --git a/processing/post/addCurl.py b/processing/post/addCurl.py index 01499b672..c197c5664 100755 --- a/processing/post/addCurl.py +++ b/processing/post/addCurl.py @@ -139,12 +139,7 @@ for name in filenames: table.data_readArray() - coords = [np.unique(table.data[:,coordCol+i]) for i in range(3)] - mincorner = np.array(map(min,coords)) - maxcorner = np.array(map(max,coords)) - grid = np.array(map(len,coords),'i') - size = grid/np.maximum(np.ones(3,'d'), grid-1.0) * (maxcorner-mincorner) # size from edge to edge = dim * n/(n-1) - size = np.where(grid > 1, size, min(size[grid > 1]/grid[grid > 1])) # spacing for grid==1 equal to smallest among other ones + grid,size = damask.util.gridAndSize(table.data[:,table.label_indexrange(options.pos)]) # ------------------------------------------ process value field ----------------------------------- diff --git a/processing/post/addDivergence.py b/processing/post/addDivergence.py index d7905630e..98916f56c 100755 --- a/processing/post/addDivergence.py +++ b/processing/post/addDivergence.py @@ -136,12 +136,7 @@ for name in filenames: table.data_readArray() - coords = [np.unique(table.data[:,coordCol+i]) for i in range(3)] - mincorner = np.array(map(min,coords)) - maxcorner = np.array(map(max,coords)) - grid = np.array(map(len,coords),'i') - size = grid/np.maximum(np.ones(3,'d'), grid-1.0) * (maxcorner-mincorner) # size from edge to edge = dim * n/(n-1) - size = np.where(grid > 1, size, min(size[grid > 1]/grid[grid > 1])) # spacing for grid==1 equal to smallest among other ones + grid,size = damask.util.coordGridAndSize(table.data[:,table.label_indexrange(options.pos)]) # ------------------------------------------ process value field ----------------------------------- diff --git a/processing/post/addGaussian.py b/processing/post/addGaussian.py index bc5599d49..c198ef62f 100755 --- a/processing/post/addGaussian.py +++ b/processing/post/addGaussian.py @@ -18,7 +18,7 @@ scriptID = ' '.join([scriptName,damask.version]) parser = OptionParser(option_class=damask.extendableOption, usage='%prog option(s) [ASCIItable(s)]', description = """ Add column(s) containing Gaussian filtered values of requested column(s). Operates on periodic and non-periodic ordered three-dimensional data sets. -For Details see scipy.ndimage documentation. +For details see scipy.ndimage documentation. """, version = scriptID) @@ -43,15 +43,14 @@ parser.add_option('--sigma', parser.add_option('--periodic', dest = 'periodic', action = 'store_true', - help = 'assume periodic grain structure' - ) + help = 'assume periodic grain structure') parser.set_defaults(pos = 'pos', order = 0, sigma = 1, - periodic = False + periodic = False, ) (options,filenames) = parser.parse_args() @@ -110,12 +109,7 @@ for name in filenames: table.data_readArray() - coords = [np.unique(table.data[:,colCoord+i]) for i in range(3)] - mincorner = np.array(map(min,coords)) - maxcorner = np.array(map(max,coords)) - grid = np.array(map(len,coords),'i') - size = grid/np.maximum(np.ones(3,'d'), grid-1.0) * (maxcorner-mincorner) # size from edge to edge = dim * n/(n-1) - size = np.where(grid > 1, size, min(size[grid > 1]/grid[grid > 1])) # spacing for grid==1 equal to smallest among other ones + grid,size = damask.util.coordGridAndSize(table.data[:,table.label_indexrange(options.pos)]) # ------------------------------------------ process value field ----------------------------------- diff --git a/processing/post/addGradient.py b/processing/post/addGradient.py index 5a73af3a3..c788f5286 100755 --- a/processing/post/addGradient.py +++ b/processing/post/addGradient.py @@ -135,12 +135,7 @@ for name in filenames: table.data_readArray() - coords = [np.unique(table.data[:,coordCol+i]) for i in range(3)] - mincorner = np.array(map(min,coords)) - maxcorner = np.array(map(max,coords)) - grid = np.array(map(len,coords),'i') - size = grid/np.maximum(np.ones(3,'d'), grid-1.0) * (maxcorner-mincorner) # size from edge to edge = dim * n/(n-1) - size = np.where(grid > 1, size, min(size[grid > 1]/grid[grid > 1])) # spacing for grid==1 equal to smallest among other ones + grid,size = damask.util.coordGridAndSize(table.data[:,table.label_indexrange(options.pos)]) # ------------------------------------------ process value field ----------------------------------- From e304ce35dacb0e28dbaf63a1105a8ef6ce47e1f1 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Tue, 30 Jan 2018 12:58:43 -0500 Subject: [PATCH 18/27] forgot to rename function call to read "coordGridAndSize" --- processing/post/addCurl.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/processing/post/addCurl.py b/processing/post/addCurl.py index c197c5664..5ca851b22 100755 --- a/processing/post/addCurl.py +++ b/processing/post/addCurl.py @@ -139,7 +139,7 @@ for name in filenames: table.data_readArray() - grid,size = damask.util.gridAndSize(table.data[:,table.label_indexrange(options.pos)]) + grid,size = damask.util.coordGridAndSize(table.data[:,table.label_indexrange(options.pos)]) # ------------------------------------------ process value field ----------------------------------- From e4266c94e77daa268b5e5cc31e663b109704fba5 Mon Sep 17 00:00:00 2001 From: Test User Date: Wed, 31 Jan 2018 04:14:24 +0100 Subject: [PATCH 19/27] [skip ci] updated version information after successful test of v2.0.1-1029-ge304ce3 --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index 3990ae29a..3c048e8c8 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.1-1015-g532d669 +v2.0.1-1029-ge304ce3 From 60440b1a5474f2995b89903d2f7ae7f1d74f21d6 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Wed, 31 Jan 2018 12:09:14 -0500 Subject: [PATCH 20/27] polished output and streamlined code in test-class, added new consistency check in PRIVATE --- PRIVATE | 2 +- lib/damask/test/test.py | 53 ++++++++++++++++++++++------------------- 2 files changed, 30 insertions(+), 25 deletions(-) diff --git a/PRIVATE b/PRIVATE index bdbf2da71..5fc3188c8 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit bdbf2da71cd9e0825d17f673ec2fbabc2c8027f8 +Subproject commit 5fc3188c86ea1f4159db87529ac3e3169fb56e5d diff --git a/lib/damask/test/test.py b/lib/damask/test/test.py index 341a85c97..3869d45f2 100644 --- a/lib/damask/test/test.py +++ b/lib/damask/test/test.py @@ -322,8 +322,10 @@ class Test(): cur1Name = self.fileInCurrent(cur1) return self.compare_Array(cur0Name,cur1Name) - def compare_Table(self,headings0,file0,headings1,file1,normHeadings='',normType=None, - absoluteTolerance=False,perLine=False,skipLines=[]): + def compare_Table(self,headings0,file0, + headings1,file1, + normHeadings='',normType=None, + absoluteTolerance=False,perLine=False,skipLines=[]): import numpy as np logging.info('\n '.join(['comparing ASCII Tables',file0,file1])) @@ -337,7 +339,7 @@ class Test(): data = [[] for i in range(dataLength)] maxError = [0.0 for i in range(dataLength)] absTol = [absoluteTolerance for i in range(dataLength)] - column = [[1 for i in range(dataLength)] for j in range(2)] + column = [[1 for i in range(dataLength)] for j in range(2)] norm = [[] for i in range(dataLength)] normLength = [1 for i in range(dataLength)] @@ -368,11 +370,11 @@ class Test(): key1 = ('1_' if length[i]>1 else '') + headings1[i]['label'] normKey = ('1_' if normLength[i]>1 else '') + normHeadings[i]['label'] if key0 not in table0.labels(raw = True): - raise Exception('column {} not found in 1. table...\n'.format(key0)) + raise Exception('column "{}" not found in first table...\n'.format(key0)) elif key1 not in table1.labels(raw = True): - raise Exception('column {} not found in 2. table...\n'.format(key1)) + raise Exception('column "{}" not found in second table...\n'.format(key1)) elif normKey not in table0.labels(raw = True): - raise Exception('column {} not found in 1. table...\n'.format(normKey)) + raise Exception('column "{}" not found in first table...\n'.format(normKey)) else: column[0][i] = table0.label_index(key0) column[1][i] = table1.label_index(key1) @@ -400,9 +402,9 @@ class Test(): norm[i] = [1.0 for j in range(line0-len(skipLines))] absTol[i] = True if perLine: - logging.warning('At least one norm of {} in 1. table is 0.0, using absolute tolerance'.format(headings0[i]['label'])) + logging.warning('At least one norm of "{}" in first table is 0.0, using absolute tolerance'.format(headings0[i]['label'])) else: - logging.warning('Maximum norm of {} in 1. table is 0.0, using absolute tolerance'.format(headings0[i]['label'])) + logging.warning('Maximum norm of "{}" in first table is 0.0, using absolute tolerance'.format(headings0[i]['label'])) line1 = 0 while table1.data_read(): # read next data line of ASCII table @@ -414,7 +416,7 @@ class Test(): norm[i][line1-len(skipLines)]) line1 +=1 - if (line0 != line1): raise Exception('found {} lines in 1. table but {} in 2. table'.format(line0,line1)) + if (line0 != line1): raise Exception('found {} lines in first table but {} in second table'.format(line0,line1)) logging.info(' ********') for i in range(dataLength): @@ -561,25 +563,28 @@ class Test(): return allclose - def compare_TableRefCur(self,headingsRef,ref,headingsCur='',cur='',normHeadings='',normType=None,\ - absoluteTolerance=False,perLine=False,skipLines=[]): + def compare_TableRefCur(self,headingsRef,ref,headingsCur='',cur='', + normHeadings='',normType=None, + absoluteTolerance=False,perLine=False,skipLines=[]): - if cur == '': cur = ref - if headingsCur == '': headingsCur = headingsRef - refName = self.fileInReference(ref) - curName = self.fileInCurrent(cur) - return self.compare_Table(headingsRef,refName,headingsCur,curName,normHeadings,normType, - absoluteTolerance,perLine,skipLines) + return self.compare_Table(headingsRef, + self.fileInReference(ref), + headingsRef if headingsCur == '' else headingsCur, + self.fileInReference(ref if cur == '' else cur), + normHeadings,normType, + absoluteTolerance,perLine,skipLines) - def compare_TableCurCur(self,headingsCur0,Cur0,Cur1,headingsCur1='',normHeadings='',normType=None,\ - absoluteTolerance=False,perLine=False,skipLines=[]): + def compare_TableCurCur(self,headingsCur0,Cur0,Cur1, + headingsCur1='', + normHeadings='',normType=None, + absoluteTolerance=False,perLine=False,skipLines=[]): - if headingsCur1 == '': headingsCur1 = headingsCur0 - cur0Name = self.fileInCurrent(Cur0) - cur1Name = self.fileInCurrent(Cur1) - return self.compare_Table(headingsCur0,cur0Name,headingsCur1,cur1Name,normHeadings,normType, - absoluteTolerance,perLine,skipLines) + return self.compare_Table(headingsCur0, + self.fileInCurrent(Cur0), + headingsCur0 if headingsCur1 == '' else headingsCur1, + self.fileInCurrent(Cur1), + normHeadings,normType,absoluteTolerance,perLine,skipLines) def report_Success(self,culprit): From 85147605adcdf78fd1850670caaf4305a48a8c4d Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Wed, 31 Jan 2018 12:35:50 -0500 Subject: [PATCH 21/27] bugfix: use fileInCurrent instead of fileInRef --- lib/damask/test/test.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/damask/test/test.py b/lib/damask/test/test.py index 3869d45f2..0e1a0284c 100644 --- a/lib/damask/test/test.py +++ b/lib/damask/test/test.py @@ -570,7 +570,7 @@ class Test(): return self.compare_Table(headingsRef, self.fileInReference(ref), headingsRef if headingsCur == '' else headingsCur, - self.fileInReference(ref if cur == '' else cur), + self.fileInCurrent(ref if cur == '' else cur), normHeadings,normType, absoluteTolerance,perLine,skipLines) From 8776b525619c3f9c2ca22cb5a61034179482408b Mon Sep 17 00:00:00 2001 From: Test User Date: Thu, 1 Feb 2018 07:34:53 +0100 Subject: [PATCH 22/27] [skip ci] updated version information after successful test of v2.0.1-1033-g8514760 --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index 3c048e8c8..8f9e4d0c0 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.1-1029-ge304ce3 +v2.0.1-1033-g8514760 From d80a25573680a25307f374058a13e485fa8e8971 Mon Sep 17 00:00:00 2001 From: Franz Roters Date: Fri, 2 Feb 2018 15:06:13 +0100 Subject: [PATCH 23/27] new Marc2017 file format finally working! --- src/IO.f90 | 6 ++--- src/mesh.f90 | 64 ++++++++++++++++++++++++++++++++++++---------------- 2 files changed, 47 insertions(+), 23 deletions(-) diff --git a/src/IO.f90 b/src/IO.f90 index a6c3a7da8..d877379c7 100644 --- a/src/IO.f90 +++ b/src/IO.f90 @@ -1268,7 +1268,7 @@ integer(pInt) function IO_countNumericalDataLines(fileUnit) line = IO_read(fileUnit) chunkPos = IO_stringPos(line) tmp = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) - if (scan(tmp,"abcdefghijklmnopqrstuvwxyz")/=0) then ! found keyword + if (verify(trim(tmp) ,"0123456789")/=0) then ! found keyword line = IO_read(fileUnit, .true.) ! reset IO_read exit else @@ -1839,12 +1839,12 @@ integer(pInt) function IO_verifyIntValue (string,validChars,myName) invalidWhere = verify(string,validChars) if (invalidWhere == 0_pInt) then read(UNIT=string,iostat=readStatus,FMT=*) IO_verifyIntValue ! no offending chars found - if (readStatus /= 0_pInt) & ! error during string to float conversion + if (readStatus /= 0_pInt) & ! error during string to integer conversion call IO_warning(203_pInt,ext_msg=myName//'"'//string//'"') else call IO_warning(202_pInt,ext_msg=myName//'"'//string//'"') ! complain about offending characters read(UNIT=string(1_pInt:invalidWhere-1_pInt),iostat=readStatus,FMT=*) IO_verifyIntValue ! interpret remaining string - if (readStatus /= 0_pInt) & ! error during string to float conversion + if (readStatus /= 0_pInt) & ! error during string to integer conversion call IO_warning(203_pInt,ext_msg=myName//'"'//string(1_pInt:invalidWhere-1_pInt)//'"') endif diff --git a/src/mesh.f90 b/src/mesh.f90 index 827875d2f..d7d0f8c06 100644 --- a/src/mesh.f90 +++ b/src/mesh.f90 @@ -1642,11 +1642,12 @@ subroutine mesh_marc_get_matNumber(fileUnit) if (len(trim(line))/=0_pInt) then chunkPos = IO_stringPos(line) data_blocks = IO_intValue(line,chunkPos,1_pInt) - endif + endif + allocate(Marc_matNumber(data_blocks)) do i=1_pInt,data_blocks ! read all data blocks read (fileUnit,610,END=620) line chunkPos = IO_stringPos(line) - Marc_matNumber = (/Marc_matNumber, IO_intValue(line,chunkPos,1_pInt)/) + Marc_matNumber(i) = IO_intValue(line,chunkPos,1_pInt) do j=1_pint,2_pInt + hypoelasticTableStyle ! read 2 or 3 remaining lines of data block read (fileUnit,610,END=620) line enddo @@ -1809,12 +1810,11 @@ subroutine mesh_marc_count_cpElements(fileUnit) do i=1_pInt,3_pInt+hypoelasticTableStyle ! Skip 3 or 4 lines read (fileUnit,610,END=620) line enddo - mesh_NcpElems = mesh_NcpElems + IO_countContinuousIntValues(fileUnit) ! why not simply mesh_NcpElems = IO_countContinuousIntValues(fileUnit)? keyword hypoelastic might appear several times so the next line probably should not be there + mesh_NcpElems = mesh_NcpElems + IO_countContinuousIntValues(fileUnit) ! why not simply mesh_NcpElems = IO_countContinuousIntValues(fileUnit)? not fully correct as hypoelastic can have multiple data fields, needs update exit endif enddo - else ! Marc2017 and later - call IO_error(error_ID=701_pInt) + else ! Marc2017 and later do read (fileUnit,610,END=620) line chunkPos = IO_stringPos(line) @@ -1839,6 +1839,7 @@ subroutine mesh_marc_map_elements(fileUnit) use math, only: math_qsort use IO, only: IO_lc, & + IO_intValue, & IO_stringValue, & IO_stringPos, & IO_continuousIntValues @@ -1847,7 +1848,8 @@ subroutine mesh_marc_map_elements(fileUnit) integer(pInt), intent(in) :: fileUnit integer(pInt), allocatable, dimension(:) :: chunkPos - character(len=300) :: line + character(len=300) :: line, & + tmp integer(pInt), dimension (1_pInt+mesh_NcpElems) :: contInts integer(pInt) :: i,cpElem = 0_pInt @@ -1856,25 +1858,47 @@ subroutine mesh_marc_map_elements(fileUnit) 610 FORMAT(A300) + contInts = 0_pInt rewind(fileUnit) do read (fileUnit,610,END=660) line chunkPos = IO_stringPos(line) - if( IO_lc(IO_stringValue(line,chunkPos,1_pInt)) == 'hypoelastic' ) then - do i=1_pInt,3_pInt+hypoelasticTableStyle ! skip three (or four if new table style!) lines - read (fileUnit,610,END=660) line - enddo - contInts = IO_continuousIntValues(fileUnit,mesh_NcpElems,mesh_nameElemSet,& + if (MarcVersion < 13) then ! Marc 2016 or earlier + if( IO_lc(IO_stringValue(line,chunkPos,1_pInt)) == 'hypoelastic' ) then + do i=1_pInt,3_pInt+hypoelasticTableStyle ! skip three (or four if new table style!) lines + read (fileUnit,610,END=660) line + enddo + contInts = IO_continuousIntValues(fileUnit,mesh_NcpElems,mesh_nameElemSet,& mesh_mapElemSet,mesh_NelemSets) - do i = 1_pInt,contInts(1) - cpElem = cpElem+1_pInt - mesh_mapFEtoCPelem(1,cpElem) = contInts(1_pInt+i) - mesh_mapFEtoCPelem(2,cpElem) = cpElem - enddo - endif - enddo - -660 call math_qsort(mesh_mapFEtoCPelem,1_pInt,int(size(mesh_mapFEtoCPelem,2_pInt),pInt)) ! should be mesh_NcpElems + exit + endif + else ! Marc2017 and later + if ( IO_lc(IO_stringValue(line,chunkPos,1_pInt)) == 'connectivity') then + read (fileUnit,610,END=660) line + chunkPos = IO_stringPos(line) + if(any(Marc_matNumber==IO_intValue(line,chunkPos,6_pInt))) then + do + read (fileUnit,610,END=660) line + chunkPos = IO_stringPos(line) + tmp = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) + if (verify(trim(tmp),"0123456789")/=0) then ! found keyword + exit + else + contInts(1) = contInts(1) + 1_pInt + read (tmp,*) contInts(contInts(1)+1) + endif + enddo + endif + endif + endif + enddo +660 do i = 1_pInt,contInts(1) + cpElem = cpElem+1_pInt + mesh_mapFEtoCPelem(1,cpElem) = contInts(1_pInt+i) + mesh_mapFEtoCPelem(2,cpElem) = cpElem + enddo + +call math_qsort(mesh_mapFEtoCPelem,1_pInt,int(size(mesh_mapFEtoCPelem,2_pInt),pInt)) ! should be mesh_NcpElems end subroutine mesh_marc_map_elements From 911fb84e81a1154e77df472f014bc49516a917d9 Mon Sep 17 00:00:00 2001 From: Test User Date: Sat, 3 Feb 2018 03:57:56 +0100 Subject: [PATCH 24/27] [skip ci] updated version information after successful test of v2.0.1-1035-gd80a255 --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index 8f9e4d0c0..35191bcab 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.1-1033-g8514760 +v2.0.1-1035-gd80a255 From 190a2baf9f928141891618c53d6945b5a403ce29 Mon Sep 17 00:00:00 2001 From: Fengbo Han Date: Wed, 7 Feb 2018 11:35:16 +0100 Subject: [PATCH 25/27] when using yield stop criteria, if rotation of the load frame is specified, the output results in .yield and .stressstrain files are also rotated --- src/spectral_utilities.f90 | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/spectral_utilities.f90 b/src/spectral_utilities.f90 index 742068c9a..7417560a5 100755 --- a/src/spectral_utilities.f90 +++ b/src/spectral_utilities.f90 @@ -1078,15 +1078,15 @@ subroutine utilities_calcPlasticity(yieldStress, plasticStrain, eqStress, eqTota real(pReal), intent(inout) :: eqStress, eqPlasticStrain, plasticWork real(pReal), intent(out) :: eqTotalStrain real(pReal), dimension(3,3),intent(out) :: yieldStress, plasticStrain - real(pReal), intent(in), dimension(3,3) :: rotation_BC !< rotation of load frame - real(pReal), dimension(3,3) :: cauchy, P_av, F_av, Ve_av !< average + real(pReal), intent(in), dimension(3,3) :: rotation_BC !< rotation of load frame + real(pReal), dimension(3,3) :: cauchy, P_av, F_av, Ve_av !< average real(pReal), dimension(3) :: Values, S real(pReal), dimension(3,3) :: Vectors, diag real(pReal), dimension(3,3) :: & Vp, F_temp, U, VT, R, V, V_total real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & Be, Ve, Fe - real(pReal), dimension(15) :: WORK !< previous deformation gradient + real(pReal), dimension(15) :: WORK !< previous deformation gradient integer(pInt) :: INFO, i, j, k, l, ierr real(pReal) :: wgtm real(pReal) :: eqStressOld, eqPlasticStrainOld, plasticWorkOld @@ -1105,6 +1105,7 @@ subroutine utilities_calcPlasticity(yieldStress, plasticStrain, eqStress, eqTota F_av = sum(sum(sum(crystallite_subF,dim=5),dim=4),dim=3) * wgtm call MPI_Allreduce(MPI_IN_PLACE,F_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) + F_av = math_rotate_forward33(F_av,rotation_BC) cauchy = 1.0/math_det33(F_av)*math_mul33x33(P_av,transpose(F_av)) yieldStress = cauchy @@ -1143,7 +1144,8 @@ subroutine utilities_calcPlasticity(yieldStress, plasticStrain, eqStress, eqTota do k = 1_pInt, mesh_NcpElems; do j = 1_pInt, mesh_maxNips; do i = 1_pInt,homogenization_maxNgrains Fe(1:3,1:3,i,j,k) = crystallite_Fe(1:3,1:3,i,j,k) - Be(1:3,1:3,i,j,k) = math_mul33x33(Fe(1:3,1:3,i,j,k),math_transpose33(Fe(1:3,1:3,i,j,k))) ! plastic part of left Cauchy–Green deformation tensor + Fe(1:3,1:3,i,j,k) = math_rotate_forward33(Fe(1:3,1:3,i,j,k),rotation_BC) + Be(1:3,1:3,i,j,k) = math_mul33x33(Fe(1:3,1:3,i,j,k),math_transpose33(Fe(1:3,1:3,i,j,k))) ! elastic part of left Cauchy–Green deformation tensor Ve(1:3,1:3,i,j,k) = math_eigenvectorBasisSym33_log(Be(1:3,1:3,i,j,k)) enddo; enddo; enddo From b834b2a00da82d370f5b1927b4dad9684a4bec49 Mon Sep 17 00:00:00 2001 From: Fengbo Han Date: Wed, 7 Feb 2018 13:37:26 +0100 Subject: [PATCH 26/27] removed unnecessary lines --- src/DAMASK_spectral.f90 | 18 ++++++++---------- 1 file changed, 8 insertions(+), 10 deletions(-) diff --git a/src/DAMASK_spectral.f90 b/src/DAMASK_spectral.f90 index d70a818a6..beae759fc 100755 --- a/src/DAMASK_spectral.f90 +++ b/src/DAMASK_spectral.f90 @@ -159,20 +159,18 @@ program DAMASK_spectral PETScFinalize !-------------------------------------------------------------------------------------------------- ! variables related to stop criterion for yielding - real(pReal) :: plasticWorkOld, plasticWorkNew, & ! plastic work - eqTotalStrainOld, eqTotalStrainNew, & ! total equivalent strain - eqPlasticStrainOld, eqPlasticStrainNew, & ! total equivalent plastic strain - eqStressOld, eqStressNew , & ! equivalent stress + real(pReal) :: plasticWorkOld, plasticWorkNew, & ! plastic work + eqTotalStrainOld, eqTotalStrainNew, & ! total equivalent strain + eqPlasticStrainOld, eqPlasticStrainNew, & ! total equivalent plastic strain + eqStressOld, eqStressNew , & ! equivalent stress yieldStopValue - real(pReal), dimension(3,3) :: yieldStress,yieldStressOld,yieldStressNew, plasticStrainOld, plasticStrainNew, plasticStrainRate + real(pReal), dimension(3,3) :: yieldStress,yieldStressOld,yieldStressNew, & + plasticStrainOld, plasticStrainNew, plasticStrainRate integer(pInt) :: yieldResUnit = 0_pInt integer(pInt) :: stressstrainUnit = 0_pInt character(len=13) :: stopFlag logical :: yieldStop, yieldStopSatisfied - ! logical :: & - ! stop_totalStrain, & !< stop criterion - ! stop_plasticStrain, & - ! stop_plasticWork + !-------------------------------------------------------------------------------------------------- ! init DAMASK (all modules) call CPFEM_initAll(el = 1_pInt, ip = 1_pInt) @@ -321,7 +319,7 @@ program DAMASK_spectral enddo; enddo close(FILEUNIT) - if(yieldStop) then ! initialize variables related to yield stop + if(yieldStop) then ! initialize variables related to yield stop yieldStressNew = 0.0_pReal plasticStrainNew = 0.0_pReal eqStressNew = 0.0_pReal From 208c4affa452c0f3f033f24eceaa211fc684b1e7 Mon Sep 17 00:00:00 2001 From: Fengbo Han Date: Wed, 7 Feb 2018 17:11:43 +0100 Subject: [PATCH 27/27] using math_equivStrain33 and math_equivStress33 instead of Mises --- src/DAMASK_spectral.f90 | 6 +++--- src/spectral_utilities.f90 | 36 +++++++----------------------------- 2 files changed, 10 insertions(+), 32 deletions(-) diff --git a/src/DAMASK_spectral.f90 b/src/DAMASK_spectral.f90 index beae759fc..7b4cddd77 100755 --- a/src/DAMASK_spectral.f90 +++ b/src/DAMASK_spectral.f90 @@ -303,15 +303,15 @@ program DAMASK_spectral temp_valueVector(j) = IO_floatValue(line,chunkPos,i+j) enddo loadCases(currentLoadCase)%rotation = math_plain9to33(temp_valueVector) - case('totalStrain', 'totalstrain') + case('totalstrain') yieldStop = .True. stopFlag = 'totalStrain' yieldStopValue = IO_floatValue(line,chunkPos,i+1_pInt) - case('plasticStrain', 'plasticstrain') + case('plasticstrain') yieldStop = .True. stopFlag = 'plasticStrain' yieldStopValue = IO_floatValue(line,chunkPos,i+1_pInt) - case('plasticWork', 'plasticwork') + case('plasticwork') yieldStop = .True. stopFlag = 'plasticWork' yieldStopValue = IO_floatValue(line,chunkPos,i+1_pInt) diff --git a/src/spectral_utilities.f90 b/src/spectral_utilities.f90 index 1efaf8879..ea7f50a31 100755 --- a/src/spectral_utilities.f90 +++ b/src/spectral_utilities.f90 @@ -1068,6 +1068,8 @@ subroutine utilities_calcPlasticity(yieldStress, plasticStrain, eqStress, eqTota math_mul33x33, & math_trace33, & math_transpose33, & + math_equivStrain33, & + math_equivStress33, & math_rotate_forward33, & math_identity2nd, & math_crossproduct, & @@ -1099,7 +1101,7 @@ subroutine utilities_calcPlasticity(yieldStress, plasticStrain, eqStress, eqTota eqStressOld = eqStress eqPlasticStrainOld = eqPlasticStrain plasticWorkOld = plasticWork - wgtm = 1.0/real(mesh_NcpElems*mesh_maxNips*homogenization_maxNgrains,pReal) + wgtm = 1.0_pReal/real(mesh_NcpElems*mesh_maxNips*homogenization_maxNgrains,pReal) diag = 0.0_pReal P_av = sum(sum(sum(crystallite_P,dim=5),dim=4),dim=3) * wgtm @@ -1110,9 +1112,9 @@ subroutine utilities_calcPlasticity(yieldStress, plasticStrain, eqStress, eqTota call MPI_Allreduce(MPI_IN_PLACE,F_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) F_av = math_rotate_forward33(F_av,rotation_BC) - cauchy = 1.0/math_det33(F_av)*math_mul33x33(P_av,transpose(F_av)) + cauchy = 1.0_pReal/math_det33(F_av)*math_mul33x33(P_av,transpose(F_av)) yieldStress = cauchy - eqStress = Mises(cauchy, 'stress') + eqStress = math_equivStress33(cauchy) F_temp = F_av call dgesvd ('A', 'A', 3, 3, F_temp, 3, S, U, 3, VT, 3, WORK, 15, INFO) ! singular value decomposition @@ -1143,7 +1145,7 @@ subroutine utilities_calcPlasticity(yieldStress, plasticStrain, eqStress, eqTota endif V_total = REAL(math_mul33x33(Vectors, math_mul33x33(diag, transpose(Vectors)))) - eqTotalStrain = Mises(V_total, 'strain') + eqTotalStrain = math_equivStrain33(V_total) do k = 1_pInt, mesh_NcpElems; do j = 1_pInt, mesh_maxNips; do i = 1_pInt,homogenization_maxNgrains Fe(1:3,1:3,i,j,k) = crystallite_Fe(1:3,1:3,i,j,k) @@ -1157,7 +1159,7 @@ subroutine utilities_calcPlasticity(yieldStress, plasticStrain, eqStress, eqTota Vp = V_total - Ve_av - eqPlasticStrain = Mises(Vp, 'strain') + eqPlasticStrain = math_equivStrain33(Vp) plasticStrain = Vp @@ -1165,30 +1167,6 @@ subroutine utilities_calcPlasticity(yieldStress, plasticStrain, eqStress, eqTota end subroutine utilities_calcPlasticity -!-------------------------------------------------------------------------------------------------- -!> @brief vonMises equivalent values for symmetric part of requested strains and/or stresses. -!-------------------------------------------------------------------------------------------------- -function Mises(m, mType) - use math, only: & - math_trace33, & - math_identity2nd - - implicit none - real(pReal), intent(in), dimension(3,3) :: m - character(len=*), intent(in) :: mType - real(pReal) :: trace, Mises - real(pReal), dimension(3,3) :: eye, dev, symdev - - eye = math_identity2nd(3) - trace = math_trace33(m) - dev = m - trace/3.0*eye - symdev = 0.5*(dev+transpose(dev)) - if (mType == 'stress') then - Mises = sqrt(3.0/2.0*sum(symdev*transpose(symdev))) - else - Mises = sqrt(2.0/3.0*sum(symdev*transpose(symdev))) - endif -end function Mises !-------------------------------------------------------------------------------------------------- !> @brief calculates forward rate, either guessing or just add delta/timeinc