diff --git a/src/CPFEM.f90 b/src/CPFEM.f90 index 76522fb16..a19a70432 100644 --- a/src/CPFEM.f90 +++ b/src/CPFEM.f90 @@ -153,7 +153,7 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS H integer(pInt) elCP, & ! crystal plasticity element number - i, j, k, l, m, n, ph, homog, mySource + i, j, k, l, m, n, ph, homog, mySource,ma real(pReal), parameter :: ODD_STRESS = 1e15_pReal, & !< return value for stress if terminallyIll ODD_JACOBIAN = 1e50_pReal !< return value for jacobian if terminallyIll @@ -161,6 +161,8 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS elCP = mesh_FEM2DAMASK_elem(elFE) + ma = (elCP-1) * discretization_nIPs + ip + if (debugCPFEM%basic .and. elCP == debugCPFEM%element .and. ip == debugCPFEM%ip) then print'(/,a)', '#############################################' print'(a1,a22,1x,i8,a13)', '#','element', elCP, '#' @@ -184,8 +186,8 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS temperature(material_homogenizationAt(elCP))%p(material_homogenizationMemberAt(ip,elCP)) = & temperature_inp end select chosenThermal1 - homogenization_F0(1:3,1:3,ip,elCP) = ffn - homogenization_F(1:3,1:3,ip,elCP) = ffn1 + homogenization_F0(1:3,1:3,ma) = ffn + homogenization_F(1:3,1:3,ma) = ffn1 if (iand(mode, CPFEM_CALCRESULTS) /= 0_pInt) then @@ -212,17 +214,17 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS else terminalIllness ! translate from P to sigma - 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)) + Kirchhoff = matmul(homogenization_P(1:3,1:3,ma), transpose(homogenization_F(1:3,1:3,ma))) + J_inverse = 1.0_pReal / math_det33(homogenization_F(1:3,1:3,ma)) 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) & - + 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) & + + homogenization_F(j,m,ma) * homogenization_F(l,n,ma) & + * homogenization_dPdF(i,m,k,n,ma) & + - math_delta(j,l) * homogenization_F(i,m,ma) * homogenization_P(k,m,ma) & + 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 741ce404a..cdf806b35 100644 --- a/src/grid/grid_mech_FEM.f90 +++ b/src/grid/grid_mech_FEM.f90 @@ -238,7 +238,7 @@ subroutine grid_mech_FEM_init F = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) endif restartRead - homogenization_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,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent call utilities_updateCoords(F) call utilities_constitutiveResponse(P_current,P_av,C_volAvg,devNull, & ! stress field, stress avg, global average of stiffness and (min+max)/2 F, & ! target F @@ -359,7 +359,7 @@ subroutine grid_mech_FEM_forward(cutBack,guess,Delta_t,Delta_t_old,t_remaining,& F_lastInc = F - homogenization_F0 = reshape(F, [3,3,1,product(grid(1:2))*grid3]) + homogenization_F0 = reshape(F, [3,3,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)*(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 + matmul(HGMat,x_elem)*(homogenization_dPdF(1,1,1,1,ele) + & + homogenization_dPdF(2,2,2,2,ele) + & + homogenization_dPdF(3,3,3,3,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*(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(1 :8 ,1 :8 ) = HGMat*(homogenization_dPdF(1,1,1,1,ele) + & + homogenization_dPdF(2,2,2,2,ele) + & + homogenization_dPdF(3,3,3,3,ele))/3.0_pReal + K_ele(9 :16,9 :16) = HGMat*(homogenization_dPdF(1,1,1,1,ele) + & + homogenization_dPdF(2,2,2,2,ele) + & + homogenization_dPdF(3,3,3,3,ele))/3.0_pReal + K_ele(17:24,17:24) = HGMat*(homogenization_dPdF(1,1,1,1,ele) + & + homogenization_dPdF(2,2,2,2,ele) + & + homogenization_dPdF(3,3,3,3,ele))/3.0_pReal K_ele = K_ele + & matmul(transpose(BMatFull), & - matmul(reshape(reshape(homogenization_dPdF(1:3,1:3,1:3,1:3,1,ele), & + matmul(reshape(reshape(homogenization_dPdF(1:3,1:3,1:3,1:3,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 d87a22fb7..ebaaf3b55 100644 --- a/src/grid/grid_mech_spectral_basic.f90 +++ b/src/grid/grid_mech_spectral_basic.f90 @@ -199,7 +199,7 @@ subroutine grid_mech_spectral_basic_init F = reshape(F_lastInc,[9,grid(1),grid(2),grid3]) endif restartRead - homogenization_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,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent call utilities_updateCoords(reshape(F,shape(F_lastInc))) call utilities_constitutiveResponse(P,P_av,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2 reshape(F,shape(F_lastInc)), & ! target F @@ -319,7 +319,7 @@ subroutine grid_mech_spectral_basic_forward(cutBack,guess,Delta_t,Delta_t_old,t_ rotation_BC%rotate(F_aimDot,active=.true.)) F_lastInc = reshape(F,[3,3,grid(1),grid(2),grid3]) - homogenization_F0 = reshape(F,[3,3,1,product(grid(1:2))*grid3]) + homogenization_F0 = reshape(F,[3,3,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 9bbb40a53..9f2a17c97 100644 --- a/src/grid/grid_mech_spectral_polarisation.f90 +++ b/src/grid/grid_mech_spectral_polarisation.f90 @@ -225,7 +225,7 @@ subroutine grid_mech_spectral_polarisation_init F_tau_lastInc = 2.0_pReal*F_lastInc endif restartRead - homogenization_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,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent call utilities_updateCoords(reshape(F,shape(F_lastInc))) call utilities_constitutiveResponse(P,P_av,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2 reshape(F,shape(F_lastInc)), & ! target F @@ -359,7 +359,7 @@ subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,Delta_t,Delta_t F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid3]) F_tau_lastInc = reshape(F_tau,[3,3,grid(1),grid(2),grid3]) - homogenization_F0 = reshape(F,[3,3,1,product(grid(1:2))*grid3]) + homogenization_F0 = reshape(F,[3,3,product(grid(1:2))*grid3]) endif !-------------------------------------------------------------------------------------------------- @@ -604,7 +604,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(homogenization_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,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 989448dc3..c0c84233d 100644 --- a/src/grid/spectral_utilities.f90 +++ b/src/grid/spectral_utilities.f90 @@ -810,7 +810,7 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,& print'(/,a)', ' ... evaluating constitutive response ......................................' flush(IO_STDOUT) - homogenization_F = reshape(F,[3,3,1,product(grid(1:2))*grid3]) ! set materialpoint target F to estimated field + homogenization_F = reshape(F,[3,3,product(grid(1:2))*grid3]) ! set materialpoint target F to estimated field call materialpoint_stressAndItsTangent(timeinc) ! calculate P field @@ -829,13 +829,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(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) + if (dPdF_norm_max < sum(homogenization_dPdF(1:3,1:3,1:3,1:3,i)**2.0_pReal)) then + dPdF_max = homogenization_dPdF(1:3,1:3,1:3,1:3,i) + dPdF_norm_max = sum(homogenization_dPdF(1:3,1:3,1:3,1:3,i)**2.0_pReal) endif - 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) + if (dPdF_norm_min > sum(homogenization_dPdF(1:3,1:3,1:3,1:3,i)**2.0_pReal)) then + dPdF_min = homogenization_dPdF(1:3,1:3,1:3,1:3,i) + dPdF_norm_min = sum(homogenization_dPdF(1:3,1:3,1:3,1:3,i)**2.0_pReal) endif end do @@ -853,7 +853,7 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,& C_minmaxAvg = 0.5_pReal*(dPdF_max + dPdF_min) - C_volAvg = sum(sum(homogenization_dPdF,dim=6),dim=5) + C_volAvg = sum(homogenization_dPdF,dim=5) call MPI_Allreduce(MPI_IN_PLACE,C_volAvg,81,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) if (ierr /= 0) error stop 'MPI error' C_volAvg = C_volAvg * wgt diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 5958f35fa..57478e039 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -30,12 +30,12 @@ module homogenization !-------------------------------------------------------------------------------------------------- ! General variables for the homogenization at a material point - real(pReal), dimension(:,:,:,:), allocatable, public :: & + real(pReal), dimension(:,:,:), allocatable, public :: & 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 :: & ! Issue with ifort + real(pReal), dimension(:,:,:), allocatable, public :: & !, protected :: & Issue with ifort homogenization_P !< first P--K stress of IP - real(pReal), dimension(:,:,:,:,:,:), allocatable, public :: & !, protected :: & + real(pReal), dimension(:,:,:,:,:), allocatable, public :: & !, protected :: & homogenization_dPdF !< tangent of first P--K stress at IP @@ -193,6 +193,7 @@ subroutine materialpoint_stressAndItsTangent(dt) converged logical, dimension(2,discretization_nIPs,discretization_Nelems) :: & doneAndHappy + integer :: m !-------------------------------------------------------------------------------------------------- @@ -227,7 +228,7 @@ subroutine materialpoint_stressAndItsTangent(dt) any(subStep(FEsolving_execIP(1):FEsolving_execIP(2),& FEsolving_execElem(1):FEsolving_execElem(2)) > num%subStepMinHomog)) - !$OMP PARALLEL DO + !$OMP PARALLEL DO PRIVATE(m) elementLooping1: do e = FEsolving_execElem(1),FEsolving_execElem(2) myNgrains = homogenization_Nconstituents(material_homogenizationAt(e)) IpLooping1: do i = FEsolving_execIP(1),FEsolving_execIP(2) @@ -297,13 +298,14 @@ subroutine materialpoint_stressAndItsTangent(dt) !-------------------------------------------------------------------------------------------------- ! deformation partitioning - !$OMP PARALLEL DO PRIVATE(myNgrains) + !$OMP PARALLEL DO PRIVATE(myNgrains,m) elementLooping2: do e = FEsolving_execElem(1),FEsolving_execElem(2) myNgrains = homogenization_Nconstituents(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 mech_partition(homogenization_F0(1:3,1:3,i,e) & - + (homogenization_F(1:3,1:3,i,e)-homogenization_F0(1:3,1:3,i,e))& + m = (e-1)*discretization_nIPs + i + call mech_partition(homogenization_F0(1:3,1:3,m) & + + (homogenization_F(1:3,1:3,m)-homogenization_F0(1:3,1:3,m))& *(subStep(i,e)+subFrac(i,e)), & i,e) crystallite_dt(1:myNgrains,i,e) = dt*subStep(i,e) ! propagate materialpoint dt to grains @@ -321,16 +323,17 @@ subroutine materialpoint_stressAndItsTangent(dt) !-------------------------------------------------------------------------------------------------- ! state update - !$OMP PARALLEL DO + !$OMP PARALLEL DO PRIVATE(m) elementLooping3: do e = FEsolving_execElem(1),FEsolving_execElem(2) IpLooping3: do i = FEsolving_execIP(1),FEsolving_execIP(2) if (requested(i,e) .and. .not. doneAndHappy(1,i,e)) then if (.not. converged(i,e)) then doneAndHappy(1:2,i,e) = [.true.,.false.] else + m = (e-1)*discretization_nIPs + i doneAndHappy(1:2,i,e) = updateState(dt*subStep(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)) & + homogenization_F0(1:3,1:3,m) & + + (homogenization_F(1:3,1:3,m)-homogenization_F0(1:3,1:3,m)) & *(subStep(i,e)+subFrac(i,e)), & i,e) converged(i,e) = all(doneAndHappy(1:2,i,e)) ! converged if done and happy diff --git a/src/homogenization_mech.f90 b/src/homogenization_mech.f90 index 40d26df51..b0641be07 100644 --- a/src/homogenization_mech.f90 +++ b/src/homogenization_mech.f90 @@ -73,10 +73,10 @@ module subroutine mech_init(num_homog) print'(/,a)', ' <<<+- homogenization_mech init -+>>>' - allocate(homogenization_dPdF(3,3,3,3,discretization_nIPs,discretization_Nelems), source=0.0_pReal) - homogenization_F0 = spread(spread(math_I3,3,discretization_nIPs),4,discretization_Nelems) ! initialize to identity - homogenization_F = homogenization_F0 ! initialize to identity - allocate(homogenization_P(3,3,discretization_nIPs,discretization_Nelems), source=0.0_pReal) + allocate(homogenization_dPdF(3,3,3,3,discretization_nIPs*discretization_Nelems), source=0.0_pReal) + homogenization_F0 = spread(math_I3,3,discretization_nIPs*discretization_Nelems) ! initialize to identity + homogenization_F = homogenization_F0 ! initialize to identity + allocate(homogenization_P(3,3,discretization_nIPs*discretization_Nelems), source=0.0_pReal) num_homogMech => num_homog%get('mech',defaultVal=emptyDict) if (any(homogenization_type == HOMOGENIZATION_NONE_ID)) call mech_none_init @@ -127,23 +127,24 @@ module subroutine mech_homogenize(ip,el) integer, intent(in) :: & ip, & !< integration point el !< element number - integer :: c + integer :: c,m real(pReal) :: dPdFs(3,3,3,3,homogenization_Nconstituents(material_homogenizationAt(el))) + m = (el-1)* discretization_nIPs + ip chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el))) case (HOMOGENIZATION_NONE_ID) chosenHomogenization - 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) + homogenization_P(1:3,1:3,m) = crystallite_P(1:3,1:3,1,ip,el) + homogenization_dPdF(1:3,1:3,1:3,1:3,m) = crystallite_stressTangent(1,ip,el) case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization do c = 1, homogenization_Nconstituents(material_homogenizationAt(el)) dPdFs(:,:,:,:,c) = crystallite_stressTangent(c,ip,el) enddo call mech_isostrain_averageStressAndItsTangent(& - homogenization_P(1:3,1:3,ip,el), & - homogenization_dPdF(1:3,1:3,1:3,1:3,ip,el),& + homogenization_P(1:3,1:3,m), & + homogenization_dPdF(1:3,1:3,1:3,1:3,m),& crystallite_P(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), & dPdFs, & homogenization_typeInstance(material_homogenizationAt(el))) @@ -153,8 +154,8 @@ module subroutine mech_homogenize(ip,el) dPdFs(:,:,:,:,c) = crystallite_stressTangent(c,ip,el) enddo call mech_RGC_averageStressAndItsTangent(& - homogenization_P(1:3,1:3,ip,el), & - homogenization_dPdF(1:3,1:3,1:3,1:3,ip,el),& + homogenization_P(1:3,1:3,m), & + homogenization_dPdF(1:3,1:3,1:3,1:3,m),& crystallite_P(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), & dPdFs, & homogenization_typeInstance(material_homogenizationAt(el))) diff --git a/src/material.f90 b/src/material.f90 index 574da0d51..bb5f484f6 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -140,7 +140,6 @@ contains subroutine material_init(restart) logical, intent(in) :: restart - integer :: myHomog print'(/,a)', ' <<<+- material init -+>>>'; flush(IO_STDOUT) diff --git a/src/mesh/FEM_utilities.f90 b/src/mesh/FEM_utilities.f90 index 118735e89..cb81f1f0c 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(homogenization_P,dim=4),dim=3) * wgt ! average of P + P_av = sum(homogenization_P,dim=3) * wgt 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 eb5f862c2..e19c35998 100644 --- a/src/mesh/mesh_mech_FEM.f90 +++ b/src/mesh/mesh_mech_FEM.f90 @@ -316,16 +316,16 @@ subroutine FEM_mech_formResidual(dm_local,xx_local,f_local,dummy,ierr) Vec :: x_local, f_local, xx_local PetscSection :: section PetscScalar, dimension(:), pointer :: x_scal, pf_scal - PetscScalar, target :: f_scal(cellDof) - PetscReal :: detJ, IcellJMat(dimPlex,dimPlex) - PetscReal, pointer,dimension(:) :: pV0, pCellJ, pInvcellJ, basisField, basisFieldDer + PetscScalar, dimension(cellDof), target :: f_scal + PetscReal :: IcellJMat(dimPlex,dimPlex) + PetscReal, dimension(:),pointer :: pV0, pCellJ, pInvcellJ, basisField, basisFieldDer PetscInt :: cellStart, cellEnd, cell, field, face, & qPt, basis, comp, cidx, & - numFields - PetscReal :: detFAvg - PetscReal :: BMat(dimPlex*dimPlex,cellDof) + numFields, & + bcSize,m + PetscReal :: detFAvg, detJ + PetscReal, dimension(dimPlex*dimPlex,cellDof) :: BMat - PetscInt :: bcSize IS :: bcPoints @@ -366,6 +366,7 @@ subroutine FEM_mech_formResidual(dm_local,xx_local,f_local,dummy,ierr) CHKERRQ(ierr) IcellJMat = reshape(pInvcellJ,shape=[dimPlex,dimPlex]) do qPt = 0, nQuadrature-1 + m = cell*nQuadrature + qPt+1 BMat = 0.0 do basis = 0, nBasis-1 do comp = 0, dimPlex-1 @@ -375,15 +376,14 @@ subroutine FEM_mech_formResidual(dm_local,xx_local,f_local,dummy,ierr) (((qPt*nBasis + basis)*dimPlex + comp)*dimPlex+comp+1)*dimPlex)) enddo enddo - homogenization_F(1:dimPlex,1:dimPlex,qPt+1,cell+1) = & - reshape(matmul(BMat,x_scal),shape=[dimPlex,dimPlex], order=[2,1]) + homogenization_F(1:dimPlex,1:dimPlex,m) = reshape(matmul(BMat,x_scal),shape=[dimPlex,dimPlex], order=[2,1]) enddo if (num%BBarStabilisation) then - detFAvg = math_det33(sum(homogenization_F(1:3,1:3,1:nQuadrature,cell+1),dim=3)/real(nQuadrature)) - do qPt = 1, nQuadrature - 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)) + detFAvg = math_det33(sum(homogenization_F(1:3,1:3,cell*nQuadrature+1:(cell+1)*nQuadrature),dim=3)/real(nQuadrature)) + do qPt = 0, nQuadrature-1 + m = cell*nQuadrature + qPt+1 + homogenization_F(1:dimPlex,1:dimPlex,m) = homogenization_F(1:dimPlex,1:dimPlex,m) & + * (detFAvg/math_det33(homogenization_F(1:3,1:3,m)))**(1.0/real(dimPlex)) enddo endif @@ -407,6 +407,7 @@ subroutine FEM_mech_formResidual(dm_local,xx_local,f_local,dummy,ierr) IcellJMat = reshape(pInvcellJ,shape=[dimPlex,dimPlex]) f_scal = 0.0 do qPt = 0, nQuadrature-1 + m = cell*nQuadrature + qPt+1 BMat = 0.0 do basis = 0, nBasis-1 do comp = 0, dimPlex-1 @@ -418,7 +419,7 @@ subroutine FEM_mech_formResidual(dm_local,xx_local,f_local,dummy,ierr) enddo f_scal = f_scal + & matmul(transpose(BMat), & - reshape(transpose(homogenization_P(1:dimPlex,1:dimPlex,qPt+1,cell+1)), & + reshape(transpose(homogenization_P(1:dimPlex,1:dimPlex,m)), & shape=[dimPlex*dimPlex]))*qWeights(qPt+1) enddo f_scal = f_scal*abs(detJ) @@ -463,7 +464,7 @@ subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr) K_eB PetscInt :: cellStart, cellEnd, cell, field, face, & - qPt, basis, comp, cidx,bcSize + qPt, basis, comp, cidx,bcSize, m IS :: bcPoints @@ -506,6 +507,7 @@ subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr) FAvg = 0.0 BMatAvg = 0.0 do qPt = 0, nQuadrature-1 + m = cell*nQuadrature + qPt + 1 BMat = 0.0 do basis = 0, nBasis-1 do comp = 0, dimPlex-1 @@ -516,7 +518,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(homogenization_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,m), & shape=[dimPlex,dimPlex,dimPlex,dimPlex], order=[2,1,4,3]), & shape=[dimPlex*dimPlex,dimPlex*dimPlex]),BMat)*qWeights(qPt+1) if (num%BBarStabilisation) then @@ -524,12 +526,11 @@ 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(homogenization_F(1:dimPlex,1:dimPlex,qPt+1,cell+1), & - shape=[dimPlex*dimPlex,1]), & + matmul(transpose(matmul(reshape(homogenization_F(1:dimPlex,1:dimPlex,m),shape=[dimPlex*dimPlex,1]), & matmul(reshape(FInv(1:dimPlex,1:dimPlex), & shape=[1,dimPlex*dimPlex],order=[2,1]),BMat))),MatA) - MatB = MatB + & - matmul(reshape(homogenization_F(1:dimPlex,1:dimPlex,qPt+1,cell+1),shape=[1,dimPlex*dimPlex]),MatA) + MatB = MatB & + + matmul(reshape(homogenization_F(1:dimPlex,1:dimPlex,m),shape=[1,dimPlex*dimPlex]),MatA) FAvg = FAvg + F BMatAvg = BMatAvg + BMat else