From 7e683ca7c3741f9f01acb05545d83d0f891e1b9e Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 14 Nov 2012 14:38:10 +0000 Subject: [PATCH] removed coordinates from call to CPFEM_general and introduced direct storage of coordinates to mesh_ipCoordinates --- code/CPFEM.f90 | 30 +-- code/DAMASK_abaqus_exp.f | 12 +- code/DAMASK_abaqus_std.f | 7 +- code/DAMASK_marc.f90 | 22 +- code/DAMASK_spectral.f90 | 22 +- code/DAMASK_spectral_solverAL.f90 | 23 +- code/DAMASK_spectral_solverBasic.f90 | 32 ++- code/DAMASK_spectral_solverBasicPETSc.f90 | 49 ++--- code/DAMASK_spectral_utilities.f90 | 7 +- code/mesh.f90 | 253 +++++++++++++++++----- 10 files changed, 295 insertions(+), 162 deletions(-) diff --git a/code/CPFEM.f90 b/code/CPFEM.f90 index 20668d740..3d5aedeb8 100644 --- a/code/CPFEM.f90 +++ b/code/CPFEM.f90 @@ -230,15 +230,14 @@ end subroutine CPFEM_init !*** perform initialization at first call, update variables and *** !*** call the actual material model *** !*********************************************************************** -subroutine CPFEM_general(mode, coords, ffn, ffn1, Temperature, dt, element, IP, cauchyStress,& +subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchyStress,& & jacobian, pstress, dPdF) ! note: cauchyStress = Cauchy stress cs(6) and jacobian = Consistent tangent dcs/dE !*** variables and functions from other modules ***! use prec, only: pInt use numerics, only: defgradTolerance, & - iJacoStiffness, & - numerics_unitlength + iJacoStiffness use debug, only: debug_level, & debug_CPFEM, & debug_levelBasic, & @@ -277,12 +276,6 @@ subroutine CPFEM_general(mode, coords, ffn, ffn1, Temperature, dt, element, IP, mesh_NcpElems, & mesh_maxNips, & mesh_element, & - mesh_node0, & - mesh_node, & - mesh_ipCoordinates, & - mesh_build_subNodeCoords, & - mesh_build_ipVolumes, & - mesh_build_ipCoordinates, & FE_Nips, & FE_Nnodes use material, only: homogenization_maxNgrains, & @@ -324,9 +317,6 @@ subroutine CPFEM_general(mode, coords, ffn, ffn1, Temperature, dt, element, IP, IP ! FE integration point number real(pReal), intent(inout) :: Temperature ! temperature real(pReal), intent(in) :: dt ! time increment - real(pReal), dimension (3,*), intent(in) :: coords ! MARC: displacements for each node of the current element - ! ABAQUS: coordinates of the current material point (IP) - ! SPECTRAL: coordinates of the current material point (IP) real(pReal), dimension (3,3), intent(in) :: ffn, & ! deformation gradient for t=t0 ffn1 ! deformation gradient for t=t1 integer(pInt), intent(in) :: mode ! computation mode 1: regular computation plus aging of results @@ -356,9 +346,6 @@ subroutine CPFEM_general(mode, coords, ffn, ffn1, Temperature, dt, element, IP, jacobian3333 ! jacobian in Matrix notation integer(pInt) cp_en, & ! crystal plasticity element number i, j, k, l, m, n, e -#ifdef Marc - integer(pInt):: node, FEnodeID -#endif logical updateJaco ! flag indicating if JAcobian has to be updated @@ -535,11 +522,6 @@ subroutine CPFEM_general(mode, coords, ffn, ffn1, Temperature, dt, element, IP, write(6,'(a,i8,a,i8)') '<< CPFEM >> Calculation for elements ',FEsolving_execElem(1),' to ',FEsolving_execElem(2) !$OMP END CRITICAL (write2out) endif -#ifdef Marc -! marc returns nodal coordinates, whereas Abaqus and spectral solver return ip coordinates. So for marc we have to calculate the ip coordinates from the nodal coordinates. - call mesh_build_subNodeCoords() ! update subnodal coordinates - call mesh_build_ipCoordinates() ! update ip coordinates -#endif if (iand(debug_level(debug_CPFEM), debug_levelExtensive) /= 0_pInt) then !$OMP CRITICAL (write2out) write(6,'(a,i8,a,i8)') '<< CPFEM >> Start stress and tangent ',FEsolving_execElem(1),' to ',FEsolving_execElem(2) @@ -611,14 +593,6 @@ subroutine CPFEM_general(mode, coords, ffn, ffn1, Temperature, dt, element, IP, CPFEM_cs(1:6,IP,cp_en) = rnd * CPFEM_odd_stress CPFEM_dcsde(1:6,1:6,IP,cp_en) = CPFEM_odd_jacobian * math_identity2nd(6) CPFEM_calc_done = .false. -#ifndef Marc - mesh_ipCoordinates(1:3,IP,cp_en) = numerics_unitlength * coords(1:3,1) -#else - do node = 1,FE_Nnodes(mesh_element(2,cp_en)) - FEnodeID = mesh_FEasCP('node',mesh_element(4+node,cp_en)) - mesh_node(1:3,FEnodeID) = mesh_node0(1:3,FEnodeID) + numerics_unitlength * coords(1:3,node) - enddo -#endif ! --+>> RECYCLING OF FORMER RESULTS (MARC SPECIALTY) <<+-- diff --git a/code/DAMASK_abaqus_exp.f b/code/DAMASK_abaqus_exp.f index 3d4149206..91ff6bdd3 100644 --- a/code/DAMASK_abaqus_exp.f +++ b/code/DAMASK_abaqus_exp.f @@ -170,6 +170,7 @@ subroutine vumat (jblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal, & use prec, only: pReal, & pInt + use numerics, only: numerics_unitlength use FEsolving, only: cycleCounter, & theTime, & outdatedByNewInc, & @@ -182,7 +183,8 @@ subroutine vumat (jblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal, & debug_levelBasic, & debug_level, & debug_abaqus - use mesh, only: mesh_FEasCP + use mesh, only: mesh_FEasCP, & + mesh_ipCoordinates use CPFEM, only: CPFEM_general,CPFEM_init_done, CPFEM_initAll use homogenization, only: materialpoint_sizeResults, materialpoint_results @@ -211,7 +213,7 @@ subroutine vumat (jblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal, & real(pReal), dimension(6) :: stress real(pReal), dimension(6,6) :: ddsdde real(pReal) temp, timeInc - integer(pInt) computationMode, n, i + integer(pInt) computationMode, n, i, cp_en logical :: cutBack do n = 1,nblock ! loop over vector of IPs @@ -289,8 +291,10 @@ subroutine vumat (jblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal, & defgrd1(3,2) = defgradNew(n,8) endif - coordinates = coordMp(n,1:3) - call CPFEM_general(computationMode,coordinates,defgrd0,defgrd1,temp,timeInc,nElement(n),nMatPoint(n),stress,ddsdde, pstress, dPdF) + cp_en = mesh_FEasCP('elem',nElement(n)) + mesh_ipCoordinates(1:3,n,cp_en) = numerics_unitlength * coordMp(n,1:3) + + call CPFEM_general(computationMode,defgrd0,defgrd1,temp,timeInc,cp_en,nMatPoint(n),stress,ddsdde, pstress, dPdF) ! Mandel: 11, 22, 33, SQRT(2)*12, SQRT(2)*23, SQRT(2)*13 ! straight: 11, 22, 33, 12, 23, 13 diff --git a/code/DAMASK_abaqus_std.f b/code/DAMASK_abaqus_std.f index 76d788ab7..c301175f6 100644 --- a/code/DAMASK_abaqus_std.f +++ b/code/DAMASK_abaqus_std.f @@ -130,6 +130,7 @@ subroutine UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,& use prec, only: pReal, & pInt + use numerics, only: numerics_unitlength use FEsolving, only: cycleCounter, & theInc, & calcMode, & @@ -147,7 +148,8 @@ subroutine UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,& debug_levelBasic, & debug_level, & debug_abaqus - use mesh, only: mesh_FEasCP + use mesh, only: mesh_FEasCP, & + mesh_ipCoordinates use CPFEM, only: CPFEM_general,CPFEM_init_done, CPFEM_initAll use homogenization, only: materialpoint_sizeResults, materialpoint_results @@ -253,6 +255,7 @@ subroutine UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,& else computationMode = 3 ! plain collect endif + mesh_ipCoordinates(1:3,npt,cp_en) = numerics_unitlength * COORDS endif theTime = time(2) ! record current starting time @@ -266,7 +269,7 @@ subroutine UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,& !$OMP END CRITICAL (write2out) endif - call CPFEM_general(computationMode,COORDS,dfgrd0,dfgrd1,temp,dtime,noel,npt,stress_h,ddsdde_h, pstress, dPdF) + call CPFEM_general(computationMode,dfgrd0,dfgrd1,temp,dtime,noel,npt,stress_h,ddsdde_h, pstress, dPdF) ! Mandel: 11, 22, 33, SQRT(2)*12, SQRT(2)*23, SQRT(2)*13 ! straight: 11, 22, 33, 12, 23, 13 diff --git a/code/DAMASK_marc.f90 b/code/DAMASK_marc.f90 index e6b6add79..59cb8f2f5 100644 --- a/code/DAMASK_marc.f90 +++ b/code/DAMASK_marc.f90 @@ -221,6 +221,7 @@ subroutine hypela2(& use prec, only: pReal, & pInt + use numerics, only: numerics_unitlength use FEsolving, only: cycleCounter, & theInc, & calcMode, & @@ -235,7 +236,13 @@ subroutine hypela2(& use math, only: invnrmMandel use debug, only: debug_info, & debug_reset - use mesh, only: mesh_FEasCP + use mesh, only: mesh_FEasCP, & + mesh_element, & + mesh_node0, & + mesh_node, & + mesh_build_subNodeCoords, & + mesh_build_ipCoordinates, & + FE_Nnodes use CPFEM, only: CPFEM_initAll,CPFEM_general,CPFEM_init_done !$ use numerics, only: DAMASK_NumThreadsInt ! number of threads set by DAMASK_NUM_THREADS @@ -265,7 +272,8 @@ subroutine hypela2(& real(pReal), dimension (3,3) :: pstress ! dummy argument for call of cpfem_general (used by mpie_spectral) real(pReal), dimension (3,3,3,3) :: dPdF ! dummy argument for call of cpfem_general (used by mpie_spectral) - integer(pInt) computationMode, i, cp_en + integer(pInt) computationMode, i, cp_en + integer(pInt) node, FEnodeID ! OpenMP variable !$ integer(pInt) defaultNumThreadsInt ! default value set by Marc @@ -351,6 +359,10 @@ subroutine hypela2(& else computationMode = 3 ! plain collect endif + do node = 1,FE_Nnodes(mesh_element(2,cp_en)) + FEnodeID = mesh_FEasCP('node',mesh_element(4+node,cp_en)) + mesh_node(1:3,FEnodeID) = mesh_node0(1:3,FEnodeID) + numerics_unitlength * dispt(1:3,node) + enddo endif theTime = cptim ! record current starting time @@ -359,7 +371,11 @@ subroutine hypela2(& lastMode = calcMode(nn,cp_en) ! record calculationMode endif - call CPFEM_general(computationMode,dispt,ffn,ffn1,t(1),timinc,n(1),nn,stress,ddsdde, pstress, dPdF) + ! marc returns nodal coordinates. So for marc we have to calculate the ip coordinates from the nodal coordinates. + call mesh_build_subNodeCoords() ! update subnodal coordinates + call mesh_build_ipCoordinates() ! update ip coordinates + + call CPFEM_general(computationMode,ffn,ffn1,t(1),timinc,n(1),nn,stress,ddsdde, pstress, dPdF) ! Mandel: 11, 22, 33, SQRT(2)*12, SQRT(2)*23, SQRT(2)*13 ! Marc: 11, 22, 33, 12, 23, 13 diff --git a/code/DAMASK_spectral.f90 b/code/DAMASK_spectral.f90 index cd5346251..5420219f9 100644 --- a/code/DAMASK_spectral.f90 +++ b/code/DAMASK_spectral.f90 @@ -84,7 +84,8 @@ program DAMASK_spectral wgt, & geomdim, & virt_dim, & - deformed_FFT + deformed_FFT, & + mesh_ipCoordinates use CPFEM, only: & CPFEM_general, & @@ -519,15 +520,20 @@ program DAMASK_spectral if (debugRestart) write(6,'(a)') 'Data read in' endif ielem = 0_pInt + do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) + ielem = ielem +1_pInt + mesh_ipCoordinates(1:3,1,ielem) = coordinates(i,j,k,1:3) + enddo; enddo; enddo + ielem = 0_pInt do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) ielem = ielem + 1_pInt - call CPFEM_general(3_pInt,coordinates(i,j,k,1:3),F(i,j,k,1:3,1:3),F(i,j,k,1:3,1:3),temperature(i,j,k),& + call CPFEM_general(3_pInt,F(i,j,k,1:3,1:3),F(i,j,k,1:3,1:3),temperature(i,j,k),& 0.0_pReal,ielem,1_pInt,sigma,dsde,P_real(i,j,k,1:3,1:3),dPdF) enddo; enddo; enddo ielem = 0_pInt do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) ielem = ielem + 1_pInt - call CPFEM_general(2_pInt,coordinates(i,j,k,1:3),F(i,j,k,1:3,1:3),F(i,j,k,1:3,1:3),temperature(i,j,k),& + call CPFEM_general(2_pInt,F(i,j,k,1:3,1:3),F(i,j,k,1:3,1:3),temperature(i,j,k),& 0.0_pReal,ielem,1_pInt,sigma,dsde,P_real(i,j,k,1:3,1:3),dPdF) C = C + dPdF enddo; enddo; enddo @@ -686,7 +692,11 @@ program DAMASK_spectral enddo; enddo; enddo call deformed_fft(res,geomdim,math_rotate_backward33(F_aim,bc(loadcase)%rotation),& ! calculate current coordinates 1.0_pReal,F_lastInc,coordinates) - + ielem = 0_pInt + do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) + ielem = ielem +1_pInt + mesh_ipCoordinates(1:3,1,ielem) = coordinates(i,j,k,1:3) + enddo; enddo; enddo !-------------------------------------------------------------------------------------------------- ! calculate reduced compliance if(size_reduced > 0_pInt) then ! calculate compliance in case stress BC is applied @@ -753,7 +763,7 @@ program DAMASK_spectral do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) ielem = ielem + 1_pInt call CPFEM_general(3_pInt,& ! collect cycle - coordinates(i,j,k,1:3), F_lastInc(i,j,k,1:3,1:3),F(i,j,k,1:3,1:3), & + F_lastInc(i,j,k,1:3,1:3),F(i,j,k,1:3,1:3), & temperature(i,j,k),timeinc,ielem,1_pInt,sigma,dsde,& P_real(i,j,k,1:3,1:3),dPdF) enddo; enddo; enddo @@ -765,7 +775,7 @@ program DAMASK_spectral do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) ielem = ielem + 1_pInt call CPFEM_general(CPFEM_mode,& ! first element in first iteration retains CPFEM_mode 1, - coordinates(i,j,k,1:3),F_lastInc(i,j,k,1:3,1:3), F(i,j,k,1:3,1:3), & ! others get 2 (saves winding forward effort) + F_lastInc(i,j,k,1:3,1:3), F(i,j,k,1:3,1:3), & ! others get 2 (saves winding forward effort) temperature(i,j,k),timeinc,ielem,1_pInt,sigma,dsde, & P_real(i,j,k,1:3,1:3),dPdF) CPFEM_mode = 2_pInt diff --git a/code/DAMASK_spectral_solverAL.f90 b/code/DAMASK_spectral_solverAL.f90 index 43d9d51ae..a7f0eae2a 100644 --- a/code/DAMASK_spectral_solverAL.f90 +++ b/code/DAMASK_spectral_solverAL.f90 @@ -46,7 +46,6 @@ module DAMASK_spectral_solverAL F_lambda_lastInc, & !< field of previous incompatible deformation gradient Fdot, & !< field of assumed rate of compatible deformation gradient F_lambdaDot !< field of assumed rate of incopatible deformation gradient - real(pReal), private, dimension(:,:,:,:), allocatable :: coordinates real(pReal), private :: temperature !< temperature, no spatial quantity at the moment !-------------------------------------------------------------------------------------------------- @@ -94,7 +93,8 @@ subroutine AL_init() use mesh, only: & res, & geomdim, & - mesh_NcpElems + mesh_NcpElems, & + mesh_ipCoordinates use math, only: & math_invSym3333 @@ -119,7 +119,6 @@ subroutine AL_init() ! allocate (Fdot,source = F_lastInc) somethin like that should be possible allocate (F_lambda_lastInc(3,3, res(1), res(2),res(3)), source = 0.0_pReal) allocate (F_lambdaDot(3,3, res(1), res(2),res(3)), source = 0.0_pReal) - allocate (coordinates( res(1), res(2),res(3),3), source = 0.0_pReal) !-------------------------------------------------------------------------------------------------- ! PETSc Init @@ -151,10 +150,6 @@ subroutine AL_init() F_lambda_lastInc = F_lastInc F = reshape(F_lastInc,[9,res(1),res(2),res(3)]) F_lambda = F - do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) - coordinates(i,j,k,1:3) = geomdim/real(res,pReal)*real([i,j,k],pReal) & - - geomdim/real(2_pInt*res,pReal) - enddo; enddo; enddo elseif (restartInc > 1_pInt) then ! using old values from file if (debugRestart) write(6,'(a,i6,a)') 'Reading values of increment ',& restartInc - 1_pInt,' from file' @@ -180,9 +175,10 @@ subroutine AL_init() call IO_read_jobBinaryFile(777,'F_aim_lastInc',trim(getSolverJobName()),size(F_aim_lastInc)) read (777,rec=1) F_aim_lastInc close (777) - coordinates = 0.0 ! change it later!!! endif - call Utilities_constitutiveResponse(coordinates,F,F,temperature,0.0_pReal,P,C,P_av,.false.,math_I3) + mesh_ipCoordinates = 0.0_pReal !reshape(mesh_deformedCoordsFFT(geomdim,& + !reshape(F,[3,3,res(1),res(2),res(3)])),[3,1,mesh_NcpElems]) + call Utilities_constitutiveResponse(F,F,temperature,0.0_pReal,P,C,P_av,.false.,math_I3) call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr) CHKERRQ(ierr) @@ -219,7 +215,7 @@ type(tSolutionState) function & use mesh, only: & res,& geomdim,& - deformed_fft + mesh_ipCoordinates use IO, only: & IO_write_JobBinaryFile use DAMASK_spectral_Utilities, only: & @@ -292,7 +288,8 @@ type(tSolutionState) function & !-------------------------------------------------------------------------------------------------- ! update coordinates and rate and forward last inc - + mesh_ipCoordinates = 0.0_pReal !reshape(mesh_deformedCoordsFFT(geomdim,& + !reshape(F,[3,3,res(1),res(2),res(3)])),[3,1,mesh_NcpElems]) Fdot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), & timeinc,timeinc_old,guess,F_lastInc,reshape(F,[3,3,res(1),res(2),res(3)])) F_lambdaDot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), & @@ -300,8 +297,6 @@ type(tSolutionState) function & F_lastInc = reshape(F,[3,3,res(1),res(2),res(3)]) F_lambda_lastInc = reshape(F_lambda,[3,3,res(1),res(2),res(3)]) - call deformed_fft(res,geomdim,math_rotate_backward33(F_aim_lastInc,rotation_BC), & - 1.0_pReal,F_lastInc,coordinates) endif F_aim = F_aim + f_aimDot * timeinc @@ -406,7 +401,7 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr) !-------------------------------------------------------------------------------------------------- ! evaluate constitutive response - call Utilities_constitutiveResponse(coordinates,F_lastInc,F,temperature,params%timeinc, & + call Utilities_constitutiveResponse(F_lastInc,F,temperature,params%timeinc, & residual_F,C,P_av,ForwardData,params%rotation_BC) ForwardData = .False. diff --git a/code/DAMASK_spectral_solverBasic.f90 b/code/DAMASK_spectral_solverBasic.f90 index b5d9b3409..f0330fbda 100644 --- a/code/DAMASK_spectral_solverBasic.f90 +++ b/code/DAMASK_spectral_solverBasic.f90 @@ -27,7 +27,6 @@ module DAMASK_spectral_SolverBasic F, & !< deformation gradient field F_lastInc, & !< deformation gradient field last increment Fdot !< assumed rate for F n to F n+1 - real(pReal), private, dimension(:,:,:,:), allocatable :: coordinates real(pReal), private :: temperature !< not pointwise !-------------------------------------------------------------------------------------------------- @@ -63,7 +62,10 @@ subroutine basic_init() debugRestart use mesh, only: & res, & - geomdim + geomdim, & + mesh_ipCoordinates, & + mesh_NcpElems, & + mesh_deformedCoordsFFT implicit none real(pReal), dimension(3,3,res(1),res(2),res(3)) :: P @@ -86,18 +88,12 @@ subroutine basic_init() allocate (F ( 3,3,res(1), res(2),res(3)), source = 0.0_pReal) allocate (F_lastInc ( 3,3,res(1), res(2),res(3)), source = 0.0_pReal) allocate (Fdot ( 3,3,res(1), res(2),res(3)), source = 0.0_pReal) - allocate (coordinates( res(1), res(2),res(3),3),source = 0.0_pReal) !-------------------------------------------------------------------------------------------------- ! init fields and average quantities if (restartInc == 1_pInt) then ! no deformation (no restart) F = spread(spread(spread(math_I3,3,res(1)),4,res(2)),5,res(3)) ! initialize to identity F_lastInc = F - do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) - coordinates(i,j,k,1:3) = geomdim/real(res,pReal)*real([i,j,k],pReal) & - - geomdim/real(2_pInt*res,pReal) - enddo; enddo; enddo - elseif (restartInc > 1_pInt) then ! using old values from file if (debugRestart) write(6,'(a,'//IO_intOut(restartInc-1_pInt)//',a)') & 'Reading values of increment', restartInc - 1_pInt, 'from file' @@ -127,11 +123,10 @@ subroutine basic_init() call IO_read_jobBinaryFile(777,'C_ref',trim(getSolverJobName()),size(temp3333_Real)) read (777,rec=1) temp3333_Real close (777) - coordinates = 0.0 ! change it later!!! endif - call Utilities_constitutiveResponse(coordinates,F,F,temperature,0.0_pReal,P,C,temp33_Real,.false.,math_I3) ! constitutive response with no deformation in no time to get reference stiffness - !no rotation bc call deformed_fft(res,geomdim,math_rotate_backward33(F_aim,rotation_BC),1.0_pReal,F_lastInc,coordinates) - + mesh_ipCoordinates = 0.0_pReal !reshape(mesh_deformedCoordsFFT(geomdim,& + !reshape(F,[3,3,res(1),res(2),res(3)])),[3,1,mesh_NcpElems]) + call Utilities_constitutiveResponse(F,F,temperature,0.0_pReal,P,C,temp33_Real,.false.,math_I3) ! constitutive response with no deformation in no time to get reference stiffness if (restartInc == 1_pInt) then ! use initial stiffness as reference stiffness temp3333_Real = C endif @@ -158,8 +153,10 @@ type(tSolutionState) function & use mesh, only: & res,& geomdim, & - deformed_fft, & - wgt + wgt, & + mesh_ipCoordinates,& + mesh_NcpElems, & + mesh_deformedCoordsFFT use IO, only: & IO_write_JobBinaryFile, & IO_intOut @@ -252,6 +249,8 @@ type(tSolutionState) function & C = C_lastInc else C_lastInc = C + mesh_ipCoordinates = 0.0_pReal !reshape(mesh_deformedCoordsFFT(geomdim,& + !reshape(F,[3,3,res(1),res(2),res(3)])),[3,1,mesh_NcpElems]) if (F_BC%myType=='l') then ! calculate f_aimDot from given L and current F f_aimDot = F_BC%maskFloat * math_mul33x33(F_BC%values, F_aim) elseif(F_BC%myType=='fdot') then ! f_aimDot is prescribed @@ -259,9 +258,6 @@ type(tSolutionState) function & endif if (guess) f_aimDot = f_aimDot + P_BC%maskFloat * (F_aim - F_aim_lastInc)/timeinc_old F_aim_lastInc = F_aim - - call deformed_fft(res,geomdim,math_rotate_backward33(F_aim_lastInc,rotation_BC), & ! update coordinates and rate and forward last inc - 1.0_pReal,F_lastInc,coordinates) Fdot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), & timeinc,timeinc_old,guess,F_lastInc,F) F_lastInc = F @@ -291,7 +287,7 @@ type(tSolutionState) function & ! evaluate constitutive response F_aim_lab_lastIter = math_rotate_backward33(F_aim,rotation_BC) basic_solution%termIll = .false. - call Utilities_constitutiveResponse(coordinates,F_lastInc,F,temperature,timeinc,& + call Utilities_constitutiveResponse(F_lastInc,F,temperature,timeinc,& P,C,P_av,ForwardData,rotation_BC) basic_solution%termIll = terminallyIll terminallyIll = .false. diff --git a/code/DAMASK_spectral_solverBasicPETSc.f90 b/code/DAMASK_spectral_solverBasicPETSc.f90 index d54173f7b..4e28f6425 100644 --- a/code/DAMASK_spectral_solverBasicPETSc.f90 +++ b/code/DAMASK_spectral_solverBasicPETSc.f90 @@ -10,10 +10,8 @@ module DAMASK_spectral_SolverBasicPETSc use prec, only: & pInt, & pReal - use math, only: & math_I3 - use DAMASK_spectral_Utilities, only: & tSolutionState @@ -44,7 +42,6 @@ module DAMASK_spectral_SolverBasicPETSc !-------------------------------------------------------------------------------------------------- ! common pointwise data real(pReal), private, dimension(:,:,:,:,:), allocatable :: F_lastInc, Fdot - real(pReal), private, dimension(:,:,:,:), allocatable :: coordinates real(pReal) :: temperature !-------------------------------------------------------------------------------------------------- @@ -72,33 +69,27 @@ module DAMASK_spectral_SolverBasicPETSc !> @brief allocates all neccessary fields and fills them with data, potentially from restart info !-------------------------------------------------------------------------------------------------- subroutine basicPETSc_init() - use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran >4.6 at the moment) - use IO, only: & IO_read_JobBinaryFile, & IO_write_JobBinaryFile - use FEsolving, only: & restartInc - use DAMASK_interface, only: & getSolverJobName - use DAMASK_spectral_Utilities, only: & Utilities_init, & Utilities_constitutiveResponse, & Utilities_updateGamma, & debugRestart - use numerics, only: & petsc_options - use mesh, only: & res, & geomdim, & - mesh_NcpElems - + mesh_NcpElems, & + mesh_ipCoordinates, & + mesh_deformedCoordsFFT use math, only: & math_invSym3333 @@ -106,7 +97,7 @@ subroutine basicPETSc_init() #include #include integer(pInt) :: i,j,k - real(pReal), dimension(3,3, res(1), res(2),res(3)) :: P + real(pReal), dimension(3,3,res(1),res(2),res(3)) :: P PetscScalar, dimension(:,:,:,:), pointer :: F PetscErrorCode :: ierr PetscObject :: dummy @@ -119,9 +110,8 @@ subroutine basicPETSc_init() !-------------------------------------------------------------------------------------------------- ! allocate global fields - allocate (F_lastInc (3,3, res(1), res(2),res(3)), source = 0.0_pReal) - allocate (Fdot (3,3, res(1), res(2),res(3)), source = 0.0_pReal) - allocate (coordinates( res(1), res(2),res(3),3), source = 0.0_pReal) + allocate (F_lastInc(3,3,res(1),res(2),res(3)), source = 0.0_pReal) + allocate (Fdot (3,3,res(1),res(2),res(3)), source = 0.0_pReal) !-------------------------------------------------------------------------------------------------- ! initialize solver specific parts of PETSc @@ -150,10 +140,6 @@ subroutine basicPETSc_init() if (restartInc == 1_pInt) then ! no deformation (no restart) F_lastInc = spread(spread(spread(math_I3,3,res(1)),4,res(2)),5,res(3)) ! initialize to identity F = reshape(F_lastInc,[9,res(1),res(2),res(3)]) - do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) - coordinates(i,j,k,1:3) = geomdim/real(res,pReal)*real([i,j,k],pReal) & - - geomdim/real(2_pInt*res,pReal) - enddo; enddo; enddo elseif (restartInc > 1_pInt) then ! using old values from file if (debugRestart) write(6,'(a,i6,a)') 'Reading values of increment ',& restartInc - 1_pInt,' from file' @@ -172,11 +158,10 @@ subroutine basicPETSc_init() call IO_read_jobBinaryFile(777,'F_aim_lastInc',trim(getSolverJobName()),size(F_aim_lastInc)) read (777,rec=1) F_aim_lastInc close (777) - - coordinates = 0.0 ! change it later!!! endif - - call Utilities_constitutiveResponse(coordinates,& + mesh_ipCoordinates = 0.0_pReal !reshape(mesh_deformedCoordsFFT(geomdim,& + !reshape(F,[3,3,res(1),res(2),res(3)])),[3,1,mesh_NcpElems]) + call Utilities_constitutiveResponse(& reshape(F(0:8,0:res(1)-1_pInt,0:res(2)-1_pInt,0:res(3)-1_pInt),[3,3,res(1),res(2),res(3)]),& reshape(F(0:8,0:res(1)-1_pInt,0:res(2)-1_pInt,0:res(3)-1_pInt),[3,3,res(1),res(2),res(3)]),& temperature,0.0_pReal,P,C,P_av,.false.,math_I3) @@ -211,7 +196,9 @@ type(tSolutionState) function & use mesh, only: & res,& geomdim,& - deformed_fft + mesh_ipCoordinates,& + mesh_NcpElems, & + mesh_deformedCoordsFFT use IO, only: & IO_write_JobBinaryFile use DAMASK_spectral_Utilities, only: & @@ -260,14 +247,16 @@ type(tSolutionState) function & close(777) endif call DMDAVecGetArrayF90(da,solution_vec,F,ierr) - + mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomdim,reshape(F,[3,3,res(1),res(2),res(3)])),& + [3,1,mesh_NcpElems]) if ( cutBack) then F_aim = F_aim_lastInc F = reshape(F_lastInc,[9,res(1),res(2),res(3)]) C = C_lastInc else C_lastInc = C - + mesh_ipCoordinates = 0.0_pReal !reshape(mesh_deformedCoordsFFT(geomdim,& + !reshape(F,[3,3,res(1),res(2),res(3)])),[3,1,mesh_NcpElems]) !-------------------------------------------------------------------------------------------------- ! calculate rate for aim if (F_BC%myType=='l') then ! calculate f_aimDot from given L and current F @@ -280,8 +269,6 @@ type(tSolutionState) function & !-------------------------------------------------------------------------------------------------- ! update coordinates and rate and forward last inc - call deformed_fft(res,geomdim,math_rotate_backward33(F_aim_lastInc,rotation_BC), & - 1.0_pReal,F_lastInc,coordinates) Fdot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), & timeinc,timeinc_old,guess,F_lastInc,reshape(F,[3,3,res(1),res(2),res(3)])) F_lastInc = reshape(F,[3,3,res(1),res(2),res(3)]) @@ -292,7 +279,6 @@ type(tSolutionState) function & F = reshape(Utilities_forwardField(timeinc,F_aim,F_lastInc,Fdot),[9,res(1),res(2),res(3)]) call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr) CHKERRQ(ierr) - call deformed_fft(res,geomdim,math_rotate_backward33(F_aim,rotation_BC),1.0_pReal,F_lastInc,coordinates) !-------------------------------------------------------------------------------------------------- ! update stiffness (and gamma operator) @@ -322,7 +308,6 @@ end function BasicPETSc_solution !> @brief forms the AL residual vector !-------------------------------------------------------------------------------------------------- subroutine BasicPETSC_formResidual(myIn,x_scal,f_scal,dummy,ierr) - use numerics, only: & itmax, & itmin @@ -366,7 +351,7 @@ subroutine BasicPETSC_formResidual(myIn,x_scal,f_scal,dummy,ierr) !-------------------------------------------------------------------------------------------------- ! evaluate constitutive response - call Utilities_constitutiveResponse(coordinates,F_lastInc,x_scal,temperature,params%timeinc, & + call Utilities_constitutiveResponse(F_lastInc,x_scal,temperature,params%timeinc, & f_scal,C,P_av,ForwardData,params%rotation_BC) ForwardData = .false. diff --git a/code/DAMASK_spectral_utilities.f90 b/code/DAMASK_spectral_utilities.f90 index 00f736331..df2d6a537 100644 --- a/code/DAMASK_spectral_utilities.f90 +++ b/code/DAMASK_spectral_utilities.f90 @@ -630,7 +630,7 @@ end function utilities_maskedCompliance !-------------------------------------------------------------------------------------------------- !> @brief calculates constitutive response !-------------------------------------------------------------------------------------------------- -subroutine utilities_constitutiveResponse(coordinates,F_lastInc,F,temperature,timeinc,& +subroutine utilities_constitutiveResponse(F_lastInc,F,temperature,timeinc,& P,C,P_av,forwardData,rotation_BC) use debug, only: & debug_reset, & @@ -648,7 +648,6 @@ subroutine utilities_constitutiveResponse(coordinates,F_lastInc,F,temperature,ti implicit none real(pReal), intent(inout) :: temperature !< temperature (no field) - real(pReal), intent(in), dimension(res(1),res(2),res(3),3) :: coordinates !< coordinates field real(pReal), intent(in), dimension(3,3,res(1),res(2),res(3)) :: & F_lastInc, & !< target deformation gradient F !< previous deformation gradient @@ -702,7 +701,7 @@ subroutine utilities_constitutiveResponse(coordinates,F_lastInc,F,temperature,ti do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) ielem = ielem + 1_pInt call CPFEM_general(collectMode,& ! collect cycle - coordinates(i,j,k,1:3), F_lastInc(1:3,1:3,i,j,k),F(1:3,1:3,i,j,k), & + F_lastInc(1:3,1:3,i,j,k),F(1:3,1:3,i,j,k), & temperature,timeinc,ielem,1_pInt,sigma,dsde,P(1:3,1:3,i,j,k),dPdF) collectMode = 3_pInt enddo; enddo; enddo @@ -714,7 +713,7 @@ subroutine utilities_constitutiveResponse(coordinates,F_lastInc,F,temperature,ti do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) ielem = ielem + 1_pInt call CPFEM_general(calcMode,& ! first element in first iteration retains CPFEM_mode 1, - coordinates(i,j,k,1:3),F_lastInc(1:3,1:3,i,j,k), F(1:3,1:3,i,j,k), & ! others get 2 (saves winding forward effort) + F_lastInc(1:3,1:3,i,j,k), F(1:3,1:3,i,j,k), & ! others get 2 (saves winding forward effort) temperature,timeinc,ielem,1_pInt,sigma,dsde,P(1:3,1:3,i,j,k),dPdF) calcMode = 2_pInt C = C + dPdF diff --git a/code/mesh.f90 b/code/mesh.f90 index d53f82a81..d70b62be1 100644 --- a/code/mesh.f90 +++ b/code/mesh.f90 @@ -1,7 +1,7 @@ -! Copyright 2011 Max-Planck-Institut für Eisenforschung GmbH +! Copyright 2011 Max-Planck-Institut für Eisenforschung GmbH ! ! This file is part of DAMASK, -! the Düsseldorf Advanced MAterial Simulation Kit. +! the Düsseldorf Advanced Material Simulation Kit. ! ! DAMASK is free software: you can redistribute it and/or modify ! it under the terms of the GNU General Public License as published by @@ -19,11 +19,11 @@ !-------------------------------------------------------------------------------------------------- !* $Id$ !-------------------------------------------------------------------------------------------------- -!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH -!! Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH -!! Christoph Koords, Max-Planck-Institut für Eisenforschung GmbH -!! Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH -!! Krishna Komerla, Max-Planck-Institut für Eisenforschung GmbH +!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH +!! Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH +!! Christoph Koords, Max-Planck-Institut für Eisenforschung GmbH +!! Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH +!! Krishna Komerla, Max-Planck-Institut für Eisenforschung GmbH !> @brief Sets up the mesh for the solvers MSC.Marc, Abaqus and the spectral solver !-------------------------------------------------------------------------------------------------- @@ -286,6 +286,7 @@ module mesh mesh_regular_grid, & deformed_linear, & deformed_fft, & + mesh_deformedCoordsFFT, & volume_compare, & shape_compare #endif @@ -1184,7 +1185,7 @@ function mesh_regrid(adaptive,resNewInput,minRes) integer(pInt), dimension(:,:,:), allocatable :: & material_phase, material_phaseNew, & sizeStateConst - + write(6,*) 'Regridding geometry' if (adaptive) then write(6,*) 'adaptive resolution determination' @@ -1355,7 +1356,7 @@ function mesh_regrid(adaptive,resNewInput,minRes) allocate(F_Linear(3,3,mesh_NcpElems)) allocate(F_Linear_New(3,3,NpointsNew)) allocate(FNew(resNew(1),resNew(2),resNew(3),3,3)) - + ielem = 0_pInt do k=1_pInt,res(3); do j=1_pInt, res(2); do i=1_pInt, res(1) ielem = ielem + 1_pInt @@ -1365,7 +1366,7 @@ function mesh_regrid(adaptive,resNewInput,minRes) do i=1_pInt, NpointsNew F_Linear_New(1:3,1:3,i) = F_Linear(1:3,1:3,indices(i)) ! -- mapping old to new ...based on indices enddo - + ielem = 0_pInt do k=1_pInt,resNew(3); do j=1_pInt, resNew(2); do i=1_pInt, resNew(1) ielem = ielem + 1_pInt @@ -1377,11 +1378,11 @@ function mesh_regrid(adaptive,resNewInput,minRes) enddo; enddo deltaF = Favg - FavgNew - + do k=1_pInt,resNew(3); do j=1_pInt, resNew(2); do i=1_pInt, resNew(1) FNew(i,j,k,1:3,1:3) = FNew(i,j,k,1:3,1:3) + deltaF enddo; enddo; enddo - + call IO_write_jobBinaryFile(777,'convergedSpectralDefgrad',size(FNew)) write (777,rec=1) FNew close (777) @@ -1394,12 +1395,12 @@ function mesh_regrid(adaptive,resNewInput,minRes) allocate(F_lastIncNew(resNew(1),resNew(2),resNew(3),3,3)) allocate(F_Linear(3,3,mesh_NcpElems)) allocate(F_Linear_New(3,3,NpointsNew)) - + call IO_read_jobBinaryFile(777,'convergedSpectralDefgrad_lastInc', & trim(getSolverJobName()),size(F_lastInc)) read (777,rec=1) F_lastInc close (777) - + call IO_read_jobBinaryFile(777,'F_aim_lastInc', & trim(getSolverJobName()),size(Favg_LastInc)) read (777,rec=1) Favg_LastInc @@ -1433,7 +1434,7 @@ function mesh_regrid(adaptive,resNewInput,minRes) do k=1_pInt,resNew(3); do j=1_pInt, resNew(2); do i=1_pInt, resNew(1) F_LastIncNew(i,j,k,1:3,1:3) = F_LastIncNew(i,j,k,1:3,1:3) + deltaF_lastInc enddo; enddo; enddo - + call IO_write_jobBinaryFile(777,'convergedSpectralDefgrad_lastInc',size(F_LastIncNew)) write (777,rec=1) F_LastIncNew close (777) @@ -1441,7 +1442,7 @@ function mesh_regrid(adaptive,resNewInput,minRes) deallocate(F_Linear_New) deallocate(F_lastInc) deallocate(F_lastIncNew) - + ! relocating data of material subroutine --------------------------------------------------------- allocate(material_phase (1,1, mesh_NcpElems)) allocate(material_phaseNew (1,1, NpointsNew)) @@ -1924,6 +1925,156 @@ subroutine deformed_fft(res,geomdim,defgrad_av,scaling,defgrad,coords) end subroutine deformed_fft +!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +function mesh_deformedCoordsFFT(geomdim,F,scalingIn,FavgIn) +!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +! Routine to calculate coordinates in current configuration for given defgrad +! using integration in Fourier space (more accurate than deformed(...)) +! + use IO, only: & + IO_error + use numerics, only: & + fftw_timelimit, & + fftw_planner_flag + use debug, only: & + debug_mesh, & + debug_level, & + debug_levelBasic + use math, only: & + PI + + implicit none + + real(pReal), intent(in), dimension(3) :: geomdim + real(pReal), intent(in), dimension(:,:,:,:,:) :: F + real(pReal), intent(in), dimension(3,3), optional :: FavgIn + real(pReal), intent(in), optional :: scalingIn +! function + real(pReal), dimension(3,size(F,3),size(F,4),size(F,5)) :: mesh_deformedCoordsFFT +! allocatable arrays for fftw c routines + type(C_PTR) :: fftw_forth, fftw_back + type(C_PTR) :: coords_fftw, defgrad_fftw + real(pReal), dimension(:,:,:,:,:), pointer :: F_real + complex(pReal), dimension(:,:,:,:,:), pointer :: F_fourier + real(pReal), dimension(:,:,:,:), pointer :: coords_real + complex(pReal), dimension(:,:,:,:), pointer :: coords_fourier + ! other variables + integer(pInt) :: i, j, k, m, res1_red + integer(pInt), dimension(3) :: k_s, res + real(pReal), dimension(3) :: step, offset_coords, integrator + real(pReal), dimension(3,3) :: Favg + real(pReal) :: scaling + + if (present(scalingIn)) then + if (scalingIn < 0.0_pReal) then !the f2py way to tell it is not present + scaling = 1.0_pReal + else + scaling = scalingIn + endif + else + scaling = 1.0_pReal + endif + + res = [size(F,3),size(F,4),size(F,5)] + integrator = geomdim / 2.0_pReal / pi ! see notes where it is used + + if (iand(debug_level(debug_mesh),debug_levelBasic) /= 0_pInt) then + print*, 'Restore geometry using FFT-based integration' + print '(a,3(e12.5))', ' Dimension: ', geomdim + print '(a,3(i5))', ' Resolution:', res + endif + + res1_red = res(1)/2_pInt + 1_pInt ! size of complex array in first dimension (c2r, r2c) + step = geomdim/real(res, pReal) + + if (pReal /= C_DOUBLE .or. pInt /= C_INT) call IO_error(error_ID=808_pInt) + call fftw_set_timelimit(fftw_timelimit) + defgrad_fftw = fftw_alloc_complex(int(res1_red *res(2)*res(3)*9_pInt,C_SIZE_T)) !C_SIZE_T is of type integer(8) + call c_f_pointer(defgrad_fftw, F_real, [res(1)+2_pInt,res(2),res(3),3_pInt,3_pInt]) + call c_f_pointer(defgrad_fftw, F_fourier,[res1_red ,res(2),res(3),3_pInt,3_pInt]) + coords_fftw = fftw_alloc_complex(int(res1_red *res(2)*res(3)*3_pInt,C_SIZE_T)) !C_SIZE_T is of type integer(8) + call c_f_pointer(coords_fftw, coords_real, [res(1)+2_pInt,res(2),res(3),3_pInt]) + call c_f_pointer(coords_fftw, coords_fourier, [res1_red ,res(2),res(3),3_pInt]) + fftw_forth = fftw_plan_many_dft_r2c(3_pInt,(/res(3),res(2) ,res(1)/),9_pInt,& ! dimensions , length in each dimension in reversed order + F_real,(/res(3),res(2) ,res(1)+2_pInt/),& ! input data , physical length in each dimension in reversed order + 1_pInt, res(3)*res(2)*(res(1)+2_pInt),& ! striding , product of physical lenght in the 3 dimensions + F_fourier,(/res(3),res(2) ,res1_red/),& + 1_pInt, res(3)*res(2)* res1_red,fftw_planner_flag) + + fftw_back = fftw_plan_many_dft_c2r(3_pInt,(/res(3),res(2) ,res(1)/),3_pInt,& + coords_fourier,(/res(3),res(2) ,res1_red/),& + 1_pInt, res(3)*res(2)* res1_red,& + coords_real,(/res(3),res(2) ,res(1)+2_pInt/),& + 1_pInt, res(3)*res(2)*(res(1)+2_pInt),fftw_planner_flag) + + + do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) + F_real(i,j,k,1:3,1:3) = F(1:3,1:3,i,j,k) ! ensure that data is aligned properly (fftw_alloc) + enddo; enddo; enddo + + call fftw_execute_dft_r2c(fftw_forth, F_real, F_fourier) + + if (present(FavgIn)) then + if (all(FavgIn < 0.0_pReal)) then + Favg = real(F_fourier(1,1,1,1:3,1:3)*real((res(1)*res(2)*res(3)),pReal),pReal) !the f2py way to tell it is not present + else + Favg = FavgIn + endif + else + Favg = real(F_fourier(1,1,1,1:3,1:3)*real((res(1)*res(2)*res(3)),pReal),pReal) + endif + + !remove highest frequency in each direction + if(res(1)>1_pInt) & + F_fourier( res(1)/2_pInt+1_pInt,1:res(2) ,1:res(3) ,& + 1:3,1:3) = cmplx(0.0_pReal,0.0_pReal,pReal) + if(res(2)>1_pInt) & + F_fourier(1:res1_red ,res(2)/2_pInt+1_pInt,1:res(3) ,& + 1:3,1:3) = cmplx(0.0_pReal,0.0_pReal,pReal) + if(res(3)>1_pInt) & + F_fourier(1:res1_red ,1:res(2) ,res(3)/2_pInt+1_pInt,& + 1:3,1:3) = cmplx(0.0_pReal,0.0_pReal,pReal) + + coords_fourier = cmplx(0.0_pReal,0.0_pReal,pReal) + do k = 1_pInt, res(3) + k_s(3) = k-1_pInt + if(k > res(3)/2_pInt+1_pInt) k_s(3) = k_s(3)-res(3) + do j = 1_pInt, res(2) + k_s(2) = j-1_pInt + if(j > res(2)/2_pInt+1_pInt) k_s(2) = k_s(2)-res(2) + do i = 1_pInt, res1_red + k_s(1) = i-1_pInt + do m = 1_pInt,3_pInt + coords_fourier(i,j,k,m) = sum(F_fourier(i,j,k,m,1:3)*cmplx(0.0_pReal,real(k_s,pReal)*integrator,pReal)) + enddo + if (k_s(3) /= 0_pInt .or. k_s(2) /= 0_pInt .or. k_s(1) /= 0_pInt) & + coords_fourier(i,j,k,1:3) = coords_fourier(i,j,k,1:3) / real(-sum(k_s*k_s),pReal) + enddo; enddo; enddo + + call fftw_execute_dft_c2r(fftw_back,coords_fourier,coords_real) + coords_real = coords_real/real(res(1)*res(2)*res(3),pReal) + + do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) + mesh_deformedCoordsFFT(1:3,i,j,k) = coords_real(i,j,k,1:3) ! ensure that data is aligned properly (fftw_alloc) + enddo; enddo; enddo + + offset_coords = matmul(F(1:3,1:3,1,1,1),step/2.0_pReal) - scaling*mesh_deformedCoordsFFT(1:3,1,1,1) + do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) + mesh_deformedCoordsFFT(1:3,i,j,k) = scaling*mesh_deformedCoordsFFT(1:3,i,j,k) & + + offset_coords + matmul(Favg,& + (/step(1)*real(i-1_pInt,pReal),& + step(2)*real(j-1_pInt,pReal),& + step(3)*real(k-1_pInt,pReal)/)) + + enddo; enddo; enddo + + call fftw_destroy_plan(fftw_forth) + call fftw_destroy_plan(fftw_back) + call fftw_free(defgrad_fftw) + call fftw_free(coords_fftw) + +end function mesh_deformedCoordsFFT + !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ subroutine volume_compare(res,geomdim,defgrad,nodes,volume_mismatch) !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ @@ -3711,28 +3862,28 @@ do e = 1_pInt,mesh_NcpElems enddo !$OMP CRITICAL (write2out) if (iand(myDebug,debug_levelBasic) /= 0_pInt) then - write (6,*) - write (6,*) 'Input Parser: STATISTICS' - write (6,*) - write (6,*) mesh_Nelems, ' : total number of elements in mesh' - write (6,*) mesh_NcpElems, ' : total number of CP elements in mesh' - write (6,*) mesh_Nnodes, ' : total number of nodes in mesh' - write (6,*) mesh_maxNnodes, ' : max number of nodes in any CP element' - write (6,*) mesh_maxNips, ' : max number of IPs in any CP element' - write (6,*) mesh_maxNipNeighbors, ' : max number of IP neighbors in any CP element' - write (6,*) mesh_maxNsubNodes, ' : max number of (additional) subnodes in any CP element' - write (6,*) mesh_maxNsharedElems, ' : max number of CP elements sharing a node' - write (6,*) - write (6,*) 'Input Parser: HOMOGENIZATION/MICROSTRUCTURE' - write (6,*) - write (6,*) mesh_maxValStateVar(1), ' : maximum homogenization index' - write (6,*) mesh_maxValStateVar(2), ' : maximum microstructure index' - write (6,*) + write(6,*) + write(6,*) 'Input Parser: STATISTICS' + write(6,*) + write(6,*) mesh_Nelems, ' : total number of elements in mesh' + write(6,*) mesh_NcpElems, ' : total number of CP elements in mesh' + write(6,*) mesh_Nnodes, ' : total number of nodes in mesh' + write(6,*) mesh_maxNnodes, ' : max number of nodes in any CP element' + write(6,*) mesh_maxNips, ' : max number of IPs in any CP element' + write(6,*) mesh_maxNipNeighbors, ' : max number of IP neighbors in any CP element' + write(6,*) mesh_maxNsubNodes, ' : max number of (additional) subnodes in any CP element' + write(6,*) mesh_maxNsharedElems, ' : max number of CP elements sharing a node' + write(6,*) + write(6,*) 'Input Parser: HOMOGENIZATION/MICROSTRUCTURE' + write(6,*) + write(6,*) mesh_maxValStateVar(1), ' : maximum homogenization index' + write(6,*) mesh_maxValStateVar(2), ' : maximum microstructure index' + write(6,*) write (myFmt,'(a,i32.32,a)') '(9x,a2,1x,',mesh_maxValStateVar(2),'(i8))' - write (6,myFmt) '+-',math_range(mesh_maxValStateVar(2)) + write(6,myFmt) '+-',math_range(mesh_maxValStateVar(2)) write (myFmt,'(a,i32.32,a)') '(i8,1x,a2,1x,',mesh_maxValStateVar(2),'(i8))' do i=1_pInt,mesh_maxValStateVar(1) ! loop over all (possibly assigned) homogenizations - write (6,myFmt) i,'| ',mesh_HomogMicro(i,:) ! loop over all (possibly assigned) microstructures + write(6,myFmt) i,'| ',mesh_HomogMicro(i,:) ! loop over all (possibly assigned) microstructures enddo write(6,*) write(6,*) 'Input Parser: ADDITIONAL MPIE OPTIONS' @@ -3743,9 +3894,9 @@ enddo endif if (iand(myDebug,debug_levelExtensive) /= 0_pInt) then - write (6,*) - write (6,*) 'Input Parser: SUBNODE COORDINATES' - write (6,*) + write(6,*) + write(6,*) 'Input Parser: SUBNODE COORDINATES' + write(6,*) write(6,'(a8,1x,a5,1x,2(a15,1x),a20,3(1x,a12))')& 'elem','IP','IP neighbor','IPFaceNodes','subNodeOnIPFace','x','y','z' do e = 1_pInt,mesh_NcpElems ! loop over cpElems @@ -3770,28 +3921,28 @@ enddo if (iand(myDebug,debug_levelSelective) /= 0_pInt .and. debug_e /= e) cycle do i = 1_pInt,FE_Nips(mesh_element(2,e)) if (iand(myDebug,debug_levelSelective) /= 0_pInt .and. debug_i /= i) cycle - write (6,'(i8,1x,i5,3(1x,f12.8))') e, i, mesh_ipCoordinates(:,i,e) + write(6,'(i8,1x,i5,3(1x,f12.8))') e, i, mesh_ipCoordinates(:,i,e) enddo enddo - write (6,*) - write (6,*) 'Input Parser: ELEMENT VOLUME' - write (6,*) - write (6,'(a13,1x,e15.8)') 'total volume', sum(mesh_ipVolume) - write (6,*) - write (6,'(a8,1x,a5,1x,a15,1x,a5,1x,a15,1x,a16)') 'elem','IP','volume','face','area','-- normal --' + write(6,*) + write(6,*) 'Input Parser: ELEMENT VOLUME' + write(6,*) + write(6,'(a13,1x,e15.8)') 'total volume', sum(mesh_ipVolume) + write(6,*) + write(6,'(a8,1x,a5,1x,a15,1x,a5,1x,a15,1x,a16)') 'elem','IP','volume','face','area','-- normal --' do e = 1_pInt,mesh_NcpElems if (iand(myDebug,debug_levelSelective) /= 0_pInt .and. debug_e /= e) cycle do i = 1_pInt,FE_Nips(mesh_element(2,e)) if (iand(myDebug,debug_levelSelective) /= 0_pInt .and. debug_i /= i) cycle - write (6,'(i8,1x,i5,1x,e15.8)') e,i,mesh_IPvolume(i,e) + write(6,'(i8,1x,i5,1x,e15.8)') e,i,mesh_IPvolume(i,e) do f = 1_pInt,FE_NipNeighbors(mesh_element(2,e)) - write (6,'(i33,1x,e15.8,1x,3(f6.3,1x))') f,mesh_ipArea(f,i,e),mesh_ipAreaNormal(:,f,i,e) + write(6,'(i33,1x,e15.8,1x,3(f6.3,1x))') f,mesh_ipArea(f,i,e),mesh_ipAreaNormal(:,f,i,e) enddo enddo enddo - write (6,*) - write (6,*) 'Input Parser: NODE TWINS' - write (6,*) + write(6,*) + write(6,*) 'Input Parser: NODE TWINS' + write(6,*) write(6,'(a6,3(3x,a6))') ' node','twin_x','twin_y','twin_z' do n = 1_pInt,mesh_Nnodes ! loop over cpNodes if (debug_e <= mesh_NcpElems) then @@ -3810,7 +3961,7 @@ enddo do i = 1_pInt,FE_Nips(t) ! loop over IPs of elem if (iand(myDebug,debug_levelSelective) /= 0_pInt .and. debug_i /= i) cycle do n = 1_pInt,FE_NipNeighbors(t) ! loop over neighbors of IP - write (6,'(i8,1x,i10,1x,i10,1x,a3,1x,i13,1x,i13)') e,i,n,'-->',mesh_ipNeighborhood(1,n,i,e),mesh_ipNeighborhood(2,n,i,e) + write(6,'(i8,1x,i10,1x,i10,1x,a3,1x,i13,1x,i13)') e,i,n,'-->',mesh_ipNeighborhood(1,n,i,e),mesh_ipNeighborhood(2,n,i,e) enddo enddo enddo