diff --git a/src/crystallite.f90 b/src/crystallite.f90 index aad4a32ac..ce82e803f 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -20,6 +20,7 @@ module crystallite use FEsolving use material use constitutive + use discretization use lattice use future use plastic_nonlocal @@ -174,8 +175,8 @@ subroutine crystallite_init write(6,'(/,a)') ' <<<+- crystallite init -+>>>' cMax = homogenization_maxNgrains - iMax = theMesh%elem%nIPs - eMax = theMesh%nElems + iMax = discretization_nIP + eMax = discretization_nElem allocate(crystallite_S0(3,3,cMax,iMax,eMax), source=0.0_pReal) allocate(crystallite_partionedS0(3,3,cMax,iMax,eMax), source=0.0_pReal) @@ -426,7 +427,6 @@ subroutine crystallite_init write(6,'(a42,1x,i10)') ' # of elements: ', eMax write(6,'(a42,1x,i10)') 'max # of integration points/element: ', iMax write(6,'(a42,1x,i10)') 'max # of constituents/integration point: ', cMax - write(6,'(a42,1x,i10)') 'max # of neigbours/integration point: ', theMesh%elem%nIPneighbors write(6,'(a42,1x,i10)') ' # of nonlocal constituents: ',count(.not. crystallite_localPlasticity) flush(6) endif @@ -443,7 +443,7 @@ end subroutine crystallite_init !-------------------------------------------------------------------------------------------------- function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC) - logical, dimension(theMesh%elem%nIPs,theMesh%Nelems) :: crystallite_stress + logical, dimension(discretization_nIP,discretization_nElem) :: crystallite_stress real(pReal), intent(in), optional :: & dummyArgumentToPreventInternalCompilerErrorWithGCC real(pReal) :: & @@ -512,7 +512,7 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC) endIP = startIP else singleRun startIP = 1 - endIP = theMesh%elem%nIPs + endIP = discretization_nIP endif singleRun NiterationCrystallite = 0 @@ -1744,11 +1744,11 @@ subroutine integrateStateAdaptiveEuler ! ToDo: MD: once all constitutives use allocate state, attach residuum arrays to the state in case of adaptive Euler real(pReal), dimension(constitutive_plasticity_maxSizeDotState, & - homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%Nelems) :: & + homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: & residuum_plastic real(pReal), dimension(constitutive_source_maxSizeDotState,& maxval(phase_Nsources), & - homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%Nelems) :: & + homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: & residuum_source !-------------------------------------------------------------------------------------------------- @@ -1919,11 +1919,11 @@ subroutine integrateStateRKCK45 ! ToDo: MD: once all constitutives use allocate state, attach residuum arrays to the state in case of RKCK45 real(pReal), dimension(constitutive_plasticity_maxSizeDotState, & - homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%Nelems) :: & + homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: & residuum_plastic ! relative residuum from evolution in microstructure real(pReal), dimension(constitutive_source_maxSizeDotState, & maxval(phase_Nsources), & - homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%Nelems) :: & + homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: & residuum_source ! relative residuum from evolution in microstructure diff --git a/src/discretization.f90 b/src/discretization.f90 index 1a4666e1b..f55fe0fcb 100644 --- a/src/discretization.f90 +++ b/src/discretization.f90 @@ -10,15 +10,15 @@ module discretization implicit none private - integer, public, protected :: & + integer, public, protected :: & discretization_nIP, & discretization_nElem - integer, dimension(:), allocatable :: & + integer, public, protected, dimension(:), allocatable :: & discretization_homogenizationAt, & discretization_microstructureAt - real(pReal), dimension(:,:), allocatable :: & + real(pReal), public, protected, dimension(:,:), allocatable :: & discretization_IPcoords0, & discretization_NodeCoords0, & discretization_IPcoords, & diff --git a/src/grid/grid_mech_FEM.f90 b/src/grid/grid_mech_FEM.f90 index 2feec76df..dd5028d48 100644 --- a/src/grid/grid_mech_FEM.f90 +++ b/src/grid/grid_mech_FEM.f90 @@ -11,13 +11,18 @@ module grid_mech_FEM use HDF5_utilities use PETScdmda use PETScsnes - use prec, only: & - pReal - use math, only: & - math_I3 - use spectral_utilities, only: & - tSolutionState, & - tSolutionParams + use prec + use CPFEM2 + use IO + use debug + use FEsolving + use numerics + use homogenization + use DAMASK_interface + use spectral_utilities + use discretization + use mesh + use math implicit none private @@ -74,30 +79,6 @@ contains !> @brief allocates all necessary fields and fills them with data, potentially from restart info !-------------------------------------------------------------------------------------------------- subroutine grid_mech_FEM_init - use IO, only: & - IO_intOut, & - IO_error, & - IO_open_jobFile_binary - use FEsolving, only: & - restartInc - use numerics, only: & - worldrank, & - worldsize, & - petsc_options - use homogenization, only: & - materialpoint_F0 - use DAMASK_interface, only: & - getSolverJobName - use spectral_utilities, only: & - utilities_constitutiveResponse, & - utilities_updateIPcoords, & - wgt - use mesh, only: & - geomSize, & - grid, & - grid3 - use math, only: & - math_invSym3333 real(pReal) :: HGCoeff = 0e-2_pReal PetscInt, dimension(:), allocatable :: localK @@ -243,14 +224,6 @@ end subroutine grid_mech_FEM_init !> @brief solution for the FEM scheme with internal iterations !-------------------------------------------------------------------------------------------------- function grid_mech_FEM_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation_BC) result(solution) - use IO, only: & - IO_error - use spectral_utilities, only: & - tBoundaryCondition, & - utilities_maskedCompliance - use FEsolving, only: & - restartWrite, & - terminallyIll !-------------------------------------------------------------------------------------------------- ! input data for solution @@ -304,25 +277,6 @@ end function grid_mech_FEM_solution !> possibly writing restart information, triggering of state increment in DAMASK, and updating of IPcoordinates !-------------------------------------------------------------------------------------------------- subroutine grid_mech_FEM_forward(guess,timeinc,timeinc_old,loadCaseTime,deformation_BC,stress_BC,rotation_BC) - use math, only: & - math_rotate_backward33 - use numerics, only: & - worldrank - use homogenization, only: & - materialpoint_F0 - use mesh, only: & - grid, & - grid3 - use CPFEM2, only: & - CPFEM_age - use spectral_utilities, only: & - utilities_updateIPcoords, & - tBoundaryCondition, & - cutBack - use IO, only: & - IO_open_jobFile_binary - use FEsolving, only: & - restartWrite logical, intent(in) :: & guess @@ -422,17 +376,6 @@ end subroutine grid_mech_FEM_forward !> @brief convergence check !-------------------------------------------------------------------------------------------------- subroutine converged(snes_local,PETScIter,devNull1,devNull2,fnorm,reason,dummy,ierr) - use mesh - use spectral_utilities - use numerics, only: & - itmax, & - itmin, & - err_div_tolRel, & - err_div_tolAbs, & - err_stress_tolRel, & - err_stress_tolAbs - use FEsolving, only: & - terminallyIll SNES :: snes_local PetscInt, intent(in) :: PETScIter @@ -481,28 +424,6 @@ end subroutine converged !-------------------------------------------------------------------------------------------------- subroutine formResidual(da_local,x_local, & f_local,dummy,ierr) - use numerics, only: & - itmax, & - itmin - use numerics, only: & - worldrank - use mesh, only: & - grid - use math, only: & - math_rotate_backward33, & - math_mul3333xx33 - use debug, only: & - debug_level, & - debug_spectral, & - debug_spectralRotation - use spectral_utilities, only: & - utilities_constitutiveResponse - use IO, only: & - IO_intOut - use FEsolving, only: & - terminallyIll - use homogenization, only: & - materialpoint_dPdF DM :: da_local Vec :: x_local, f_local @@ -617,12 +538,7 @@ end subroutine formResidual !> @brief forms the FEM stiffness matrix !-------------------------------------------------------------------------------------------------- subroutine formJacobian(da_local,x_local,Jac_pre,Jac,dummy,ierr) - use mesh, only: & - mesh_ipCoordinates - use homogenization, only: & - materialpoint_dPdF - - + DM :: da_local Vec :: x_local, coordinates Mat :: Jac_pre, Jac @@ -699,7 +615,7 @@ subroutine formJacobian(da_local,x_local,Jac_pre,Jac,dummy,ierr) ele = 0 do k = zstart, zend; do j = ystart, yend; do i = xstart, xend ele = ele + 1 - x_scal(0:2,i,j,k) = mesh_ipCoordinates(1:3,1,ele) + x_scal(0:2,i,j,k) = discretization_IPcoords(1:3,ele) enddo; enddo; enddo call DMDAVecRestoreArrayF90(da_local,coordinates,x_scal,ierr);CHKERRQ(ierr) ! initialize to undeformed coordinates (ToDo: use ip coordinates) call MatNullSpaceCreateRigidBody(coordinates,matnull,ierr);CHKERRQ(ierr) ! get rigid body deformation modes diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 9cfbb804e..99518ab31 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -14,8 +14,9 @@ module homogenization use numerics use constitutive use crystallite - use mesh use FEsolving + use mesh + use discretization use thermal_isothermal use thermal_adiabatic use thermal_conduction @@ -236,20 +237,20 @@ subroutine homogenization_init !-------------------------------------------------------------------------------------------------- ! allocate and initialize global variables - allocate(materialpoint_dPdF(3,3,3,3,theMesh%elem%nIPs,theMesh%nElems), source=0.0_pReal) - allocate(materialpoint_F0(3,3,theMesh%elem%nIPs,theMesh%nElems), source=0.0_pReal) - materialpoint_F0 = spread(spread(math_I3,3,theMesh%elem%nIPs),4,theMesh%nElems) ! initialize to identity - allocate(materialpoint_F(3,3,theMesh%elem%nIPs,theMesh%nElems), source=0.0_pReal) + allocate(materialpoint_dPdF(3,3,3,3,discretization_nIP,discretization_nElem), source=0.0_pReal) + allocate(materialpoint_F0(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 + allocate(materialpoint_F(3,3,discretization_nIP,discretization_nElem), source=0.0_pReal) materialpoint_F = materialpoint_F0 ! initialize to identity - allocate(materialpoint_subF0(3,3,theMesh%elem%nIPs,theMesh%nElems), source=0.0_pReal) - allocate(materialpoint_subF(3,3,theMesh%elem%nIPs,theMesh%nElems), source=0.0_pReal) - allocate(materialpoint_P(3,3,theMesh%elem%nIPs,theMesh%nElems), source=0.0_pReal) - allocate(materialpoint_subFrac(theMesh%elem%nIPs,theMesh%nElems), source=0.0_pReal) - allocate(materialpoint_subStep(theMesh%elem%nIPs,theMesh%nElems), source=0.0_pReal) - allocate(materialpoint_subdt(theMesh%elem%nIPs,theMesh%nElems), source=0.0_pReal) - allocate(materialpoint_requested(theMesh%elem%nIPs,theMesh%nElems), source=.false.) - allocate(materialpoint_converged(theMesh%elem%nIPs,theMesh%nElems), source=.true.) - allocate(materialpoint_doneAndHappy(2,theMesh%elem%nIPs,theMesh%nElems), source=.true.) + allocate(materialpoint_subF0(3,3,discretization_nIP,discretization_nElem), source=0.0_pReal) + allocate(materialpoint_subF(3,3,discretization_nIP,discretization_nElem), source=0.0_pReal) + allocate(materialpoint_P(3,3,discretization_nIP,discretization_nElem), source=0.0_pReal) + allocate(materialpoint_subFrac(discretization_nIP,discretization_nElem), source=0.0_pReal) + allocate(materialpoint_subStep(discretization_nIP,discretization_nElem), source=0.0_pReal) + allocate(materialpoint_subdt(discretization_nIP,discretization_nElem), source=0.0_pReal) + allocate(materialpoint_requested(discretization_nIP,discretization_nElem), source=.false.) + allocate(materialpoint_converged(discretization_nIP,discretization_nElem), source=.true.) + allocate(materialpoint_doneAndHappy(2,discretization_nIP,discretization_nElem), source=.true.) !-------------------------------------------------------------------------------------------------- ! allocate and initialize global state and postresutls variables @@ -266,7 +267,7 @@ subroutine homogenization_init + homogenization_maxNgrains * (1 + crystallite_maxSizePostResults & ! crystallite size & crystallite results + 1 + constitutive_plasticity_maxSizePostResults & ! constitutive size & constitutive results + constitutive_source_maxSizePostResults) - allocate(materialpoint_results(materialpoint_sizeResults,theMesh%elem%nIPs,theMesh%nElems)) + allocate(materialpoint_results(materialpoint_sizeResults,discretization_nIP,discretization_nElem)) write(6,'(/,a)') ' <<<+- homogenization init -+>>>' @@ -286,7 +287,7 @@ subroutine homogenization_init endif flush(6) - if (debug_g < 1 .or. debug_g > homogenization_Ngrains(mesh_element(3,debug_e))) & + if (debug_g < 1 .or. debug_g > homogenization_Ngrains(material_homogenizationAt(debug_e))) & call IO_error(602,ext_msg='constituent', el=debug_e, g=debug_g) end subroutine homogenization_init @@ -601,7 +602,7 @@ subroutine materialpoint_postResults !$OMP PARALLEL DO PRIVATE(myNgrains,myCrystallite,thePos,theSize) elementLooping: do e = FEsolving_execElem(1),FEsolving_execElem(2) myNgrains = homogenization_Ngrains(material_homogenizationAt(e)) - myCrystallite = microstructure_crystallite(mesh_element(4,e)) + myCrystallite = microstructure_crystallite(discretization_microstructureAt(e)) IpLooping: do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) thePos = 0