diff --git a/src/CPFEM.f90 b/src/CPFEM.f90 index 767e86851..af6ee854a 100644 --- a/src/CPFEM.f90 +++ b/src/CPFEM.f90 @@ -5,9 +5,21 @@ !-------------------------------------------------------------------------------------------------- module CPFEM use prec + use numerics + use debug + use FEsolving + use math + use mesh + use material + use config + use crystallite + use homogenization + use IO + use DAMASK_interface implicit none private + real(pReal), parameter, private :: & CPFEM_odd_stress = 1e15_pReal, & !< return value for stress in case of ping pong dummy cycle CPFEM_odd_jacobian = 1e50_pReal !< return value for jacobian in case of ping pong dummy cycle @@ -238,86 +250,6 @@ end subroutine CPFEM_init !> @brief perform initialization at first call, update variables and call the actual material model !-------------------------------------------------------------------------------------------------- subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyStress, jacobian) - use numerics, only: & - defgradTolerance, & - iJacoStiffness - use debug, only: & - debug_level, & - debug_CPFEM, & - debug_levelBasic, & - debug_levelExtensive, & - debug_levelSelective, & - debug_stressMaxLocation, & - debug_stressMinLocation, & - debug_jacobianMaxLocation, & - debug_jacobianMinLocation, & - debug_stressMax, & - debug_stressMin, & - debug_jacobianMax, & - debug_jacobianMin, & - debug_e, & - debug_i - use FEsolving, only: & - terminallyIll, & - FEsolving_execElem, & - FEsolving_execIP, & - restartWrite - use math, only: & - math_identity2nd, & - math_det33, & - math_delta, & - math_sym3333to66, & - math_66toSym3333, & - math_sym33to6, & - math_6toSym33 - use mesh, only: & - mesh_FEasCP, & - theMesh, & - mesh_element - use material, only: & - microstructure_elemhomo, & - plasticState, & - sourceState, & - homogState, & - thermalState, & - damageState, & - phaseAt, phasememberAt, & - material_phase, & - phase_plasticity, & - temperature, & - thermalMapping, & - thermal_type, & - THERMAL_conduction_ID, & - phase_Nsources, & - material_homogenizationAt - use config, only: & - material_Nhomogenization - use crystallite, only: & - crystallite_partionedF,& - crystallite_F0, & - crystallite_Fp0, & - crystallite_Fp, & - crystallite_Fi0, & - crystallite_Fi, & - crystallite_Lp0, & - crystallite_Lp, & - crystallite_Li0, & - crystallite_Li, & - crystallite_dPdF, & - crystallite_S0, & - crystallite_S - use homogenization, only: & - materialpoint_F, & - materialpoint_F0, & - materialpoint_P, & - materialpoint_dPdF, & - materialpoint_results, & - materialpoint_sizeResults, & - materialpoint_stressAndItsTangent, & - materialpoint_postResults - use IO, only: & - IO_warning - use DAMASK_interface integer(pInt), intent(in) :: elFE, & !< FE element number ip !< integration point number @@ -520,15 +452,12 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt if (.not. parallelExecution) then FEsolving_execElem(1) = elCP FEsolving_execElem(2) = elCP - if (.not. microstructure_elemhomo(mesh_element(4,elCP)) .or. & ! calculate unless homogeneous - (microstructure_elemhomo(mesh_element(4,elCP)) .and. ip == 1_pInt)) then ! and then only first ip - FEsolving_execIP(1,elCP) = ip - FEsolving_execIP(2,elCP) = ip - if (iand(debug_level(debug_CPFEM), debug_levelExtensive) /= 0_pInt) & - write(6,'(a,i8,1x,i2)') '<< CPFEM >> calculation for elFE ip ',elFE,ip - call materialpoint_stressAndItsTangent(updateJaco, dt) ! calculate stress and its tangent - call materialpoint_postResults() - endif + FEsolving_execIP(1,elCP) = ip + FEsolving_execIP(2,elCP) = ip + if (iand(debug_level(debug_CPFEM), debug_levelExtensive) /= 0_pInt) & + write(6,'(a,i8,1x,i2)') '<< CPFEM >> calculation for elFE ip ',elFE,ip + call materialpoint_stressAndItsTangent(updateJaco, dt) ! calculate stress and its tangent + call materialpoint_postResults() !* parallel computation and calulation not yet done @@ -551,13 +480,6 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt else terminalIllness - if (microstructure_elemhomo(mesh_element(4,elCP)) .and. ip > 1_pInt) then ! me homogenous? --> copy from first ip - materialpoint_P(1:3,1:3,ip,elCP) = materialpoint_P(1:3,1:3,1,elCP) - materialpoint_F(1:3,1:3,ip,elCP) = materialpoint_F(1:3,1:3,1,elCP) - materialpoint_dPdF(1:3,1:3,1:3,1:3,ip,elCP) = materialpoint_dPdF(1:3,1:3,1:3,1:3,1,elCP) - materialpoint_results(1:materialpoint_sizeResults,ip,elCP) = & - materialpoint_results(1:materialpoint_sizeResults,1,elCP) - endif ! translate from P to CS Kirchhoff = matmul(materialpoint_P(1:3,1:3,ip,elCP), transpose(materialpoint_F(1:3,1:3,ip,elCP))) diff --git a/src/geometry_plastic_nonlocal.f90 b/src/geometry_plastic_nonlocal.f90 index 0b63b7f9c..a96b6525f 100644 --- a/src/geometry_plastic_nonlocal.f90 +++ b/src/geometry_plastic_nonlocal.f90 @@ -10,8 +10,6 @@ module geometry_plastic_nonlocal implicit none private - logical, dimension(3), public, parameter :: & - geometry_plastic_nonlocal_periodicSurface = .true. !< flag indicating periodic outer surfaces (used for fluxes) NEEDED? integer, dimension(:,:,:,:), allocatable, public, protected :: & geometry_plastic_nonlocal_IPneighborhood !< 6 or less neighboring IPs as [element_num, IP_index, neighbor_index that points to me] diff --git a/src/material.f90 b/src/material.f90 index fd8f52ba9..7a9108d97 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -152,7 +152,6 @@ module material logical, dimension(:), allocatable, public, protected :: & microstructure_active, & - microstructure_elemhomo, & !< flag to indicate homogeneous microstructure distribution over element's IPs phase_localPlasticity !< flags phases with local constitutive law integer, private :: & @@ -323,12 +322,11 @@ subroutine material_init do h = 1,size(config_homogenization) write(6,'(1x,a32,1x,a16,1x,i6)') homogenization_name(h),homogenization_type(h),homogenization_Ngrains(h) enddo - write(6,'(/,a14,18x,1x,a11,1x,a12,1x,a13)') 'microstructure','crystallite','constituents','homogeneous' + write(6,'(/,a14,18x,1x,a11,1x,a12,1x,a13)') 'microstructure','crystallite','constituents' do m = 1,size(config_microstructure) - write(6,'(1x,a32,1x,i11,1x,i12,1x,l13)') microstructure_name(m), & - microstructure_crystallite(m), & - microstructure_Nconstituents(m), & - microstructure_elemhomo(m) + write(6,'(1x,a32,1x,i11,1x,i12)') microstructure_name(m), & + microstructure_crystallite(m), & + microstructure_Nconstituents(m) if (microstructure_Nconstituents(m) > 0) then do c = 1,microstructure_Nconstituents(m) write(6,'(a1,1x,a32,1x,a32,1x,f7.4)') '>',phase_name(microstructure_phase(c,m)),& @@ -538,7 +536,6 @@ subroutine material_parseMicrostructure allocate(microstructure_crystallite(size(config_microstructure)), source=0) allocate(microstructure_Nconstituents(size(config_microstructure)), source=0) allocate(microstructure_active(size(config_microstructure)), source=.false.) - allocate(microstructure_elemhomo(size(config_microstructure)), source=.false.) if(any(theMesh%microstructureAt > size(config_microstructure))) & call IO_error(155,ext_msg='More microstructures in geometry than sections in material.config') @@ -549,7 +546,6 @@ subroutine material_parseMicrostructure do m=1, size(config_microstructure) microstructure_Nconstituents(m) = config_microstructure(m)%countKeys('(constituent)') microstructure_crystallite(m) = config_microstructure(m)%getInt('crystallite') - microstructure_elemhomo(m) = config_microstructure(m)%keyExists('/elementhomogeneous/') enddo microstructure_maxNconstituents = maxval(microstructure_Nconstituents) diff --git a/src/mesh_marc.f90 b/src/mesh_marc.f90 index 683bbf3b2..a54be311e 100644 --- a/src/mesh_marc.f90 +++ b/src/mesh_marc.f90 @@ -229,7 +229,6 @@ integer, dimension(:,:), allocatable, private :: & private :: & - mesh_get_damaskOptions, & mesh_build_cellconnectivity, & mesh_build_ipAreas, & FE_mapElemtype, & @@ -341,8 +340,6 @@ subroutine mesh_init(ip,el) call mesh_marc_build_elements(FILEUNIT) if (myDebug) write(6,'(a)') ' Built elements'; flush(6) - call mesh_get_damaskOptions(mesh_periodicSurface,FILEUNIT) - if (myDebug) write(6,'(a)') ' Got DAMASK options'; flush(6) close (FILEUNIT) @@ -808,8 +805,8 @@ subroutine mesh_marc_build_elements(fileUnit) exit endif enddo - 620 rewind(fileUnit) ! just in case "initial state" appears before "connectivity" + call calcCells(theMesh,mesh_element(5:,:)) read (fileUnit,'(A300)',END=630) line do chunkPos = IO_stringPos(line) @@ -847,45 +844,6 @@ subroutine mesh_marc_build_elements(fileUnit) 630 end subroutine mesh_marc_build_elements -!-------------------------------------------------------------------------------------------------- -!> @brief get any additional damask options from input file, sets mesh_periodicSurface -!-------------------------------------------------------------------------------------------------- -subroutine mesh_get_damaskOptions(periodic_surface,fileUnit) - - integer, intent(in) :: fileUnit - - integer, allocatable, dimension(:) :: chunkPos - character(len=300) :: line - integer :: myStat - integer :: chunk, Nchunks - character(len=300) :: v - logical, dimension(3) :: periodic_surface - - - periodic_surface = .false. - myStat = 0 - rewind(fileUnit) - do while(myStat == 0) - read (fileUnit,'(a300)',iostat=myStat) line - chunkPos = IO_stringPos(line) - Nchunks = chunkPos(1) - if (IO_lc(IO_stringValue(line,chunkPos,1)) == '$damask' .and. Nchunks > 1) then ! found keyword for damask option and there is at least one more chunk to read - select case(IO_lc(IO_stringValue(line,chunkPos,2))) - case('periodic') ! damask Option that allows to specify periodic fluxes - do chunk = 3,Nchunks ! loop through chunks (skipping the keyword) - v = IO_lc(IO_stringValue(line,chunkPos,chunk)) ! chunk matches keyvalues x,y, or z? - mesh_periodicSurface(1) = mesh_periodicSurface(1) .or. v == 'x' - mesh_periodicSurface(2) = mesh_periodicSurface(2) .or. v == 'y' - mesh_periodicSurface(3) = mesh_periodicSurface(3) .or. v == 'z' - enddo - endselect - endif - enddo - -end subroutine mesh_get_damaskOptions - - - subroutine calcCells(thisMesh,connectivity_elem) class(tMesh) :: thisMesh @@ -893,11 +851,17 @@ subroutine calcCells(thisMesh,connectivity_elem) integer(pInt),dimension(:,:), allocatable :: con_elem,temp,con,parentsAndWeights,candidates_global integer(pInt),dimension(:), allocatable :: l, nodes, candidates_local integer(pInt),dimension(:,:,:), allocatable :: con_cell,connectivity_cell - integer(pInt),dimension(:,:), allocatable :: sorted,test + integer(pInt),dimension(:,:), allocatable :: sorted,test,connectivity_cell_reshape real(pReal), dimension(:,:), allocatable :: coordinates,nodes5 integer(pInt) :: e, n, c, p, s,u,i,m,j,nParentNodes,nCellNode,ierr - !nodes5 = thisMesh%node_0 +#if defined(DAMASK_HDF5) + call results_openJobFile + call HDF5_closeGroup(results_addGroup('geometry')) + call results_writeDataset('geometry',connectivity_elem,'connectivity_element',& + 'connectivity of the elements','-') +#endif + !--------------------------------------------------------------------------------------------------- ! initialize global connectivity to negative local connectivity allocate(connectivity_cell(thisMesh%elem%NcellNodesPerCell,thisMesh%elem%nIPs,thisMesh%Nelems)) @@ -1020,8 +984,13 @@ subroutine calcCells(thisMesh,connectivity_elem) if (p/=0) nodes5 = reshape([nodes5,coordinates],[3,nCellNode]) enddo thisMesh%node_0 = nodes5 - - + connectivity_cell_reshape = reshape(connectivity_cell,[thisMesh%elem%NcellNodesPerCell,thisMesh%elem%nIPs*thisMesh%Nelems]) + +#if defined(DAMASK_HDF5) + call results_writeDataset('geometry',connectivity_cell_reshape,'connectivity_cell',& + 'connectivity of the cells','-') + call results_closeJobFile +#endif end subroutine calcCells !-------------------------------------------------------------------------------------------------- diff --git a/src/plastic_nonlocal.f90 b/src/plastic_nonlocal.f90 index 66e8f8980..cddb585bb 100644 --- a/src/plastic_nonlocal.f90 +++ b/src/plastic_nonlocal.f90 @@ -8,7 +8,6 @@ module plastic_nonlocal use prec use future use geometry_plastic_nonlocal, only: & - periodicSurface => geometry_plastic_nonlocal_periodicSurface, & IPneighborhood => geometry_plastic_nonlocal_IPneighborhood, & IPvolume => geometry_plastic_nonlocal_IPvolume, & IParea => geometry_plastic_nonlocal_IParea, &