diff --git a/src/CPFEM.f90 b/src/CPFEM.f90 index fa774de66..a55b9e4ce 100644 --- a/src/CPFEM.f90 +++ b/src/CPFEM.f90 @@ -184,8 +184,8 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS temperature(material_homogenizationAt(elCP))%p(thermalMapping(material_homogenizationAt(elCP))%p(ip,elCP)) = & temperature_inp end select chosenThermal1 - materialpoint_F0(1:3,1:3,ip,elCP) = ffn - materialpoint_F(1:3,1:3,ip,elCP) = ffn1 + homogenization_F0(1:3,1:3,ip,elCP) = ffn + homogenization_F(1:3,1:3,ip,elCP) = ffn1 if (iand(mode, CPFEM_CALCRESULTS) /= 0_pInt) then @@ -212,17 +212,17 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS else terminalIllness ! translate from P to sigma - Kirchhoff = matmul(materialpoint_P(1:3,1:3,ip,elCP), transpose(materialpoint_F(1:3,1:3,ip,elCP))) - J_inverse = 1.0_pReal / math_det33(materialpoint_F(1:3,1:3,ip,elCP)) + Kirchhoff = matmul(homogenization_P(1:3,1:3,ip,elCP), transpose(homogenization_F(1:3,1:3,ip,elCP))) + J_inverse = 1.0_pReal / math_det33(homogenization_F(1:3,1:3,ip,elCP)) CPFEM_cs(1:6,ip,elCP) = math_sym33to6(J_inverse * Kirchhoff,weighted=.false.) ! translate from dP/dF to dCS/dE H = 0.0_pReal do i=1,3; do j=1,3; do k=1,3; do l=1,3; do m=1,3; do n=1,3 H(i,j,k,l) = H(i,j,k,l) & - + materialpoint_F(j,m,ip,elCP) * materialpoint_F(l,n,ip,elCP) & - * materialpoint_dPdF(i,m,k,n,ip,elCP) & - - math_delta(j,l) * materialpoint_F(i,m,ip,elCP) * materialpoint_P(k,m,ip,elCP) & + + homogenization_F(j,m,ip,elCP) * homogenization_F(l,n,ip,elCP) & + * homogenization_dPdF(i,m,k,n,ip,elCP) & + - math_delta(j,l) * homogenization_F(i,m,ip,elCP) * homogenization_P(k,m,ip,elCP) & + 0.5_pReal * ( Kirchhoff(j,l)*math_delta(i,k) + Kirchhoff(i,k)*math_delta(j,l) & + Kirchhoff(j,k)*math_delta(i,l) + Kirchhoff(i,l)*math_delta(j,k)) enddo; enddo; enddo; enddo; enddo; enddo diff --git a/src/grid/grid_mech_FEM.f90 b/src/grid/grid_mech_FEM.f90 index dfbd3f2f3..bf3d7752d 100644 --- a/src/grid/grid_mech_FEM.f90 +++ b/src/grid/grid_mech_FEM.f90 @@ -236,7 +236,7 @@ subroutine grid_mech_FEM_init F = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) endif restartRead - materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent + homogenization_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent call utilities_updateCoords(F) call utilities_constitutiveResponse(P_current,temp33_Real,C_volAvg,devNull, & ! stress field, stress avg, global average of stiffness and (min+max)/2 F, & ! target F @@ -364,7 +364,7 @@ subroutine grid_mech_FEM_forward(cutBack,guess,timeinc,timeinc_old,loadCaseTime, F_lastInc = F - materialpoint_F0 = reshape(F, [3,3,1,product(grid(1:2))*grid3]) + homogenization_F0 = reshape(F, [3,3,1,product(grid(1:2))*grid3]) endif !-------------------------------------------------------------------------------------------------- @@ -557,9 +557,9 @@ subroutine formResidual(da_local,x_local, & ii = i-xstart+1; jj = j-ystart+1; kk = k-zstart+1 ele = ele + 1 f_elem = matmul(transpose(BMat),transpose(P_current(1:3,1:3,ii,jj,kk)))*detJ + & - matmul(HGMat,x_elem)*(materialpoint_dPdF(1,1,1,1,1,ele) + & - materialpoint_dPdF(2,2,2,2,1,ele) + & - materialpoint_dPdF(3,3,3,3,1,ele))/3.0_pReal + matmul(HGMat,x_elem)*(homogenization_dPdF(1,1,1,1,1,ele) + & + homogenization_dPdF(2,2,2,2,1,ele) + & + homogenization_dPdF(3,3,3,3,1,ele))/3.0_pReal ctr = 0 do kk = 0, 1; do jj = 0, 1; do ii = 0, 1 ctr = ctr + 1 @@ -636,18 +636,18 @@ subroutine formJacobian(da_local,x_local,Jac_pre,Jac,dummy,ierr) row = col ele = ele + 1 K_ele = 0.0 - K_ele(1 :8 ,1 :8 ) = HGMat*(materialpoint_dPdF(1,1,1,1,1,ele) + & - materialpoint_dPdF(2,2,2,2,1,ele) + & - materialpoint_dPdF(3,3,3,3,1,ele))/3.0_pReal - K_ele(9 :16,9 :16) = HGMat*(materialpoint_dPdF(1,1,1,1,1,ele) + & - materialpoint_dPdF(2,2,2,2,1,ele) + & - materialpoint_dPdF(3,3,3,3,1,ele))/3.0_pReal - K_ele(17:24,17:24) = HGMat*(materialpoint_dPdF(1,1,1,1,1,ele) + & - materialpoint_dPdF(2,2,2,2,1,ele) + & - materialpoint_dPdF(3,3,3,3,1,ele))/3.0_pReal + K_ele(1 :8 ,1 :8 ) = HGMat*(homogenization_dPdF(1,1,1,1,1,ele) + & + homogenization_dPdF(2,2,2,2,1,ele) + & + homogenization_dPdF(3,3,3,3,1,ele))/3.0_pReal + K_ele(9 :16,9 :16) = HGMat*(homogenization_dPdF(1,1,1,1,1,ele) + & + homogenization_dPdF(2,2,2,2,1,ele) + & + homogenization_dPdF(3,3,3,3,1,ele))/3.0_pReal + K_ele(17:24,17:24) = HGMat*(homogenization_dPdF(1,1,1,1,1,ele) + & + homogenization_dPdF(2,2,2,2,1,ele) + & + homogenization_dPdF(3,3,3,3,1,ele))/3.0_pReal K_ele = K_ele + & matmul(transpose(BMatFull), & - matmul(reshape(reshape(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,ele), & + matmul(reshape(reshape(homogenization_dPdF(1:3,1:3,1:3,1:3,1,ele), & shape=[3,3,3,3], order=[2,1,4,3]),shape=[9,9]),BMatFull))*detJ call MatSetValuesStencil(Jac,24,row,24,col,K_ele,ADD_VALUES,ierr) CHKERRQ(ierr) diff --git a/src/grid/grid_mech_spectral_basic.f90 b/src/grid/grid_mech_spectral_basic.f90 index d2f6b40da..ec6c3a540 100644 --- a/src/grid/grid_mech_spectral_basic.f90 +++ b/src/grid/grid_mech_spectral_basic.f90 @@ -198,7 +198,7 @@ subroutine grid_mech_spectral_basic_init F = reshape(F_lastInc,[9,grid(1),grid(2),grid3]) endif restartRead - materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent + homogenization_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent call utilities_updateCoords(reshape(F,shape(F_lastInc))) call utilities_constitutiveResponse(P,temp33_Real,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2 reshape(F,shape(F_lastInc)), & ! target F @@ -324,7 +324,7 @@ subroutine grid_mech_spectral_basic_forward(cutBack,guess,timeinc,timeinc_old,lo rotation_BC%rotate(F_aimDot,active=.true.)) F_lastInc = reshape(F,[3,3,grid(1),grid(2),grid3]) - materialpoint_F0 = reshape(F, [3,3,1,product(grid(1:2))*grid3]) + homogenization_F0 = reshape(F, [3,3,1,product(grid(1:2))*grid3]) endif !-------------------------------------------------------------------------------------------------- diff --git a/src/grid/grid_mech_spectral_polarisation.f90 b/src/grid/grid_mech_spectral_polarisation.f90 index 3d495ddf0..8f9ea81b3 100644 --- a/src/grid/grid_mech_spectral_polarisation.f90 +++ b/src/grid/grid_mech_spectral_polarisation.f90 @@ -224,7 +224,7 @@ subroutine grid_mech_spectral_polarisation_init F_tau_lastInc = 2.0_pReal*F_lastInc endif restartRead - materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent + homogenization_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent call utilities_updateCoords(reshape(F,shape(F_lastInc))) call utilities_constitutiveResponse(P,temp33_Real,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2 reshape(F,shape(F_lastInc)), & ! target F @@ -364,7 +364,7 @@ subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,timeinc,timeinc F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid3]) F_tau_lastInc = reshape(F_tau,[3,3,grid(1),grid(2),grid3]) - materialpoint_F0 = reshape(F,[3,3,1,product(grid(1:2))*grid3]) + homogenization_F0 = reshape(F,[3,3,1,product(grid(1:2))*grid3]) endif !-------------------------------------------------------------------------------------------------- @@ -606,7 +606,7 @@ subroutine formResidual(in, FandF_tau, & do k = 1, grid3; do j = 1, grid(2); do i = 1, grid(1) e = e + 1 residual_F(1:3,1:3,i,j,k) = & - math_mul3333xx33(math_invSym3333(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,e) + C_scale), & + math_mul3333xx33(math_invSym3333(homogenization_dPdF(1:3,1:3,1:3,1:3,1,e) + C_scale), & residual_F(1:3,1:3,i,j,k) - matmul(F(1:3,1:3,i,j,k), & math_mul3333xx33(C_scale,F_tau(1:3,1:3,i,j,k) - F(1:3,1:3,i,j,k) - math_I3))) & + residual_F_tau(1:3,1:3,i,j,k) diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index 77047a317..fddd1885f 100644 --- a/src/grid/spectral_utilities.f90 +++ b/src/grid/spectral_utilities.f90 @@ -802,7 +802,7 @@ end subroutine utilities_fourierTensorDivergence !-------------------------------------------------------------------------------------------------- -!> @brief calculate constitutive response from materialpoint_F0 to F during timeinc +!> @brief calculate constitutive response from homogenization_F0 to F during timeinc !-------------------------------------------------------------------------------------------------- subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,& F,timeinc,rotation_BC) @@ -824,11 +824,11 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,& print'(/,a)', ' ... evaluating constitutive response ......................................' flush(IO_STDOUT) - materialpoint_F = reshape(F,[3,3,1,product(grid(1:2))*grid3]) ! set materialpoint target F to estimated field + homogenization_F = reshape(F,[3,3,1,product(grid(1:2))*grid3]) ! set materialpoint target F to estimated field call materialpoint_stressAndItsTangent(timeinc) ! calculate P field - P = reshape(materialpoint_P, [3,3,grid(1),grid(2),grid3]) + P = reshape(homogenization_P, [3,3,grid(1),grid(2),grid3]) P_av = sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt ! average of P call MPI_Allreduce(MPI_IN_PLACE,P_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) if (debugRotation) & @@ -845,13 +845,13 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,& dPdF_min = huge(1.0_pReal) dPdF_norm_min = huge(1.0_pReal) do i = 1, product(grid(1:2))*grid3 - if (dPdF_norm_max < sum(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,i)**2.0_pReal)) then - dPdF_max = materialpoint_dPdF(1:3,1:3,1:3,1:3,1,i) - dPdF_norm_max = sum(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,i)**2.0_pReal) + if (dPdF_norm_max < sum(homogenization_dPdF(1:3,1:3,1:3,1:3,1,i)**2.0_pReal)) then + dPdF_max = homogenization_dPdF(1:3,1:3,1:3,1:3,1,i) + dPdF_norm_max = sum(homogenization_dPdF(1:3,1:3,1:3,1:3,1,i)**2.0_pReal) endif - if (dPdF_norm_min > sum(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,i)**2.0_pReal)) then - dPdF_min = materialpoint_dPdF(1:3,1:3,1:3,1:3,1,i) - dPdF_norm_min = sum(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,i)**2.0_pReal) + if (dPdF_norm_min > sum(homogenization_dPdF(1:3,1:3,1:3,1:3,1,i)**2.0_pReal)) then + dPdF_min = homogenization_dPdF(1:3,1:3,1:3,1:3,1,i) + dPdF_norm_min = sum(homogenization_dPdF(1:3,1:3,1:3,1:3,1,i)**2.0_pReal) endif end do @@ -869,7 +869,7 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,& C_minmaxAvg = 0.5_pReal*(dPdF_max + dPdF_min) - C_volAvg = sum(sum(materialpoint_dPdF,dim=6),dim=5) + C_volAvg = sum(sum(homogenization_dPdF,dim=6),dim=5) call MPI_Allreduce(MPI_IN_PLACE,C_volAvg,81,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) C_volAvg = C_volAvg * wgt diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 75f93dc08..c717c24fd 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -31,12 +31,12 @@ module homogenization !-------------------------------------------------------------------------------------------------- ! General variables for the homogenization at a material point real(pReal), dimension(:,:,:,:), allocatable, public :: & - materialpoint_F0, & !< def grad of IP at start of FE increment - materialpoint_F !< def grad of IP to be reached at end of FE increment + homogenization_F0, & !< def grad of IP at start of FE increment + homogenization_F !< def grad of IP to be reached at end of FE increment real(pReal), dimension(:,:,:,:), allocatable, public, protected :: & - materialpoint_P !< first P--K stress of IP + homogenization_P !< first P--K stress of IP real(pReal), dimension(:,:,:,:,:,:), allocatable, public, protected :: & - materialpoint_dPdF !< tangent of first P--K stress at IP + homogenization_dPdF !< tangent of first P--K stress at IP type :: tNumerics integer :: & @@ -181,10 +181,10 @@ subroutine homogenization_init !-------------------------------------------------------------------------------------------------- ! allocate and initialize global variables - allocate(materialpoint_dPdF(3,3,3,3,discretization_nIP,discretization_nElem), source=0.0_pReal) - materialpoint_F0 = spread(spread(math_I3,3,discretization_nIP),4,discretization_nElem) ! initialize to identity - materialpoint_F = materialpoint_F0 ! initialize to identity - allocate(materialpoint_P(3,3,discretization_nIP,discretization_nElem), source=0.0_pReal) + allocate(homogenization_dPdF(3,3,3,3,discretization_nIP,discretization_nElem), source=0.0_pReal) + homogenization_F0 = spread(spread(math_I3,3,discretization_nIP),4,discretization_nElem) ! initialize to identity + homogenization_F = homogenization_F0 ! initialize to identity + allocate(homogenization_P(3,3,discretization_nIP,discretization_nElem), source=0.0_pReal) print'(/,a)', ' <<<+- homogenization init -+>>>'; flush(IO_STDOUT) @@ -330,8 +330,8 @@ subroutine materialpoint_stressAndItsTangent(dt) myNgrains = homogenization_Nconstituent(material_homogenizationAt(e)) IpLooping2: do i = FEsolving_execIP(1),FEsolving_execIP(2) if(requested(i,e) .and. .not. doneAndHappy(1,i,e)) then ! requested but not yet done - call partitionDeformation(materialpoint_F0(1:3,1:3,i,e) & - + (materialpoint_F(1:3,1:3,i,e)-materialpoint_F0(1:3,1:3,i,e))& + call partitionDeformation(homogenization_F0(1:3,1:3,i,e) & + + (homogenization_F(1:3,1:3,i,e)-homogenization_F0(1:3,1:3,i,e))& *(subStep(i,e)+subFrac(i,e)), & i,e) crystallite_dt(1:myNgrains,i,e) = dt*subStep(i,e) ! propagate materialpoint dt to grains @@ -357,8 +357,8 @@ subroutine materialpoint_stressAndItsTangent(dt) doneAndHappy(1:2,i,e) = [.true.,.false.] else doneAndHappy(1:2,i,e) = updateState(dt*subStep(i,e), & - materialpoint_F0(1:3,1:3,i,e) & - + (materialpoint_F(1:3,1:3,i,e)-materialpoint_F0(1:3,1:3,i,e)) & + homogenization_F0(1:3,1:3,i,e) & + + (homogenization_F(1:3,1:3,i,e)-homogenization_F0(1:3,1:3,i,e)) & *(subStep(i,e)+subFrac(i,e)), & i,e) converged(i,e) = all(doneAndHappy(1:2,i,e)) ! converged if done and happy @@ -492,16 +492,16 @@ subroutine averageStressAndItsTangent(ip,el) chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el))) case (HOMOGENIZATION_NONE_ID) chosenHomogenization - materialpoint_P(1:3,1:3,ip,el) = crystallite_P(1:3,1:3,1,ip,el) - materialpoint_dPdF(1:3,1:3,1:3,1:3,ip,el) = crystallite_stressTangent(1,ip,el) + homogenization_P(1:3,1:3,ip,el) = crystallite_P(1:3,1:3,1,ip,el) + homogenization_dPdF(1:3,1:3,1:3,1:3,ip,el) = crystallite_stressTangent(1,ip,el) case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization do c = 1, homogenization_Nconstituent(material_homogenizationAt(el)) dPdFs(:,:,:,:,c) = crystallite_stressTangent(c,ip,el) enddo call mech_isostrain_averageStressAndItsTangent(& - materialpoint_P(1:3,1:3,ip,el), & - materialpoint_dPdF(1:3,1:3,1:3,1:3,ip,el),& + homogenization_P(1:3,1:3,ip,el), & + homogenization_dPdF(1:3,1:3,1:3,1:3,ip,el),& crystallite_P(1:3,1:3,1:homogenization_Nconstituent(material_homogenizationAt(el)),ip,el), & dPdFs, & homogenization_typeInstance(material_homogenizationAt(el))) @@ -511,8 +511,8 @@ subroutine averageStressAndItsTangent(ip,el) dPdFs(:,:,:,:,c) = crystallite_stressTangent(c,ip,el) enddo call mech_RGC_averageStressAndItsTangent(& - materialpoint_P(1:3,1:3,ip,el), & - materialpoint_dPdF(1:3,1:3,1:3,1:3,ip,el),& + homogenization_P(1:3,1:3,ip,el), & + homogenization_dPdF(1:3,1:3,1:3,1:3,ip,el),& crystallite_P(1:3,1:3,1:homogenization_Nconstituent(material_homogenizationAt(el)),ip,el), & dPdFs, & homogenization_typeInstance(material_homogenizationAt(el))) @@ -539,10 +539,10 @@ subroutine homogenization_results group = trim(group_base)//'/generic' call results_closeGroup(results_addGroup(group)) - !temp = reshape(materialpoint_F,[3,3,discretization_nIP*discretization_nElem]) + !temp = reshape(homogenization_F,[3,3,discretization_nIP*discretization_nElem]) !call results_writeDataset(group,temp,'F',& ! 'deformation gradient','1') - !temp = reshape(materialpoint_P,[3,3,discretization_nIP*discretization_nElem]) + !temp = reshape(homogenization_P,[3,3,discretization_nIP*discretization_nElem]) !call results_writeDataset(group,temp,'P',& ! '1st Piola-Kirchhoff stress','Pa') diff --git a/src/material.f90 b/src/material.f90 index 7197506f4..a9459e42d 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -177,7 +177,7 @@ subroutine material_init(restart) if (.not. restart) then call results_openJobFile call results_mapping_constituent(material_phaseAt,material_phaseMemberAt,material_name_phase) - call results_mapping_materialpoint(material_homogenizationAt,material_homogenizationMemberAt,material_name_homogenization) + call results_mapping_homogenization(material_homogenizationAt,material_homogenizationMemberAt,material_name_homogenization) call results_closeJobFile endif diff --git a/src/mesh/FEM_utilities.f90 b/src/mesh/FEM_utilities.f90 index 4927d0c1c..4c958ee2e 100644 --- a/src/mesh/FEM_utilities.f90 +++ b/src/mesh/FEM_utilities.f90 @@ -164,7 +164,7 @@ subroutine utilities_constitutiveResponse(timeinc,P_av,forwardData) cutBack = .false. ! reset cutBack status - P_av = sum(sum(materialpoint_P,dim=4),dim=3) * wgt ! average of P + P_av = sum(sum(homogenization_P,dim=4),dim=3) * wgt ! average of P call MPI_Allreduce(MPI_IN_PLACE,P_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) end subroutine utilities_constitutiveResponse diff --git a/src/mesh/mesh_mech_FEM.f90 b/src/mesh/mesh_mech_FEM.f90 index f7d33adcf..8aa084ac8 100644 --- a/src/mesh/mesh_mech_FEM.f90 +++ b/src/mesh/mesh_mech_FEM.f90 @@ -400,15 +400,15 @@ subroutine FEM_mech_formResidual(dm_local,xx_local,f_local,dummy,ierr) (((qPt*nBasis + basis)*dimPlex + comp)*dimPlex+comp+1)*dimPlex)) enddo enddo - materialpoint_F(1:dimPlex,1:dimPlex,qPt+1,cell+1) = & + homogenization_F(1:dimPlex,1:dimPlex,qPt+1,cell+1) = & reshape(matmul(BMat,x_scal),shape=[dimPlex,dimPlex], order=[2,1]) enddo if (num%BBarStabilisation) then - detFAvg = math_det33(sum(materialpoint_F(1:3,1:3,1:nQuadrature,cell+1),dim=3)/real(nQuadrature)) + detFAvg = math_det33(sum(homogenization_F(1:3,1:3,1:nQuadrature,cell+1),dim=3)/real(nQuadrature)) do qPt = 1, nQuadrature - materialpoint_F(1:dimPlex,1:dimPlex,qPt,cell+1) = & - materialpoint_F(1:dimPlex,1:dimPlex,qPt,cell+1)* & - (detFAvg/math_det33(materialpoint_F(1:3,1:3,qPt,cell+1)))**(1.0/real(dimPlex)) + homogenization_F(1:dimPlex,1:dimPlex,qPt,cell+1) = & + homogenization_F(1:dimPlex,1:dimPlex,qPt,cell+1)* & + (detFAvg/math_det33(homogenization_F(1:3,1:3,qPt,cell+1)))**(1.0/real(dimPlex)) enddo endif @@ -443,7 +443,7 @@ subroutine FEM_mech_formResidual(dm_local,xx_local,f_local,dummy,ierr) enddo f_scal = f_scal + & matmul(transpose(BMat), & - reshape(transpose(materialpoint_P(1:dimPlex,1:dimPlex,qPt+1,cell+1)), & + reshape(transpose(homogenization_P(1:dimPlex,1:dimPlex,qPt+1,cell+1)), & shape=[dimPlex*dimPlex]))*qWeights(qPt+1) enddo f_scal = f_scal*abs(detJ) @@ -545,7 +545,7 @@ subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr) (((qPt*nBasis + basis)*dimPlex + comp)*dimPlex+comp+1)*dimPlex)) enddo enddo - MatA = matmul(reshape(reshape(materialpoint_dPdF(1:dimPlex,1:dimPlex,1:dimPlex,1:dimPlex,qPt+1,cell+1), & + MatA = matmul(reshape(reshape(homogenization_dPdF(1:dimPlex,1:dimPlex,1:dimPlex,1:dimPlex,qPt+1,cell+1), & shape=[dimPlex,dimPlex,dimPlex,dimPlex], order=[2,1,4,3]), & shape=[dimPlex*dimPlex,dimPlex*dimPlex]),BMat)*qWeights(qPt+1) if (num%BBarStabilisation) then @@ -553,12 +553,12 @@ subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr) FInv = math_inv33(F) K_eA = K_eA + matmul(transpose(BMat),MatA)*math_det33(FInv)**(1.0/real(dimPlex)) K_eB = K_eB - & - matmul(transpose(matmul(reshape(materialpoint_F(1:dimPlex,1:dimPlex,qPt+1,cell+1), & + matmul(transpose(matmul(reshape(homogenization_F(1:dimPlex,1:dimPlex,qPt+1,cell+1), & shape=[dimPlex*dimPlex,1]), & matmul(reshape(FInv(1:dimPlex,1:dimPlex), & shape=[1,dimPlex*dimPlex],order=[2,1]),BMat))),MatA) MatB = MatB + & - matmul(reshape(materialpoint_F(1:dimPlex,1:dimPlex,qPt+1,cell+1),shape=[1,dimPlex*dimPlex]),MatA) + matmul(reshape(homogenization_F(1:dimPlex,1:dimPlex,qPt+1,cell+1),shape=[1,dimPlex*dimPlex]),MatA) FAvg = FAvg + F BMatAvg = BMatAvg + BMat else @@ -630,7 +630,7 @@ subroutine FEM_mech_forward(guess,timeinc,timeinc_old,fieldBC) ! forward last inc if (guess .and. .not. cutBack) then ForwardData = .True. - materialpoint_F0 = materialpoint_F + homogenization_F0 = homogenization_F call SNESGetDM(mech_snes,dm_local,ierr); CHKERRQ(ierr) !< retrieve mesh info from mech_snes into dm_local call DMGetSection(dm_local,section,ierr); CHKERRQ(ierr) call DMGetLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr) diff --git a/src/results.f90 b/src/results.f90 index aec90d7be..1ccc6bfab 100644 --- a/src/results.f90 +++ b/src/results.f90 @@ -56,7 +56,7 @@ module results results_addAttribute, & results_removeLink, & results_mapping_constituent, & - results_mapping_materialpoint + results_mapping_homogenization contains subroutine results_init(restart) @@ -644,7 +644,7 @@ end subroutine results_mapping_constituent !-------------------------------------------------------------------------------------------------- !> @brief adds the unique mapping from spatial position and constituent ID to results !-------------------------------------------------------------------------------------------------- -subroutine results_mapping_materialpoint(homogenizationAt,memberAtLocal,label) +subroutine results_mapping_homogenization(homogenizationAt,memberAtLocal,label) integer, dimension(:), intent(in) :: homogenizationAt !< homogenization section at (element) integer, dimension(:,:), intent(in) :: memberAtLocal !< homogenization member at (IP,element) @@ -711,13 +711,13 @@ subroutine results_mapping_materialpoint(homogenizationAt,memberAtLocal,label) ! MPI settings and communication #ifdef PETSc call h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, ierr) - if (ierr < 0) call IO_error(1,ext_msg='results_mapping_materialpoint: h5pset_dxpl_mpio_f') + if (ierr < 0) call IO_error(1,ext_msg='results_mapping_homogenization: h5pset_dxpl_mpio_f') call MPI_allreduce(MPI_IN_PLACE,writeSize,worldsize,MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr) ! get output at each process - if (ierr /= 0) call IO_error(894,ext_msg='results_mapping_materialpoint: MPI_allreduce/writeSize') + if (ierr /= 0) call IO_error(894,ext_msg='results_mapping_homogenization: MPI_allreduce/writeSize') call MPI_allreduce(MPI_IN_PLACE,memberOffset,size(memberOffset),MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr)! get offset at each process - if (ierr /= 0) call IO_error(894,ext_msg='results_mapping_materialpoint: MPI_allreduce/memberOffset') + if (ierr /= 0) call IO_error(894,ext_msg='results_mapping_homogenization: MPI_allreduce/memberOffset') #endif myShape = int([writeSize(worldrank)], HSIZE_T) @@ -727,13 +727,13 @@ subroutine results_mapping_materialpoint(homogenizationAt,memberAtLocal,label) !-------------------------------------------------------------------------------------------------- ! create dataspace in memory (local shape = hyperslab) and in file (global shape) call h5screate_simple_f(1,myShape,memspace_id,ierr,myShape) - if (ierr < 0) call IO_error(1,ext_msg='results_mapping_materialpoint: h5screate_simple_f/memspace_id') + if (ierr < 0) call IO_error(1,ext_msg='results_mapping_homogenization: h5screate_simple_f/memspace_id') call h5screate_simple_f(1,totalShape,filespace_id,ierr,totalShape) - if (ierr < 0) call IO_error(1,ext_msg='results_mapping_materialpoint: h5screate_simple_f/filespace_id') + if (ierr < 0) call IO_error(1,ext_msg='results_mapping_homogenization: h5screate_simple_f/filespace_id') call h5sselect_hyperslab_f(filespace_id, H5S_SELECT_SET_F, myOffset, myShape, ierr) - if (ierr < 0) call IO_error(1,ext_msg='results_mapping_materialpoint: h5sselect_hyperslab_f') + if (ierr < 0) call IO_error(1,ext_msg='results_mapping_homogenization: h5sselect_hyperslab_f') !--------------------------------------------------------------------------------------------------- ! expand phaseAt to consider IPs (is not stored per IP) @@ -753,14 +753,14 @@ subroutine results_mapping_materialpoint(homogenizationAt,memberAtLocal,label) loc_id = results_openGroup('/mapping') call h5dcreate_f(loc_id, 'homogenization', dtype_id, filespace_id, dset_id, ierr) - if (ierr < 0) call IO_error(1,ext_msg='results_mapping_materialpoint: h5dcreate_f') + if (ierr < 0) call IO_error(1,ext_msg='results_mapping_homogenization: h5dcreate_f') call h5dwrite_f(dset_id, name_id, reshape(label(pack(homogenizationAtMaterialpoint,.true.)),myShape), & myShape, ierr, file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id) - if (ierr < 0) call IO_error(1,ext_msg='results_mapping_materialpoint: h5dwrite_f/name_id') + if (ierr < 0) call IO_error(1,ext_msg='results_mapping_homogenization: h5dwrite_f/name_id') call h5dwrite_f(dset_id, position_id, reshape(pack(memberAtGlobal,.true.),myShape), & myShape, ierr, file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id) - if (ierr < 0) call IO_error(1,ext_msg='results_mapping_materialpoint: h5dwrite_f/position_id') + if (ierr < 0) call IO_error(1,ext_msg='results_mapping_homogenization: h5dwrite_f/position_id') !-------------------------------------------------------------------------------------------------- ! close all @@ -776,7 +776,7 @@ subroutine results_mapping_materialpoint(homogenizationAt,memberAtLocal,label) ! for backward compatibility call results_setLink('/mapping/homogenization','/mapping/cellResults/materialpoint') -end subroutine results_mapping_materialpoint +end subroutine results_mapping_homogenization !--------------------------------------------------------------------------------------------------