diff --git a/CMakeLists.txt b/CMakeLists.txt index 8db6dd0c0..198528b59 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -4,7 +4,7 @@ include (FindPkgConfig REQUIRED) # Dummy project to determine compiler names and version project (Prerequisites LANGUAGES) set(ENV{PKG_CONFIG_PATH} "$ENV{PETSC_DIR}/$ENV{PETSC_ARCH}/lib/pkgconfig") -pkg_search_module (PETSC REQUIRED PETSc>3.12.0) +pkg_check_modules (PETSC REQUIRED PETSc>=3.12.0 PETSc<3.15.0) pkg_get_variable (CMAKE_Fortran_COMPILER PETSc fcompiler) pkg_get_variable (CMAKE_C_COMPILER PETSc ccompiler) diff --git a/VERSION b/VERSION index 616042312..62c706093 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v3.0.0-alpha2-97-g10bbeb561 +v3.0.0-alpha2-173-g584c7cc3a diff --git a/examples/FEM/polyXtal/material.yaml b/examples/FEM/polyXtal/material.yaml index c7d17657d..333073150 100644 --- a/examples/FEM/polyXtal/material.yaml +++ b/examples/FEM/polyXtal/material.yaml @@ -5,8 +5,8 @@ homogenization: phase: Aluminum: + lattice: cF mechanics: - lattice: cF output: [F, P, F_e, F_p, L_p] elasticity: {C_11: 106.75e9, C_12: 60.41e9, C_44: 28.34e9, type: hooke} plasticity: diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index d68cea6d3..ace2a3dba 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -367,7 +367,7 @@ class Rotation: """Intermediate representation supporting quaternion averaging.""" return np.einsum('...i,...j',quat,quat) - if not weights: + if weights is None: weights = np.ones(self.shape,dtype=float) eig, vec = np.linalg.eig(np.sum(_M(self.quaternion) * weights[...,np.newaxis,np.newaxis],axis=-3) \ diff --git a/python/tests/test_Rotation.py b/python/tests/test_Rotation.py index 707bc0210..1067e8a87 100644 --- a/python/tests/test_Rotation.py +++ b/python/tests/test_Rotation.py @@ -825,6 +825,14 @@ class TestRotation: print(f'append 3x {shape} --> {s.shape}') assert np.logical_and(s[0,...] == r[0,...], s[-1,...] == p[-1,...]).all() + @pytest.mark.parametrize('shape',[None,1,(1,),(4,2),(3,3,2)]) + def test_append_list(self,shape): + r = Rotation.from_random(shape=shape) + p = Rotation.from_random(shape=shape) + s = r.append([r,p]) + print(f'append 3x {shape} --> {s.shape}') + assert s[0,...] == r[0,...] and s[-1,...] == p[-1,...] + @pytest.mark.parametrize('quat,standardized',[ ([-1,0,0,0],[1,0,0,0]), ([-0.5,-0.5,-0.5,-0.5],[0.5,0.5,0.5,0.5]), diff --git a/src/CPFEM.f90 b/src/CPFEM.f90 index fe8c7d1b3..240688a8c 100644 --- a/src/CPFEM.f90 +++ b/src/CPFEM.f90 @@ -5,7 +5,6 @@ !-------------------------------------------------------------------------------------------------- module CPFEM use prec - use FEsolving use math use rotations use YAML_types @@ -197,11 +196,9 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS CPFEM_dcsde(1:6,1:6,ip,elCP) = ODD_JACOBIAN * math_eye(6) else validCalculation - FEsolving_execElem = elCP - FEsolving_execIP = ip if (debugCPFEM%extensive) & print'(a,i8,1x,i2)', '<< CPFEM >> calculation for elFE ip ',elFE,ip - call materialpoint_stressAndItsTangent(dt) + call materialpoint_stressAndItsTangent(dt,[ip,ip],[elCP,elCP]) terminalIllness: if (terminallyIll) then @@ -260,7 +257,7 @@ end subroutine CPFEM_general !-------------------------------------------------------------------------------------------------- subroutine CPFEM_forward - call crystallite_forward + call homogenization_forward call constitutive_forward end subroutine CPFEM_forward @@ -277,7 +274,6 @@ subroutine CPFEM_results(inc,time) call results_openJobFile call results_addIncrement(inc,time) call constitutive_results - call crystallite_results call homogenization_results call discretization_results call results_finalizeIncrement diff --git a/src/CPFEM2.f90 b/src/CPFEM2.f90 index 636962948..5a500875d 100644 --- a/src/CPFEM2.f90 +++ b/src/CPFEM2.f90 @@ -6,7 +6,6 @@ module CPFEM2 use prec use config - use FEsolving use math use rotations use YAML_types @@ -97,7 +96,7 @@ end subroutine CPFEM_restartWrite !-------------------------------------------------------------------------------------------------- subroutine CPFEM_forward - call crystallite_forward + call homogenization_forward call constitutive_forward end subroutine CPFEM_forward @@ -114,7 +113,6 @@ subroutine CPFEM_results(inc,time) call results_openJobFile call results_addIncrement(inc,time) call constitutive_results - call crystallite_results call homogenization_results call discretization_results call results_finalizeIncrement diff --git a/src/DAMASK_interface.f90 b/src/DAMASK_interface.f90 index d38020225..c43e354a2 100644 --- a/src/DAMASK_interface.f90 +++ b/src/DAMASK_interface.f90 @@ -54,12 +54,6 @@ subroutine DAMASK_interface_init =================================================================================================== -- WRONG PETSc VERSION --- WRONG PETSc VERSION --- WRONG PETSc VERSION --- WRONG PETSc VERSION -- =================================================================================================== -============ THIS VERSION OF DAMASK REQUIRES A DIFFERENT PETSc VERSION ======================== -=============== THIS VERSION OF DAMASK REQUIRES A DIFFERENT PETSc VERSION ===================== -================== THIS VERSION OF DAMASK REQUIRES A DIFFERENT PETSc VERSION ================== -=================================================================================================== --- WRONG PETSc VERSION --- WRONG PETSc VERSION --- WRONG PETSc VERSION --- WRONG PETSc VERSION -- -=================================================================================================== #endif character(len=pPathLen*3+pStringLen) :: & diff --git a/src/DAMASK_marc.f90 b/src/DAMASK_marc.f90 index ea7430c6b..0ad68445c 100644 --- a/src/DAMASK_marc.f90 +++ b/src/DAMASK_marc.f90 @@ -176,7 +176,6 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, & use DAMASK_interface use config use YAML_types - use FEsolving use discretization_marc use homogenization use CPFEM diff --git a/src/FEsolving.f90 b/src/FEsolving.f90 deleted file mode 100644 index 3fc1482d3..000000000 --- a/src/FEsolving.f90 +++ /dev/null @@ -1,15 +0,0 @@ -!-------------------------------------------------------------------------------------------------- -!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH -!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH -!> @brief global variables for flow control -!-------------------------------------------------------------------------------------------------- -module FEsolving - - implicit none - public - - integer, dimension(2) :: & - FEsolving_execElem, & !< for ping-pong scheme always whole range, otherwise one specific element - FEsolving_execIP !< for ping-pong scheme always range to max IP, otherwise one specific IP - -end module FEsolving diff --git a/src/commercialFEM_fileList.f90 b/src/commercialFEM_fileList.f90 index 08e7b9c1c..674ec4c43 100644 --- a/src/commercialFEM_fileList.f90 +++ b/src/commercialFEM_fileList.f90 @@ -11,9 +11,7 @@ #include "config.f90" #include "LAPACK_interface.f90" #include "math.f90" -#include "quaternions.f90" #include "rotations.f90" -#include "FEsolving.f90" #include "element.f90" #include "HDF5_utilities.f90" #include "results.f90" diff --git a/src/constitutive.f90 b/src/constitutive.f90 index f76b94579..696611549 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -16,55 +16,51 @@ module constitutive use parallelization use HDF5_utilities use DAMASK_interface - use FEsolving use results implicit none private - real(pReal), dimension(:,:,:), allocatable, public :: & - crystallite_dt !< requested time increment of each grain + + enum, bind(c); enumerator :: & + PLASTICITY_UNDEFINED_ID, & + PLASTICITY_NONE_ID, & + PLASTICITY_ISOTROPIC_ID, & + PLASTICITY_PHENOPOWERLAW_ID, & + PLASTICITY_KINEHARDENING_ID, & + PLASTICITY_DISLOTWIN_ID, & + PLASTICITY_DISLOTUNGSTEN_ID, & + PLASTICITY_NONLOCAL_ID, & + SOURCE_UNDEFINED_ID ,& + SOURCE_THERMAL_DISSIPATION_ID, & + SOURCE_THERMAL_EXTERNALHEAT_ID, & + SOURCE_DAMAGE_ISOBRITTLE_ID, & + SOURCE_DAMAGE_ISODUCTILE_ID, & + SOURCE_DAMAGE_ANISOBRITTLE_ID, & + SOURCE_DAMAGE_ANISODUCTILE_ID, & + KINEMATICS_UNDEFINED_ID ,& + KINEMATICS_CLEAVAGE_OPENING_ID, & + KINEMATICS_SLIPPLANE_OPENING_ID, & + KINEMATICS_THERMAL_EXPANSION_ID + end enum real(pReal), dimension(:,:,:), allocatable :: & - crystallite_subdt, & !< substepped time increment of each grain - crystallite_subStep !< size of next integration step + crystallite_subdt !< substepped time increment of each grain type(rotation), dimension(:,:,:), allocatable :: & crystallite_orientation !< current orientation real(pReal), dimension(:,:,:,:,:), allocatable :: & crystallite_F0, & !< def grad at start of FE inc - crystallite_subF, & !< def grad to be reached at end of crystallite inc - crystallite_subF0, & !< def grad at start of crystallite inc - ! crystallite_Fe, & !< current "elastic" def grad (end of converged time step) - ! - crystallite_Fp, & !< current plastic def grad (end of converged time step) - crystallite_Fp0, & !< plastic def grad at start of FE inc - crystallite_partitionedFp0,& !< plastic def grad at start of homog inc crystallite_subFp0,& !< plastic def grad at start of crystallite inc - ! crystallite_subFi0,& !< intermediate def grad at start of crystallite inc - ! crystallite_Lp0, & !< plastic velocitiy grad at start of FE inc crystallite_partitionedLp0, & !< plastic velocity grad at start of homog inc - ! crystallite_S0, & !< 2nd Piola-Kirchhoff stress vector at start of FE inc crystallite_partitionedS0 !< 2nd Piola-Kirchhoff stress vector at start of homog inc real(pReal), dimension(:,:,:,:,:), allocatable, public :: & crystallite_P, & !< 1st Piola-Kirchhoff stress per grain crystallite_Lp, & !< current plastic velocitiy grad (end of converged time step) crystallite_S, & !< current 2nd Piola-Kirchhoff stress vector (end of converged time step) - crystallite_partitionedF0 !< def grad at start of homog inc - real(pReal), dimension(:,:,:,:,:), allocatable, public :: & - crystallite_partitionedF !< def grad to be reached at end of homog inc - - logical, dimension(:,:,:), allocatable, public :: & - crystallite_requested !< used by upper level (homogenization) to request crystallite calculation - logical, dimension(:,:,:), allocatable :: & - crystallite_converged !< convergence flag - - type :: tOutput !< new requested output (per phase) - character(len=pStringLen), allocatable, dimension(:) :: & - label - end type tOutput - type(tOutput), allocatable, dimension(:) :: output_constituent + crystallite_partitionedF0, & !< def grad at start of homog inc + crystallite_F !< def grad to be reached at end of homog inc type :: tTensorContainer real(pReal), dimension(:,:,:), allocatable :: data @@ -73,10 +69,14 @@ module constitutive type(tTensorContainer), dimension(:), allocatable :: & constitutive_mech_Fi, & constitutive_mech_Fi0, & - constitutive_mech_partionedFi0, & + constitutive_mech_partitionedFi0, & constitutive_mech_Li, & constitutive_mech_Li0, & - constitutive_mech_partionedLi0 + constitutive_mech_partitionedLi0, & + constitutive_mech_Fp, & + constitutive_mech_Fp0, & + constitutive_mech_partitionedFp0 + type :: tNumerics integer :: & @@ -109,9 +109,9 @@ module constitutive type(tDebugOptions) :: debugCrystallite - procedure(integrateStateFPI), pointer :: integrateState - integer(kind(PLASTICITY_undefined_ID)), dimension(:), allocatable :: & + + integer(kind(PLASTICITY_undefined_ID)), dimension(:), allocatable, public :: & phase_plasticity !< plasticity of each phase integer(kind(SOURCE_undefined_ID)), dimension(:,:), allocatable :: & @@ -140,6 +140,7 @@ module constitutive interface +! == cleaned:begin ================================================================================= module subroutine mech_init end subroutine mech_init @@ -149,83 +150,73 @@ module constitutive module subroutine thermal_init end subroutine thermal_init - module function constitutive_collectDotState(S, FArray, Fi, FpArray, subdt, ipc, ip, el,phase,of) result(broken) - integer, intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el, & !< element - phase, & - of - real(pReal), intent(in) :: & - subdt !< timestep - real(pReal), intent(in), dimension(3,3,homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems) :: & - FArray, & !< elastic deformation gradient - FpArray !< plastic deformation gradient - real(pReal), intent(in), dimension(3,3) :: & - Fi !< intermediate deformation gradient - real(pReal), intent(in), dimension(3,3) :: & - S !< 2nd Piola Kirchhoff stress (vector notation) + module subroutine mech_results(group,ph) + character(len=*), intent(in) :: group + integer, intent(in) :: ph + end subroutine mech_results - logical :: broken -end function constitutive_collectDotState + module subroutine damage_results(group,ph) + character(len=*), intent(in) :: group + integer, intent(in) :: ph + end subroutine damage_results -module function constitutive_deltaState(S, Fi, ipc, ip, el, phase, of) result(broken) + module subroutine mech_restart_read(fileHandle) + integer(HID_T), intent(in) :: fileHandle + end subroutine mech_restart_read - integer, intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el, & !< element - phase, & - of - real(pReal), intent(in), dimension(3,3) :: & - S, & !< 2nd Piola Kirchhoff stress - Fi !< intermediate deformation gradient - logical :: & - broken + module subroutine mech_initializeRestorationPoints(ph,me) + integer, intent(in) :: ph, me + end subroutine mech_initializeRestorationPoints + module subroutine constitutive_mech_windForward(ph,me) + integer, intent(in) :: ph, me + end subroutine constitutive_mech_windForward -end function constitutive_deltaState + module subroutine constitutive_mech_forward + end subroutine constitutive_mech_forward - - module function plastic_active(plastic_label) result(active_plastic) - character(len=*), intent(in) :: plastic_label - logical, dimension(:), allocatable :: active_plastic - end function plastic_active - - module function source_active(source_label,src_length) result(active_source) - character(len=*), intent(in) :: source_label - integer, intent(in) :: src_length - logical, dimension(:,:), allocatable :: active_source - end function source_active - - module function kinematics_active(kinematics_label,kinematics_length) result(active_kinematics) - character(len=*), intent(in) :: kinematics_label - integer, intent(in) :: kinematics_length - logical, dimension(:,:), allocatable :: active_kinematics - end function kinematics_active - - - module subroutine source_damage_anisoBrittle_dotState(S, ipc, ip, el) + module subroutine mech_restore(ip,el,includeL) integer, intent(in) :: & - ipc, & !< component-ID of integration point + ip, & + el + logical, intent(in) :: & + includeL + end subroutine mech_restore + +! == cleaned:end =================================================================================== + + module function crystallite_stress(dt,co,ip,el) result(converged_) + real(pReal), intent(in) :: dt + integer, intent(in) :: co, ip, el + logical :: converged_ + end function crystallite_stress + + module function constitutive_homogenizedC(co,ip,el) result(C) + integer, intent(in) :: co, ip, el + real(pReal), dimension(6,6) :: C + end function constitutive_homogenizedC + + module subroutine source_damage_anisoBrittle_dotState(S, co, ip, el) + integer, intent(in) :: & + co, & !< component-ID of integration point ip, & !< integration point el !< element real(pReal), intent(in), dimension(3,3) :: & S end subroutine source_damage_anisoBrittle_dotState - module subroutine source_damage_anisoDuctile_dotState(ipc, ip, el) + module subroutine source_damage_anisoDuctile_dotState(co, ip, el) integer, intent(in) :: & - ipc, & !< component-ID of integration point + co, & !< component-ID of integration point ip, & !< integration point el !< element end subroutine source_damage_anisoDuctile_dotState - module subroutine source_damage_isoDuctile_dotState(ipc, ip, el) + module subroutine source_damage_isoDuctile_dotState(co, ip, el) integer, intent(in) :: & - ipc, & !< component-ID of integration point + co, & !< component-ID of integration point ip, & !< integration point el !< element end subroutine source_damage_isoDuctile_dotState @@ -261,14 +252,7 @@ end function constitutive_deltaState dTDot_dT end subroutine constitutive_thermal_getRateAndItsTangents - module function plastic_dislotwin_homogenizedC(ipc,ip,el) result(homogenizedC) - real(pReal), dimension(6,6) :: & - homogenizedC - integer, intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el !< element - end function plastic_dislotwin_homogenizedC + module subroutine plastic_nonlocal_updateCompatibility(orientation,instance,i,e) integer, intent(in) :: & @@ -291,9 +275,9 @@ end function constitutive_deltaState of end subroutine plastic_isotropic_LiAndItsTangent - module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, ip, el) + module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, co, ip, el) integer, intent(in) :: & - ipc, & !< grain number + co, & !< grain number ip, & !< integration point number el !< element number real(pReal), intent(in), dimension(3,3) :: & @@ -304,9 +288,9 @@ end function constitutive_deltaState dLd_dTstar !< derivative of Ld with respect to Tstar (4th-order tensor) end subroutine kinematics_cleavage_opening_LiAndItsTangent - module subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, ip, el) + module subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, co, ip, el) integer, intent(in) :: & - ipc, & !< grain number + co, & !< grain number ip, & !< integration point number el !< element number real(pReal), intent(in), dimension(3,3) :: & @@ -317,9 +301,9 @@ end function constitutive_deltaState dLd_dTstar !< derivative of Ld with respect to Tstar (4th-order tensor) end subroutine kinematics_slipplane_opening_LiAndItsTangent - module subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar, ipc, ip, el) + module subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar, co, ip, el) integer, intent(in) :: & - ipc, & !< grain number + co, & !< grain number ip, & !< integration point number el !< element number real(pReal), intent(out), dimension(3,3) :: & @@ -329,9 +313,9 @@ end function constitutive_deltaState end subroutine kinematics_thermal_expansion_LiAndItsTangent - module subroutine source_damage_isoBrittle_deltaState(C, Fe, ipc, ip, el) + module subroutine source_damage_isoBrittle_deltaState(C, Fe, co, ip, el) integer, intent(in) :: & - ipc, & !< component-ID of integration point + co, & !< component-ID of integration point ip, & !< integration point el !< element real(pReal), intent(in), dimension(3,3) :: & @@ -340,18 +324,11 @@ end function constitutive_deltaState C end subroutine source_damage_isoBrittle_deltaState - module subroutine plastic_results - end subroutine plastic_results - - module subroutine damage_results - end subroutine damage_results - - module subroutine constitutive_plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & - S, Fi, ipc, ip, el) + S, Fi, co, ip, el) integer, intent(in) :: & - ipc, & !< component-ID of integration point + co, & !< component-ID of integration point ip, & !< integration point el !< element real(pReal), intent(in), dimension(3,3) :: & @@ -365,21 +342,20 @@ end function constitutive_deltaState end subroutine constitutive_plastic_LpAndItsTangents - module subroutine constitutive_plastic_dependentState(F, Fp, ipc, ip, el) + module subroutine constitutive_plastic_dependentState(F, co, ip, el) integer, intent(in) :: & - ipc, & !< component-ID of integration point + co, & !< component-ID of integration point ip, & !< integration point el !< element real(pReal), intent(in), dimension(3,3) :: & - F, & !< elastic deformation gradient - Fp !< plastic deformation gradient + F !< elastic deformation gradient end subroutine constitutive_plastic_dependentState - module subroutine constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, Fe, Fi, ipc, ip, el) + module subroutine constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, Fe, Fi, co, ip, el) integer, intent(in) :: & - ipc, & !< component-ID of integration point + co, & !< component-ID of integration point ip, & !< integration point el !< element real(pReal), intent(in), dimension(3,3) :: & @@ -392,39 +368,16 @@ end function constitutive_deltaState dS_dFi !< derivative of 2nd P-K stress with respect to intermediate deformation gradient end subroutine constitutive_hooke_SandItsTangents - module subroutine integrateStateFPI(g,i,e) - integer, intent(in) :: e, i, g - end subroutine integrateStateFPI - - module subroutine integrateStateEuler(g,i,e) - integer, intent(in) :: e, i, g - end subroutine integrateStateEuler - - module subroutine integrateStateAdaptiveEuler(g,i,e) - integer, intent(in) :: e, i, g - end subroutine integrateStateAdaptiveEuler - - module subroutine integrateStateRK4(g,i,e) - integer, intent(in) :: e, i, g - end subroutine integrateStateRK4 - - module subroutine integrateStateRKCK45(g,i,e) - integer, intent(in) :: e, i, g - end subroutine integrateStateRKCK45 - end interface + type(tDebugOptions) :: debugConstitutive public :: & constitutive_init, & constitutive_homogenizedC, & constitutive_LiAndItsTangents, & - constitutive_collectDotState, & - constitutive_collectDotState_source, & - constitutive_deltaState, & - constitutive_deltaState_source, & constitutive_damage_getRateAndItsTangents, & constitutive_thermal_getRateAndItsTangents, & constitutive_results, & @@ -432,7 +385,6 @@ end function constitutive_deltaState constitutive_forward, & constitutive_restore, & plastic_nonlocal_updateCompatibility, & - plastic_active, & source_active, & kinematics_active, & converged, & @@ -441,13 +393,30 @@ end function constitutive_deltaState crystallite_stressTangent, & crystallite_orientations, & crystallite_push33ToRef, & - crystallite_results, & crystallite_restartWrite, & + integrateSourceState, & crystallite_restartRead, & - crystallite_forward, & - crystallite_initializeRestorationPoints, & - crystallite_windForward, & - crystallite_restore + constitutive_initializeRestorationPoints, & + constitutive_windForward, & + PLASTICITY_UNDEFINED_ID, & + PLASTICITY_NONE_ID, & + PLASTICITY_ISOTROPIC_ID, & + PLASTICITY_PHENOPOWERLAW_ID, & + PLASTICITY_KINEHARDENING_ID, & + PLASTICITY_DISLOTWIN_ID, & + PLASTICITY_DISLOTUNGSTEN_ID, & + PLASTICITY_NONLOCAL_ID, & + SOURCE_UNDEFINED_ID ,& + SOURCE_THERMAL_DISSIPATION_ID, & + SOURCE_THERMAL_EXTERNALHEAT_ID, & + SOURCE_DAMAGE_ISOBRITTLE_ID, & + SOURCE_DAMAGE_ISODUCTILE_ID, & + SOURCE_DAMAGE_ANISOBRITTLE_ID, & + SOURCE_DAMAGE_ANISODUCTILE_ID, & + KINEMATICS_UNDEFINED_ID ,& + KINEMATICS_CLEAVAGE_OPENING_ID, & + KINEMATICS_SLIPPLANE_OPENING_ID, & + KINEMATICS_THERMAL_EXPANSION_ID contains @@ -458,12 +427,13 @@ contains subroutine constitutive_init integer :: & - p, & !< counter in phase loop - s !< counter in source loop + ph, & !< counter in phase loop + so !< counter in source loop class (tNode), pointer :: & debug_constitutive, & phases + debug_constitutive => config_debug%get('constitutive', defaultVal=emptyList) debugConstitutive%basic = debug_constitutive%contains('basic') debugConstitutive%extensive = debug_constitutive%contains('extensive') @@ -474,26 +444,26 @@ subroutine constitutive_init !-------------------------------------------------------------------------------------------------- ! initialize constitutive laws + print'(/,a)', ' <<<+- constitutive init -+>>>'; flush(IO_STDOUT) call mech_init call damage_init call thermal_init - print'(/,a)', ' <<<+- constitutive init -+>>>'; flush(IO_STDOUT) phases => config_material%get('phase') constitutive_source_maxSizeDotState = 0 - PhaseLoop2:do p = 1,phases%length + PhaseLoop2:do ph = 1,phases%length !-------------------------------------------------------------------------------------------------- ! partition and initialize state - plasticState(p)%partitionedState0 = plasticState(p)%state0 - plasticState(p)%state = plasticState(p)%partitionedState0 - forall(s = 1:phase_Nsources(p)) - sourceState(p)%p(s)%partitionedState0 = sourceState(p)%p(s)%state0 - sourceState(p)%p(s)%state = sourceState(p)%p(s)%partitionedState0 + plasticState(ph)%partitionedState0 = plasticState(ph)%state0 + plasticState(ph)%state = plasticState(ph)%partitionedState0 + forall(so = 1:phase_Nsources(ph)) + sourceState(ph)%p(so)%partitionedState0 = sourceState(ph)%p(so)%state0 + sourceState(ph)%p(so)%state = sourceState(ph)%p(so)%partitionedState0 end forall constitutive_source_maxSizeDotState = max(constitutive_source_maxSizeDotState, & - maxval(sourceState(p)%p%sizeDotState)) + maxval(sourceState(ph)%p%sizeDotState)) enddo PhaseLoop2 constitutive_plasticity_maxSizeDotState = maxval(plasticState%sizeDotState) @@ -503,7 +473,7 @@ end subroutine constitutive_init !-------------------------------------------------------------------------------------------------- !> @brief checks if a source mechanism is active or not !-------------------------------------------------------------------------------------------------- -module function source_active(source_label,src_length) result(active_source) +function source_active(source_label,src_length) result(active_source) character(len=*), intent(in) :: source_label !< name of source mechanism integer, intent(in) :: src_length !< max. number of sources in system @@ -534,8 +504,7 @@ end function source_active !-------------------------------------------------------------------------------------------------- !> @brief checks if a kinematic mechanism is active or not !-------------------------------------------------------------------------------------------------- - -module function kinematics_active(kinematics_label,kinematics_length) result(active_kinematics) +function kinematics_active(kinematics_label,kinematics_length) result(active_kinematics) character(len=*), intent(in) :: kinematics_label !< name of kinematic mechanism integer, intent(in) :: kinematics_length !< max. number of kinematics in system @@ -563,38 +532,15 @@ module function kinematics_active(kinematics_label,kinematics_length) result(ac end function kinematics_active -!-------------------------------------------------------------------------------------------------- -!> @brief returns the homogenize elasticity matrix -!> ToDo: homogenizedC66 would be more consistent -!-------------------------------------------------------------------------------------------------- -function constitutive_homogenizedC(ipc,ip,el) - - real(pReal), dimension(6,6) :: & - constitutive_homogenizedC - integer, intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el !< element - - plasticityType: select case (phase_plasticity(material_phaseAt(ipc,el))) - case (PLASTICITY_DISLOTWIN_ID) plasticityType - constitutive_homogenizedC = plastic_dislotwin_homogenizedC(ipc,ip,el) - case default plasticityType - constitutive_homogenizedC = lattice_C66(1:6,1:6,material_phaseAt(ipc,el)) - end select plasticityType - -end function constitutive_homogenizedC - - !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the velocity gradient ! ToDo: MD: S is Mi? !-------------------------------------------------------------------------------------------------- subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & - S, Fi, ipc, ip, el) + S, Fi, co, ip, el) integer, intent(in) :: & - ipc, & !< component-ID of integration point + co, & !< component-ID of integration point ip, & !< integration point el !< element real(pReal), intent(in), dimension(3,3) :: & @@ -623,10 +569,10 @@ subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & dLi_dS = 0.0_pReal dLi_dFi = 0.0_pReal - plasticityType: select case (phase_plasticity(material_phaseAt(ipc,el))) + plasticityType: select case (phase_plasticity(material_phaseAt(co,el))) case (PLASTICITY_isotropic_ID) plasticityType - of = material_phasememberAt(ipc,ip,el) - instance = phase_plasticityInstance(material_phaseAt(ipc,el)) + of = material_phasememberAt(co,ip,el) + instance = phase_plasticityInstance(material_phaseAt(co,el)) call plastic_isotropic_LiAndItsTangent(my_Li, my_dLi_dS, S ,instance,of) case default plasticityType my_Li = 0.0_pReal @@ -636,14 +582,14 @@ subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & Li = Li + my_Li dLi_dS = dLi_dS + my_dLi_dS - KinematicsLoop: do k = 1, phase_Nkinematics(material_phaseAt(ipc,el)) - kinematicsType: select case (phase_kinematics(k,material_phaseAt(ipc,el))) + KinematicsLoop: do k = 1, phase_Nkinematics(material_phaseAt(co,el)) + kinematicsType: select case (phase_kinematics(k,material_phaseAt(co,el))) case (KINEMATICS_cleavage_opening_ID) kinematicsType - call kinematics_cleavage_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, ipc, ip, el) + call kinematics_cleavage_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, co, ip, el) case (KINEMATICS_slipplane_opening_ID) kinematicsType - call kinematics_slipplane_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, ipc, ip, el) + call kinematics_slipplane_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, co, ip, el) case (KINEMATICS_thermal_expansion_ID) kinematicsType - call kinematics_thermal_expansion_LiAndItsTangent(my_Li, my_dLi_dS, ipc, ip, el) + call kinematics_thermal_expansion_LiAndItsTangent(my_Li, my_dLi_dS, co, ip, el) case default kinematicsType my_Li = 0.0_pReal my_dLi_dS = 0.0_pReal @@ -666,70 +612,89 @@ subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & end subroutine constitutive_LiAndItsTangents - - - !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the rate of change of microstructure !-------------------------------------------------------------------------------------------------- -function constitutive_collectDotState_source(S, ipc, ip, el,phase,of) result(broken) +function constitutive_damage_collectDotState(S, co, ip, el,ph,of) result(broken) integer, intent(in) :: & - ipc, & !< component-ID of integration point + co, & !< component-ID of integration point ip, & !< integration point el, & !< element - phase, & + ph, & of real(pReal), intent(in), dimension(3,3) :: & S !< 2nd Piola Kirchhoff stress (vector notation) integer :: & - i !< counter in source loop + so !< counter in source loop logical :: broken broken = .false. - SourceLoop: do i = 1, phase_Nsources(phase) + SourceLoop: do so = 1, phase_Nsources(ph) - sourceType: select case (phase_source(i,phase)) + sourceType: select case (phase_source(so,ph)) case (SOURCE_damage_anisoBrittle_ID) sourceType - call source_damage_anisoBrittle_dotState(S, ipc, ip, el) ! correct stress? + call source_damage_anisoBrittle_dotState(S, co, ip, el) ! correct stress? case (SOURCE_damage_isoDuctile_ID) sourceType - call source_damage_isoDuctile_dotState(ipc, ip, el) + call source_damage_isoDuctile_dotState(co, ip, el) case (SOURCE_damage_anisoDuctile_ID) sourceType - call source_damage_anisoDuctile_dotState(ipc, ip, el) - - case (SOURCE_thermal_externalheat_ID) sourceType - call source_thermal_externalheat_dotState(phase,of) + call source_damage_anisoDuctile_dotState(co, ip, el) end select sourceType - broken = broken .or. any(IEEE_is_NaN(sourceState(phase)%p(i)%dotState(:,of))) + broken = broken .or. any(IEEE_is_NaN(sourceState(ph)%p(so)%dotState(:,of))) enddo SourceLoop -end function constitutive_collectDotState_source +end function constitutive_damage_collectDotState + + +!-------------------------------------------------------------------------------------------------- +!> @brief contains the constitutive equation for calculating the rate of change of microstructure +!-------------------------------------------------------------------------------------------------- +function constitutive_thermal_collectDotState(ph,me) result(broken) + + integer, intent(in) :: ph, me + logical :: broken + + integer :: i + + + broken = .false. + + SourceLoop: do i = 1, phase_Nsources(ph) + + if (phase_source(i,ph) == SOURCE_thermal_externalheat_ID) & + call source_thermal_externalheat_dotState(ph,me) + + broken = broken .or. any(IEEE_is_NaN(sourceState(ph)%p(i)%dotState(:,me))) + + enddo SourceLoop + +end function constitutive_thermal_collectDotState !-------------------------------------------------------------------------------------------------- !> @brief for constitutive models having an instantaneous change of state !> will return false if delta state is not needed/supported by the constitutive model !-------------------------------------------------------------------------------------------------- -function constitutive_deltaState_source(Fe, ipc, ip, el, phase, of) result(broken) +function constitutive_damage_deltaState(Fe, co, ip, el, ph, of) result(broken) integer, intent(in) :: & - ipc, & !< component-ID of integration point + co, & !< component-ID of integration point ip, & !< integration point el, & !< element - phase, & + ph, & of real(pReal), intent(in), dimension(3,3) :: & Fe !< elastic deformation gradient integer :: & - i, & + so, & myOffset, & mySize logical :: & @@ -738,26 +703,26 @@ function constitutive_deltaState_source(Fe, ipc, ip, el, phase, of) result(broke broken = .false. - sourceLoop: do i = 1, phase_Nsources(phase) + sourceLoop: do so = 1, phase_Nsources(ph) - sourceType: select case (phase_source(i,phase)) + sourceType: select case (phase_source(so,ph)) case (SOURCE_damage_isoBrittle_ID) sourceType - call source_damage_isoBrittle_deltaState (constitutive_homogenizedC(ipc,ip,el), Fe, & - ipc, ip, el) - broken = any(IEEE_is_NaN(sourceState(phase)%p(i)%deltaState(:,of))) + call source_damage_isoBrittle_deltaState (constitutive_homogenizedC(co,ip,el), Fe, & + co, ip, el) + broken = any(IEEE_is_NaN(sourceState(ph)%p(so)%deltaState(:,of))) if(.not. broken) then - myOffset = sourceState(phase)%p(i)%offsetDeltaState - mySize = sourceState(phase)%p(i)%sizeDeltaState - sourceState(phase)%p(i)%state(myOffset + 1: myOffset + mySize,of) = & - sourceState(phase)%p(i)%state(myOffset + 1: myOffset + mySize,of) + sourceState(phase)%p(i)%deltaState(1:mySize,of) + myOffset = sourceState(ph)%p(so)%offsetDeltaState + mySize = sourceState(ph)%p(so)%sizeDeltaState + sourceState(ph)%p(so)%state(myOffset + 1: myOffset + mySize,of) = & + sourceState(ph)%p(so)%state(myOffset + 1: myOffset + mySize,of) + sourceState(ph)%p(so)%deltaState(1:mySize,of) endif end select sourceType enddo SourceLoop -end function constitutive_deltaState_source +end function constitutive_damage_deltaState !-------------------------------------------------------------------------------------------------- @@ -774,20 +739,21 @@ subroutine constitutive_allocateState(state, & sizeDotState, & sizeDeltaState + state%sizeState = sizeState state%sizeDotState = sizeDotState state%sizeDeltaState = sizeDeltaState state%offsetDeltaState = sizeState-sizeDeltaState ! deltaState occupies latter part of state by definition - allocate(state%atol (sizeState), source=0.0_pReal) - allocate(state%state0 (sizeState,Nconstituents), source=0.0_pReal) + allocate(state%atol (sizeState), source=0.0_pReal) + allocate(state%state0 (sizeState,Nconstituents), source=0.0_pReal) allocate(state%partitionedState0(sizeState,Nconstituents), source=0.0_pReal) - allocate(state%subState0 (sizeState,Nconstituents), source=0.0_pReal) - allocate(state%state (sizeState,Nconstituents), source=0.0_pReal) + allocate(state%subState0 (sizeState,Nconstituents), source=0.0_pReal) + allocate(state%state (sizeState,Nconstituents), source=0.0_pReal) - allocate(state%dotState (sizeDotState,Nconstituents), source=0.0_pReal) + allocate(state%dotState (sizeDotState,Nconstituents), source=0.0_pReal) - allocate(state%deltaState(sizeDeltaState,Nconstituents), source=0.0_pReal) + allocate(state%deltaState (sizeDeltaState,Nconstituents), source=0.0_pReal) end subroutine constitutive_allocateState @@ -796,22 +762,27 @@ end subroutine constitutive_allocateState !-------------------------------------------------------------------------------------------------- !> @brief Restore data after homog cutback. !-------------------------------------------------------------------------------------------------- -subroutine constitutive_restore(i,e) +subroutine constitutive_restore(ip,el,includeL) + logical, intent(in) :: includeL integer, intent(in) :: & - i, & !< integration point number - e !< element number - integer :: & - c, & !< constituent number - s + ip, & !< integration point number + el !< element number - do c = 1,homogenization_Nconstituents(material_homogenizationAt(e)) - do s = 1, phase_Nsources(material_phaseAt(c,e)) - sourceState(material_phaseAt(c,e))%p(s)%state( :,material_phasememberAt(c,i,e)) = & - sourceState(material_phaseAt(c,e))%p(s)%partitionedState0(:,material_phasememberAt(c,i,e)) + integer :: & + co, & !< constituent number + so + + + do co = 1,homogenization_Nconstituents(material_homogenizationAt(el)) + do so = 1, phase_Nsources(material_phaseAt(co,el)) + sourceState(material_phaseAt(co,el))%p(so)%state( :,material_phasememberAt(co,ip,el)) = & + sourceState(material_phaseAt(co,el))%p(so)%partitionedState0(:,material_phasememberAt(co,ip,el)) enddo enddo + call mech_restore(ip,el,includeL) + end subroutine constitutive_restore @@ -823,6 +794,12 @@ subroutine constitutive_forward integer :: i, j + crystallite_F0 = crystallite_F + crystallite_Lp0 = crystallite_Lp + crystallite_S0 = crystallite_S + + call constitutive_mech_forward() + do i = 1, size(sourceState) do j = 1,phase_Nsources(i) sourceState(i)%p(j)%state0 = sourceState(i)%p(j)%state @@ -836,8 +813,21 @@ end subroutine constitutive_forward !-------------------------------------------------------------------------------------------------- subroutine constitutive_results - call plastic_results - call damage_results + integer :: ph + character(len=:), allocatable :: group + + + call results_closeGroup(results_addGroup('/current/phase/')) + + do ph = 1, size(material_name_phase) + + group = '/current/phase/'//trim(material_name_phase(ph))//'/' + call results_closeGroup(results_addGroup(group)) + + call mech_results(group,ph) + call damage_results(group,ph) + + enddo end subroutine constitutive_results @@ -849,15 +839,15 @@ subroutine crystallite_init integer :: & Nconstituents, & - p, & - m, & - c, & !< counter in integration point component loop - i, & !< counter in integration point loop - e, & !< counter in element loop + ph, & + me, & + co, & !< counter in integration point component loop + ip, & !< counter in integration point loop + el, & !< counter in element loop cMax, & !< maximum number of integration point components iMax, & !< maximum number of integration points eMax !< maximum number of elements - + class(tNode), pointer :: & num_crystallite, & @@ -866,6 +856,7 @@ subroutine crystallite_init phase, & mech + print'(/,a)', ' <<<+- crystallite init -+>>>' debug_crystallite => config_debug%get('crystallite', defaultVal=emptyList) @@ -875,28 +866,21 @@ subroutine crystallite_init iMax = discretization_nIPs eMax = discretization_Nelems - allocate(crystallite_partitionedF(3,3,cMax,iMax,eMax),source=0.0_pReal) + allocate(crystallite_F(3,3,cMax,iMax,eMax),source=0.0_pReal) allocate(crystallite_S0, & - crystallite_F0,crystallite_Fp0,crystallite_Lp0, & + crystallite_F0,crystallite_Lp0, & crystallite_partitionedS0, & - crystallite_partitionedF0,crystallite_partitionedFp0,& - crystallite_partitionedLp0, & + crystallite_partitionedF0,& + crystallite_partitionedLp0, & crystallite_S,crystallite_P, & - crystallite_Fe,crystallite_Fp,crystallite_Lp, & - crystallite_subF,crystallite_subF0, & + crystallite_Fe,crystallite_Lp, & crystallite_subFp0,crystallite_subFi0, & - source = crystallite_partitionedF) - - allocate(crystallite_dt(cMax,iMax,eMax),source=0.0_pReal) - allocate(crystallite_subdt,crystallite_subStep, & - source = crystallite_dt) + source = crystallite_F) + allocate(crystallite_subdt(cMax,iMax,eMax),source=0.0_pReal) allocate(crystallite_orientation(cMax,iMax,eMax)) - allocate(crystallite_requested(cMax,iMax,eMax), source=.false.) - allocate(crystallite_converged(cMax,iMax,eMax), source=.true.) - num_crystallite => config_numerics%get('crystallite',defaultVal=emptyDict) num%subStepMinCryst = num_crystallite%get_asFloat ('subStepMin', defaultVal=1.0e-3_pReal) @@ -927,45 +911,30 @@ subroutine crystallite_init if(num%nState < 1) call IO_error(301,ext_msg='nState') if(num%nStress< 1) call IO_error(301,ext_msg='nStress') - select case(num_crystallite%get_asString('integrator',defaultVal='FPI')) - case('FPI') - integrateState => integrateStateFPI - case('Euler') - integrateState => integrateStateEuler - case('AdaptiveEuler') - integrateState => integrateStateAdaptiveEuler - case('RK4') - integrateState => integrateStateRK4 - case('RKCK45') - integrateState => integrateStateRKCK45 - case default - call IO_error(301,ext_msg='integrator') - end select phases => config_material%get('phase') - allocate(output_constituent(phases%length)) allocate(constitutive_mech_Fi(phases%length)) allocate(constitutive_mech_Fi0(phases%length)) - allocate(constitutive_mech_partionedFi0(phases%length)) + allocate(constitutive_mech_partitionedFi0(phases%length)) + allocate(constitutive_mech_Fp(phases%length)) + allocate(constitutive_mech_Fp0(phases%length)) + allocate(constitutive_mech_partitionedFp0(phases%length)) allocate(constitutive_mech_Li(phases%length)) allocate(constitutive_mech_Li0(phases%length)) - allocate(constitutive_mech_partionedLi0(phases%length)) - do p = 1, phases%length - Nconstituents = count(material_phaseAt == p) * discretization_nIPs - phase => phases%get(p) - mech => phase%get('mechanics',defaultVal = emptyDict) -#if defined(__GFORTRAN__) - output_constituent(p)%label = output_asStrings(mech) -#else - output_constituent(p)%label = mech%get_asStrings('output',defaultVal=emptyStringArray) -#endif - allocate(constitutive_mech_Fi(p)%data(3,3,Nconstituents)) - allocate(constitutive_mech_Fi0(p)%data(3,3,Nconstituents)) - allocate(constitutive_mech_partionedFi0(p)%data(3,3,Nconstituents)) - allocate(constitutive_mech_Li(p)%data(3,3,Nconstituents)) - allocate(constitutive_mech_Li0(p)%data(3,3,Nconstituents)) - allocate(constitutive_mech_partionedLi0(p)%data(3,3,Nconstituents)) + allocate(constitutive_mech_partitionedLi0(phases%length)) + do ph = 1, phases%length + Nconstituents = count(material_phaseAt == ph) * discretization_nIPs + + allocate(constitutive_mech_Fi(ph)%data(3,3,Nconstituents)) + allocate(constitutive_mech_Fi0(ph)%data(3,3,Nconstituents)) + allocate(constitutive_mech_partitionedFi0(ph)%data(3,3,Nconstituents)) + allocate(constitutive_mech_Fp(ph)%data(3,3,Nconstituents)) + allocate(constitutive_mech_Fp0(ph)%data(3,3,Nconstituents)) + allocate(constitutive_mech_partitionedFp0(ph)%data(3,3,Nconstituents)) + allocate(constitutive_mech_Li(ph)%data(3,3,Nconstituents)) + allocate(constitutive_mech_Li0(ph)%data(3,3,Nconstituents)) + allocate(constitutive_mech_partitionedLi0(ph)%data(3,3,Nconstituents)) enddo print'(a42,1x,i10)', ' # of elements: ', eMax @@ -973,46 +942,43 @@ subroutine crystallite_init print'(a42,1x,i10)', 'max # of constituents/integration point: ', cMax flush(IO_STDOUT) - !$OMP PARALLEL DO PRIVATE(p,m) - do e = FEsolving_execElem(1),FEsolving_execElem(2) - do i = FEsolving_execIP(1), FEsolving_execIP(2); do c = 1, homogenization_Nconstituents(material_homogenizationAt(e)) + !$OMP PARALLEL DO PRIVATE(ph,me) + do el = 1, size(material_phaseMemberAt,3); do ip = 1, size(material_phaseMemberAt,2) + do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) - p = material_phaseAt(c,e) - m = material_phaseMemberAt(c,i,e) - crystallite_Fp0(1:3,1:3,c,i,e) = material_orientation0(c,i,e)%asMatrix() ! Fp reflects initial orientation (see 10.1016/j.actamat.2006.01.005) - crystallite_Fp0(1:3,1:3,c,i,e) = crystallite_Fp0(1:3,1:3,c,i,e) & - / math_det33(crystallite_Fp0(1:3,1:3,c,i,e))**(1.0_pReal/3.0_pReal) - constitutive_mech_Fi0(p)%data(1:3,1:3,m) = math_I3 + ph = material_phaseAt(co,el) + me = material_phaseMemberAt(co,ip,el) + constitutive_mech_Fp0(ph)%data(1:3,1:3,me) = material_orientation0(co,ip,el)%asMatrix() ! Fp reflects initial orientation (see 10.1016/j.actamat.2006.01.005) + constitutive_mech_Fp0(ph)%data(1:3,1:3,me) = constitutive_mech_Fp0(ph)%data(1:3,1:3,me) & + / math_det33(constitutive_mech_Fp0(ph)%data(1:3,1:3,me))**(1.0_pReal/3.0_pReal) + constitutive_mech_Fi0(ph)%data(1:3,1:3,me) = math_I3 - crystallite_F0(1:3,1:3,c,i,e) = math_I3 + crystallite_F0(1:3,1:3,co,ip,el) = math_I3 - crystallite_Fe(1:3,1:3,c,i,e) = math_inv33(matmul(constitutive_mech_Fi0(p)%data(1:3,1:3,m), & - crystallite_Fp0(1:3,1:3,c,i,e))) ! assuming that euler angles are given in internal strain free configuration - crystallite_Fp(1:3,1:3,c,i,e) = crystallite_Fp0(1:3,1:3,c,i,e) - constitutive_mech_Fi(p)%data(1:3,1:3,m) = constitutive_mech_Fi0(p)%data(1:3,1:3,m) + crystallite_Fe(1:3,1:3,co,ip,el) = math_inv33(matmul(constitutive_mech_Fi0(ph)%data(1:3,1:3,me), & + constitutive_mech_Fp0(ph)%data(1:3,1:3,me))) ! assuming that euler angles are given in internal strain free configuration + constitutive_mech_Fp(ph)%data(1:3,1:3,me) = constitutive_mech_Fp0(ph)%data(1:3,1:3,me) + constitutive_mech_Fi(ph)%data(1:3,1:3,me) = constitutive_mech_Fi0(ph)%data(1:3,1:3,me) - constitutive_mech_partionedFi0(p)%data(1:3,1:3,m) = constitutive_mech_Fi0(p)%data(1:3,1:3,m) + constitutive_mech_partitionedFi0(ph)%data(1:3,1:3,me) = constitutive_mech_Fi0(ph)%data(1:3,1:3,me) + constitutive_mech_partitionedFp0(ph)%data(1:3,1:3,me) = constitutive_mech_Fp0(ph)%data(1:3,1:3,me) - - crystallite_requested(c,i,e) = .true. - enddo; enddo - enddo + enddo + enddo; enddo !$OMP END PARALLEL DO - - crystallite_partitionedFp0 = crystallite_Fp0 crystallite_partitionedF0 = crystallite_F0 - crystallite_partitionedF = crystallite_F0 + crystallite_F = crystallite_F0 - call crystallite_orientations() - !$OMP PARALLEL DO - do e = FEsolving_execElem(1),FEsolving_execElem(2) - do i = FEsolving_execIP(1),FEsolving_execIP(2) - do c = 1,homogenization_Nconstituents(material_homogenizationAt(e)) - call constitutive_plastic_dependentState(crystallite_partitionedF0(1:3,1:3,c,i,e), & - crystallite_partitionedFp0(1:3,1:3,c,i,e), & - c,i,e) ! update dependent state variables to be consistent with basic states + !$OMP PARALLEL DO PRIVATE(ph,me) + do el = 1, size(material_phaseMemberAt,3) + do ip = 1, size(material_phaseMemberAt,2) + do co = 1,homogenization_Nconstituents(material_homogenizationAt(el)) + ph = material_phaseAt(co,el) + me = material_phaseMemberAt(co,ip,el) + call crystallite_orientations(co,ip,el) + call constitutive_plastic_dependentState(crystallite_partitionedF0(1:3,1:3,co,ip,el),co,ip,el) ! update dependent state variables to be consistent with basic states enddo enddo enddo @@ -1022,271 +988,83 @@ subroutine crystallite_init end subroutine crystallite_init -!-------------------------------------------------------------------------------------------------- -!> @brief calculate stress (P) -!-------------------------------------------------------------------------------------------------- -function crystallite_stress() - - logical, dimension(discretization_nIPs,discretization_Nelems) :: crystallite_stress - real(pReal) :: & - formerSubStep - integer :: & - NiterationCrystallite, & ! number of iterations in crystallite loop - c, & !< counter in integration point component loop - i, & !< counter in integration point loop - e, & !< counter in element loop - s, p, m - logical, dimension(homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems) :: todo !ToDo: need to set some values to false for different Ngrains - real(pReal), dimension(homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems) :: subFrac !ToDo: need to set some values to false for different Ngrains - real(pReal), dimension(:,:,:,:,:), allocatable :: & - subLp0,& !< plastic velocity grad at start of crystallite inc - subLi0 !< intermediate velocity grad at start of crystallite inc - - todo = .false. - - allocate(subLi0(3,3,homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems)) - subLp0 = crystallite_partitionedLp0 - -!-------------------------------------------------------------------------------------------------- -! initialize to starting condition - crystallite_subStep = 0.0_pReal - !$OMP PARALLEL DO PRIVATE(p,m) - elementLooping1: do e = FEsolving_execElem(1),FEsolving_execElem(2) - do i = FEsolving_execIP(1),FEsolving_execIP(2); do c = 1,homogenization_Nconstituents(material_homogenizationAt(e)) - p = material_phaseAt(c,e) - m = material_phaseMemberAt(c,i,e) - subLi0(1:3,1:3,c,i,e) = constitutive_mech_partionedLi0(p)%data(1:3,1:3,m) - homogenizationRequestsCalculation: if (crystallite_requested(c,i,e)) then - plasticState (material_phaseAt(c,e))%subState0( :,material_phaseMemberAt(c,i,e)) = & - plasticState (material_phaseAt(c,e))%partitionedState0(:,material_phaseMemberAt(c,i,e)) - - do s = 1, phase_Nsources(material_phaseAt(c,e)) - sourceState(material_phaseAt(c,e))%p(s)%subState0( :,material_phaseMemberAt(c,i,e)) = & - sourceState(material_phaseAt(c,e))%p(s)%partitionedState0(:,material_phaseMemberAt(c,i,e)) - enddo - crystallite_subFp0(1:3,1:3,c,i,e) = crystallite_partitionedFp0(1:3,1:3,c,i,e) - crystallite_subFi0(1:3,1:3,c,i,e) = constitutive_mech_partionedFi0(p)%data(1:3,1:3,m) - crystallite_subF0(1:3,1:3,c,i,e) = crystallite_partitionedF0(1:3,1:3,c,i,e) - subFrac(c,i,e) = 0.0_pReal - crystallite_subStep(c,i,e) = 1.0_pReal/num%subStepSizeCryst - todo(c,i,e) = .true. - crystallite_converged(c,i,e) = .false. ! pretend failed step of 1/subStepSizeCryst - endif homogenizationRequestsCalculation - enddo; enddo - enddo elementLooping1 - !$OMP END PARALLEL DO - - NiterationCrystallite = 0 - cutbackLooping: do while (any(todo(:,FEsolving_execIP(1):FEsolving_execIP(2),FEsolving_execELem(1):FEsolving_execElem(2)))) - NiterationCrystallite = NiterationCrystallite + 1 - -#ifdef DEBUG - if (debugCrystallite%extensive) & - print'(a,i6)', '<< CRYST stress >> crystallite iteration ',NiterationCrystallite -#endif - !$OMP PARALLEL DO PRIVATE(formerSubStep,p,m) - elementLooping3: do e = FEsolving_execElem(1),FEsolving_execElem(2) - do i = FEsolving_execIP(1),FEsolving_execIP(2) - do c = 1,homogenization_Nconstituents(material_homogenizationAt(e)) - p = material_phaseAt(c,e) - m = material_phaseMemberAt(c,i,e) -!-------------------------------------------------------------------------------------------------- -! wind forward - if (crystallite_converged(c,i,e)) then - formerSubStep = crystallite_subStep(c,i,e) - subFrac(c,i,e) = subFrac(c,i,e) + crystallite_subStep(c,i,e) - crystallite_subStep(c,i,e) = min(1.0_pReal - subFrac(c,i,e), & - num%stepIncreaseCryst * crystallite_subStep(c,i,e)) - - todo(c,i,e) = crystallite_subStep(c,i,e) > 0.0_pReal ! still time left to integrate on? - - if (todo(c,i,e)) then - crystallite_subF0 (1:3,1:3,c,i,e) = crystallite_subF(1:3,1:3,c,i,e) - subLp0(1:3,1:3,c,i,e) = crystallite_Lp (1:3,1:3,c,i,e) - subLi0(1:3,1:3,c,i,e) = constitutive_mech_Li(p)%data(1:3,1:3,m) - crystallite_subFp0(1:3,1:3,c,i,e) = crystallite_Fp (1:3,1:3,c,i,e) - crystallite_subFi0(1:3,1:3,c,i,e) = constitutive_mech_Fi(p)%data(1:3,1:3,m) - plasticState( material_phaseAt(c,e))%subState0(:,material_phaseMemberAt(c,i,e)) & - = plasticState(material_phaseAt(c,e))%state( :,material_phaseMemberAt(c,i,e)) - do s = 1, phase_Nsources(material_phaseAt(c,e)) - sourceState( material_phaseAt(c,e))%p(s)%subState0(:,material_phaseMemberAt(c,i,e)) & - = sourceState(material_phaseAt(c,e))%p(s)%state( :,material_phaseMemberAt(c,i,e)) - enddo - endif - -!-------------------------------------------------------------------------------------------------- -! cut back (reduced time and restore) - else - crystallite_subStep(c,i,e) = num%subStepSizeCryst * crystallite_subStep(c,i,e) - crystallite_Fp (1:3,1:3,c,i,e) = crystallite_subFp0(1:3,1:3,c,i,e) - constitutive_mech_Fi(p)%data(1:3,1:3,m) = crystallite_subFi0(1:3,1:3,c,i,e) - crystallite_S (1:3,1:3,c,i,e) = crystallite_S0 (1:3,1:3,c,i,e) - if (crystallite_subStep(c,i,e) < 1.0_pReal) then ! actual (not initial) cutback - crystallite_Lp (1:3,1:3,c,i,e) = subLp0(1:3,1:3,c,i,e) - constitutive_mech_Li(p)%data(1:3,1:3,m) = subLi0(1:3,1:3,c,i,e) - endif - plasticState (material_phaseAt(c,e))%state( :,material_phaseMemberAt(c,i,e)) & - = plasticState(material_phaseAt(c,e))%subState0(:,material_phaseMemberAt(c,i,e)) - do s = 1, phase_Nsources(material_phaseAt(c,e)) - sourceState( material_phaseAt(c,e))%p(s)%state( :,material_phaseMemberAt(c,i,e)) & - = sourceState(material_phaseAt(c,e))%p(s)%subState0(:,material_phaseMemberAt(c,i,e)) - enddo - - ! cant restore dotState here, since not yet calculated in first cutback after initialization - todo(c,i,e) = crystallite_subStep(c,i,e) > num%subStepMinCryst ! still on track or already done (beyond repair) - endif - -!-------------------------------------------------------------------------------------------------- -! prepare for integration - if (todo(c,i,e)) then - crystallite_subF(1:3,1:3,c,i,e) = crystallite_subF0(1:3,1:3,c,i,e) & - + crystallite_subStep(c,i,e) *( crystallite_partitionedF (1:3,1:3,c,i,e) & - -crystallite_partitionedF0(1:3,1:3,c,i,e)) - crystallite_Fe(1:3,1:3,c,i,e) = matmul(crystallite_subF(1:3,1:3,c,i,e), & - math_inv33(matmul(constitutive_mech_Fi(p)%data(1:3,1:3,m), & - crystallite_Fp(1:3,1:3,c,i,e)))) - crystallite_subdt(c,i,e) = crystallite_subStep(c,i,e) * crystallite_dt(c,i,e) - crystallite_converged(c,i,e) = .false. - call integrateState(c,i,e) - call integrateSourceState(c,i,e) - endif - - enddo - enddo - enddo elementLooping3 - !$OMP END PARALLEL DO - -!-------------------------------------------------------------------------------------------------- -! integrate --- requires fully defined state array (basic + dependent state) - where(.not. crystallite_converged .and. crystallite_subStep > num%subStepMinCryst) & ! do not try non-converged but fully cutbacked any further - todo = .true. ! TODO: again unroll this into proper elementloop to avoid N^2 for single point evaluation - enddo cutbackLooping - -! return whether converged or not - crystallite_stress = .false. - elementLooping5: do e = FEsolving_execElem(1),FEsolving_execElem(2) - do i = FEsolving_execIP(1),FEsolving_execIP(2) - crystallite_stress(i,e) = all(crystallite_converged(:,i,e)) - enddo - enddo elementLooping5 - -end function crystallite_stress - - !-------------------------------------------------------------------------------------------------- !> @brief Backup data for homog cutback. !-------------------------------------------------------------------------------------------------- -subroutine crystallite_initializeRestorationPoints(i,e) +subroutine constitutive_initializeRestorationPoints(ip,el) integer, intent(in) :: & - i, & !< integration point number - e !< element number + ip, & !< integration point number + el !< element number integer :: & - c, & !< constituent number - s,p, m + co, & !< constituent number + so,ph, me - do c = 1,homogenization_Nconstituents(material_homogenizationAt(e)) - p = material_phaseAt(c,e) - m = material_phaseMemberAt(c,i,e) - crystallite_partitionedFp0(1:3,1:3,c,i,e) = crystallite_Fp0(1:3,1:3,c,i,e) - crystallite_partitionedLp0(1:3,1:3,c,i,e) = crystallite_Lp0(1:3,1:3,c,i,e) - constitutive_mech_partionedFi0(p)%data(1:3,1:3,m) = constitutive_mech_Fi0(p)%data(1:3,1:3,m) - constitutive_mech_partionedLi0(p)%data(1:3,1:3,m) = constitutive_mech_Li0(p)%data(1:3,1:3,m) - crystallite_partitionedF0(1:3,1:3,c,i,e) = crystallite_F0(1:3,1:3,c,i,e) - crystallite_partitionedS0(1:3,1:3,c,i,e) = crystallite_S0(1:3,1:3,c,i,e) + do co = 1,homogenization_Nconstituents(material_homogenizationAt(el)) + ph = material_phaseAt(co,el) + me = material_phaseMemberAt(co,ip,el) + crystallite_partitionedLp0(1:3,1:3,co,ip,el) = crystallite_Lp0(1:3,1:3,co,ip,el) + crystallite_partitionedF0(1:3,1:3,co,ip,el) = crystallite_F0(1:3,1:3,co,ip,el) + crystallite_partitionedS0(1:3,1:3,co,ip,el) = crystallite_S0(1:3,1:3,co,ip,el) - plasticState(material_phaseAt(c,e))%partitionedState0(:,material_phasememberAt(c,i,e)) = & - plasticState(material_phaseAt(c,e))%state0( :,material_phasememberAt(c,i,e)) - do s = 1, phase_Nsources(material_phaseAt(c,e)) - sourceState(material_phaseAt(c,e))%p(s)%partitionedState0(:,material_phasememberAt(c,i,e)) = & - sourceState(material_phaseAt(c,e))%p(s)%state0( :,material_phasememberAt(c,i,e)) + call mech_initializeRestorationPoints(ph,me) + + do so = 1, phase_Nsources(material_phaseAt(co,el)) + sourceState(material_phaseAt(co,el))%p(so)%partitionedState0(:,material_phasememberAt(co,ip,el)) = & + sourceState(material_phaseAt(co,el))%p(so)%state0( :,material_phasememberAt(co,ip,el)) enddo enddo -end subroutine crystallite_initializeRestorationPoints +end subroutine constitutive_initializeRestorationPoints !-------------------------------------------------------------------------------------------------- !> @brief Wind homog inc forward. !-------------------------------------------------------------------------------------------------- -subroutine crystallite_windForward(i,e) +subroutine constitutive_windForward(ip,el) integer, intent(in) :: & - i, & !< integration point number - e !< element number - integer :: & - c, & !< constituent number - s, p, m - do c = 1,homogenization_Nconstituents(material_homogenizationAt(e)) - p = material_phaseAt(c,e) - m = material_phaseMemberAt(c,i,e) - crystallite_partitionedF0 (1:3,1:3,c,i,e) = crystallite_partitionedF(1:3,1:3,c,i,e) - crystallite_partitionedFp0(1:3,1:3,c,i,e) = crystallite_Fp (1:3,1:3,c,i,e) - crystallite_partitionedLp0(1:3,1:3,c,i,e) = crystallite_Lp (1:3,1:3,c,i,e) - constitutive_mech_partionedFi0(p)%data(1:3,1:3,m) = constitutive_mech_Fi(p)%data(1:3,1:3,m) - constitutive_mech_partionedLi0(p)%data(1:3,1:3,m) = constitutive_mech_Li(p)%data(1:3,1:3,m) - crystallite_partitionedS0 (1:3,1:3,c,i,e) = crystallite_S (1:3,1:3,c,i,e) + ip, & !< integration point number + el !< element number - plasticState (material_phaseAt(c,e))%partitionedState0(:,material_phasememberAt(c,i,e)) = & - plasticState (material_phaseAt(c,e))%state (:,material_phasememberAt(c,i,e)) - do s = 1, phase_Nsources(material_phaseAt(c,e)) - sourceState(material_phaseAt(c,e))%p(s)%partitionedState0(:,material_phasememberAt(c,i,e)) = & - sourceState(material_phaseAt(c,e))%p(s)%state (:,material_phasememberAt(c,i,e)) + integer :: & + co, & !< constituent number + so, ph, me + + + do co = 1,homogenization_Nconstituents(material_homogenizationAt(el)) + ph = material_phaseAt(co,el) + me = material_phaseMemberAt(co,ip,el) + crystallite_partitionedF0 (1:3,1:3,co,ip,el) = crystallite_F (1:3,1:3,co,ip,el) + crystallite_partitionedLp0(1:3,1:3,co,ip,el) = crystallite_Lp(1:3,1:3,co,ip,el) + crystallite_partitionedS0 (1:3,1:3,co,ip,el) = crystallite_S (1:3,1:3,co,ip,el) + + call constitutive_mech_windForward(ph,me) + do so = 1, phase_Nsources(material_phaseAt(co,el)) + sourceState(ph)%p(so)%partitionedState0(:,me) = sourceState(ph)%p(so)%state(:,me) enddo enddo -end subroutine crystallite_windForward - - -!-------------------------------------------------------------------------------------------------- -!> @brief Restore data after homog cutback. -!-------------------------------------------------------------------------------------------------- -subroutine crystallite_restore(i,e,includeL) - - integer, intent(in) :: & - i, & !< integration point number - e !< element number - logical, intent(in) :: & - includeL !< protect agains fake cutback - integer :: & - c, p, m !< constituent number - - do c = 1,homogenization_Nconstituents(material_homogenizationAt(e)) - p = material_phaseAt(c,e) - m = material_phaseMemberAt(c,i,e) - if (includeL) then - crystallite_Lp(1:3,1:3,c,i,e) = crystallite_partitionedLp0(1:3,1:3,c,i,e) - constitutive_mech_Li(p)%data(1:3,1:3,m) = constitutive_mech_partionedLi0(p)%data(1:3,1:3,m) - endif ! maybe protecting everything from overwriting makes more sense - - crystallite_Fp(1:3,1:3,c,i,e) = crystallite_partitionedFp0(1:3,1:3,c,i,e) - constitutive_mech_Fi(p)%data(1:3,1:3,m) = constitutive_mech_partionedFi0(p)%data(1:3,1:3,m) - crystallite_S (1:3,1:3,c,i,e) = crystallite_partitionedS0 (1:3,1:3,c,i,e) - - plasticState (material_phaseAt(c,e))%state( :,material_phasememberAt(c,i,e)) = & - plasticState (material_phaseAt(c,e))%partitionedState0(:,material_phasememberAt(c,i,e)) - enddo - -end subroutine crystallite_restore +end subroutine constitutive_windForward !-------------------------------------------------------------------------------------------------- !> @brief Calculate tangent (dPdF). !-------------------------------------------------------------------------------------------------- -function crystallite_stressTangent(c,i,e) result(dPdF) +function crystallite_stressTangent(co,ip,el) result(dPdF) real(pReal), dimension(3,3,3,3) :: dPdF integer, intent(in) :: & - c, & !< counter in constituent loop - i, & !< counter in integration point loop - e !< counter in element loop + co, & !< counter in constituent loop + ip, & !< counter in integration point loop + el !< counter in element loop + integer :: & o, & - p, pp, m - + p, ph, me real(pReal), dimension(3,3) :: devNull, & invSubFp0,invSubFi0,invFp,invFi, & - temp_33_1, temp_33_2, temp_33_3, temp_33_4 + temp_33_1, temp_33_2, temp_33_3 real(pReal), dimension(3,3,3,3) :: dSdFe, & dSdF, & dSdFi, & @@ -1302,100 +1080,100 @@ function crystallite_stressTangent(c,i,e) result(dPdF) real(pReal), dimension(9,9):: temp_99 logical :: error - pp = material_phaseAt(c,e) - m = material_phaseMemberAt(c,i,e) - call constitutive_hooke_SandItsTangents(devNull,dSdFe,dSdFi, & - crystallite_Fe(1:3,1:3,c,i,e), & - constitutive_mech_Fi(pp)%data(1:3,1:3,m),c,i,e) - call constitutive_LiAndItsTangents(devNull,dLidS,dLidFi, & - crystallite_S (1:3,1:3,c,i,e), & - constitutive_mech_Fi(pp)%data(1:3,1:3,m), & - c,i,e) + ph = material_phaseAt(co,el) + me = material_phaseMemberAt(co,ip,el) - invFp = math_inv33(crystallite_Fp(1:3,1:3,c,i,e)) - invFi = math_inv33(constitutive_mech_Fi(pp)%data(1:3,1:3,m)) - invSubFp0 = math_inv33(crystallite_subFp0(1:3,1:3,c,i,e)) - invSubFi0 = math_inv33(crystallite_subFi0(1:3,1:3,c,i,e)) + call constitutive_hooke_SandItsTangents(devNull,dSdFe,dSdFi, & + crystallite_Fe(1:3,1:3,co,ip,el), & + constitutive_mech_Fi(ph)%data(1:3,1:3,me),co,ip,el) + call constitutive_LiAndItsTangents(devNull,dLidS,dLidFi, & + crystallite_S (1:3,1:3,co,ip,el), & + constitutive_mech_Fi(ph)%data(1:3,1:3,me), & + co,ip,el) - if (sum(abs(dLidS)) < tol_math_check) then - dFidS = 0.0_pReal - else - lhs_3333 = 0.0_pReal; rhs_3333 = 0.0_pReal - do o=1,3; do p=1,3 - lhs_3333(1:3,1:3,o,p) = lhs_3333(1:3,1:3,o,p) & - + crystallite_subdt(c,i,e)*matmul(invSubFi0,dLidFi(1:3,1:3,o,p)) - lhs_3333(1:3,o,1:3,p) = lhs_3333(1:3,o,1:3,p) & - + invFi*invFi(p,o) - rhs_3333(1:3,1:3,o,p) = rhs_3333(1:3,1:3,o,p) & - - crystallite_subdt(c,i,e)*matmul(invSubFi0,dLidS(1:3,1:3,o,p)) - enddo; enddo - call math_invert(temp_99,error,math_3333to99(lhs_3333)) - if (error) then - call IO_warning(warning_ID=600,el=e,ip=i,g=c, & - ext_msg='inversion error in analytic tangent calculation') - dFidS = 0.0_pReal - else - dFidS = math_mul3333xx3333(math_99to3333(temp_99),rhs_3333) - endif - dLidS = math_mul3333xx3333(dLidFi,dFidS) + dLidS - endif + invFp = math_inv33(constitutive_mech_Fp(ph)%data(1:3,1:3,me)) + invFi = math_inv33(constitutive_mech_Fi(ph)%data(1:3,1:3,me)) + invSubFp0 = math_inv33(crystallite_subFp0(1:3,1:3,co,ip,el)) + invSubFi0 = math_inv33(crystallite_subFi0(1:3,1:3,co,ip,el)) - call constitutive_plastic_LpAndItsTangents(devNull,dLpdS,dLpdFi, & - crystallite_S (1:3,1:3,c,i,e), & - constitutive_mech_Fi(pp)%data(1:3,1:3,m),c,i,e) - dLpdS = math_mul3333xx3333(dLpdFi,dFidS) + dLpdS + if (sum(abs(dLidS)) < tol_math_check) then + dFidS = 0.0_pReal + else + lhs_3333 = 0.0_pReal; rhs_3333 = 0.0_pReal + do o=1,3; do p=1,3 + lhs_3333(1:3,1:3,o,p) = lhs_3333(1:3,1:3,o,p) & + + crystallite_subdt(co,ip,el)*matmul(invSubFi0,dLidFi(1:3,1:3,o,p)) + lhs_3333(1:3,o,1:3,p) = lhs_3333(1:3,o,1:3,p) & + + invFi*invFi(p,o) + rhs_3333(1:3,1:3,o,p) = rhs_3333(1:3,1:3,o,p) & + - crystallite_subdt(co,ip,el)*matmul(invSubFi0,dLidS(1:3,1:3,o,p)) + enddo; enddo + call math_invert(temp_99,error,math_3333to99(lhs_3333)) + if (error) then + call IO_warning(warning_ID=600,el=el,ip=ip,g=co, & + ext_msg='inversion error in analytic tangent calculation') + dFidS = 0.0_pReal + else + dFidS = math_mul3333xx3333(math_99to3333(temp_99),rhs_3333) + endif + dLidS = math_mul3333xx3333(dLidFi,dFidS) + dLidS + endif + + call constitutive_plastic_LpAndItsTangents(devNull,dLpdS,dLpdFi, & + crystallite_S (1:3,1:3,co,ip,el), & + constitutive_mech_Fi(ph)%data(1:3,1:3,me),co,ip,el) + dLpdS = math_mul3333xx3333(dLpdFi,dFidS) + dLpdS !-------------------------------------------------------------------------------------------------- ! calculate dSdF - temp_33_1 = transpose(matmul(invFp,invFi)) - temp_33_2 = matmul(crystallite_subF(1:3,1:3,c,i,e),invSubFp0) - temp_33_3 = matmul(matmul(crystallite_subF(1:3,1:3,c,i,e),invFp), invSubFi0) + temp_33_1 = transpose(matmul(invFp,invFi)) + temp_33_2 = matmul(crystallite_F(1:3,1:3,co,ip,el),invSubFp0) + temp_33_3 = matmul(matmul(crystallite_F(1:3,1:3,co,ip,el),invFp), invSubFi0) - do o=1,3; do p=1,3 - rhs_3333(p,o,1:3,1:3) = matmul(dSdFe(p,o,1:3,1:3),temp_33_1) - temp_3333(1:3,1:3,p,o) = matmul(matmul(temp_33_2,dLpdS(1:3,1:3,p,o)), invFi) & - + matmul(temp_33_3,dLidS(1:3,1:3,p,o)) - enddo; enddo - lhs_3333 = crystallite_subdt(c,i,e)*math_mul3333xx3333(dSdFe,temp_3333) & - + math_mul3333xx3333(dSdFi,dFidS) + do o=1,3; do p=1,3 + rhs_3333(p,o,1:3,1:3) = matmul(dSdFe(p,o,1:3,1:3),temp_33_1) + temp_3333(1:3,1:3,p,o) = matmul(matmul(temp_33_2,dLpdS(1:3,1:3,p,o)), invFi) & + + matmul(temp_33_3,dLidS(1:3,1:3,p,o)) + enddo; enddo + lhs_3333 = crystallite_subdt(co,ip,el)*math_mul3333xx3333(dSdFe,temp_3333) & + + math_mul3333xx3333(dSdFi,dFidS) - call math_invert(temp_99,error,math_eye(9)+math_3333to99(lhs_3333)) - if (error) then - call IO_warning(warning_ID=600,el=e,ip=i,g=c, & - ext_msg='inversion error in analytic tangent calculation') - dSdF = rhs_3333 - else - dSdF = math_mul3333xx3333(math_99to3333(temp_99),rhs_3333) - endif + call math_invert(temp_99,error,math_eye(9)+math_3333to99(lhs_3333)) + if (error) then + call IO_warning(warning_ID=600,el=el,ip=ip,g=co, & + ext_msg='inversion error in analytic tangent calculation') + dSdF = rhs_3333 + else + dSdF = math_mul3333xx3333(math_99to3333(temp_99),rhs_3333) + endif !-------------------------------------------------------------------------------------------------- ! calculate dFpinvdF - temp_3333 = math_mul3333xx3333(dLpdS,dSdF) - do o=1,3; do p=1,3 - dFpinvdF(1:3,1:3,p,o) = -crystallite_subdt(c,i,e) & - * matmul(invSubFp0, matmul(temp_3333(1:3,1:3,p,o),invFi)) - enddo; enddo + temp_3333 = math_mul3333xx3333(dLpdS,dSdF) + do o=1,3; do p=1,3 + dFpinvdF(1:3,1:3,p,o) = -crystallite_subdt(co,ip,el) & + * matmul(invSubFp0, matmul(temp_3333(1:3,1:3,p,o),invFi)) + enddo; enddo !-------------------------------------------------------------------------------------------------- ! assemble dPdF - temp_33_1 = matmul(crystallite_S(1:3,1:3,c,i,e),transpose(invFp)) - temp_33_2 = matmul(invFp,temp_33_1) - temp_33_3 = matmul(crystallite_subF(1:3,1:3,c,i,e),invFp) - temp_33_4 = matmul(temp_33_3,crystallite_S(1:3,1:3,c,i,e)) + temp_33_1 = matmul(crystallite_S(1:3,1:3,co,ip,el),transpose(invFp)) + temp_33_2 = matmul(crystallite_F(1:3,1:3,co,ip,el),invFp) + temp_33_3 = matmul(temp_33_2,crystallite_S(1:3,1:3,co,ip,el)) - dPdF = 0.0_pReal - do p=1,3 - dPdF(p,1:3,p,1:3) = transpose(temp_33_2) - enddo - do o=1,3; do p=1,3 - dPdF(1:3,1:3,p,o) = dPdF(1:3,1:3,p,o) & - + matmul(matmul(crystallite_subF(1:3,1:3,c,i,e), & - dFpinvdF(1:3,1:3,p,o)),temp_33_1) & - + matmul(matmul(temp_33_3,dSdF(1:3,1:3,p,o)), & - transpose(invFp)) & - + matmul(temp_33_4,transpose(dFpinvdF(1:3,1:3,p,o))) - enddo; enddo + dPdF = 0.0_pReal + do p=1,3 + dPdF(p,1:3,p,1:3) = transpose(matmul(invFp,temp_33_1)) + enddo + do o=1,3; do p=1,3 + dPdF(1:3,1:3,p,o) = dPdF(1:3,1:3,p,o) & + + matmul(matmul(crystallite_F(1:3,1:3,co,ip,el), & + dFpinvdF(1:3,1:3,p,o)),temp_33_1) & + + matmul(matmul(temp_33_2,dSdF(1:3,1:3,p,o)), & + transpose(invFp)) & + + matmul(temp_33_3,transpose(dFpinvdF(1:3,1:3,p,o))) + enddo; enddo end function crystallite_stressTangent @@ -1403,33 +1181,20 @@ end function crystallite_stressTangent !-------------------------------------------------------------------------------------------------- !> @brief calculates orientations !-------------------------------------------------------------------------------------------------- -subroutine crystallite_orientations +subroutine crystallite_orientations(co,ip,el) - integer & - c, & !< counter in integration point component loop - i, & !< counter in integration point loop - e !< counter in element loop + integer, intent(in) :: & + co, & !< counter in integration point component loop + ip, & !< counter in integration point loop + el !< counter in element loop - !$OMP PARALLEL DO - do e = FEsolving_execElem(1),FEsolving_execElem(2) - do i = FEsolving_execIP(1),FEsolving_execIP(2) - do c = 1,homogenization_Nconstituents(material_homogenizationAt(e)) - call crystallite_orientation(c,i,e)%fromMatrix(transpose(math_rotationalPart(crystallite_Fe(1:3,1:3,c,i,e)))) - enddo; enddo; enddo - !$OMP END PARALLEL DO - nonlocalPresent: if (any(plasticState%nonlocal)) then - !$OMP PARALLEL DO - do e = FEsolving_execElem(1),FEsolving_execElem(2) - if (plasticState(material_phaseAt(1,e))%nonlocal) then - do i = FEsolving_execIP(1),FEsolving_execIP(2) - call plastic_nonlocal_updateCompatibility(crystallite_orientation, & - phase_plasticityInstance(material_phaseAt(1,e)),i,e) - enddo - endif - enddo - !$OMP END PARALLEL DO - endif nonlocalPresent + call crystallite_orientation(co,ip,el)%fromMatrix(transpose(math_rotationalPart(crystallite_Fe(1:3,1:3,co,ip,el)))) + + if (plasticState(material_phaseAt(1,el))%nonlocal) & + call plastic_nonlocal_updateCompatibility(crystallite_orientation, & + phase_plasticityInstance(material_phaseAt(1,el)),ip,el) + end subroutine crystallite_orientations @@ -1437,225 +1202,104 @@ end subroutine crystallite_orientations !-------------------------------------------------------------------------------------------------- !> @brief Map 2nd order tensor to reference config !-------------------------------------------------------------------------------------------------- -function crystallite_push33ToRef(ipc,ip,el, tensor33) +function crystallite_push33ToRef(co,ip,el, tensor33) - real(pReal), dimension(3,3) :: crystallite_push33ToRef real(pReal), dimension(3,3), intent(in) :: tensor33 real(pReal), dimension(3,3) :: T integer, intent(in):: & el, & ip, & - ipc + co - T = matmul(material_orientation0(ipc,ip,el)%asMatrix(), & ! ToDo: initial orientation correct? - transpose(math_inv33(crystallite_subF(1:3,1:3,ipc,ip,el)))) + real(pReal), dimension(3,3) :: crystallite_push33ToRef + + + T = matmul(material_orientation0(co,ip,el)%asMatrix(), & ! ToDo: initial orientation correct? + transpose(math_inv33(crystallite_F(1:3,1:3,co,ip,el)))) crystallite_push33ToRef = matmul(transpose(T),matmul(tensor33,T)) end function crystallite_push33ToRef -!-------------------------------------------------------------------------------------------------- -!> @brief writes crystallite results to HDF5 output file -!-------------------------------------------------------------------------------------------------- -subroutine crystallite_results - - integer :: p,o - real(pReal), allocatable, dimension(:,:,:) :: selected_tensors - real(pReal), allocatable, dimension(:,:) :: selected_rotations - character(len=:), allocatable :: group,structureLabel - - do p=1,size(material_name_phase) - group = trim('current/phase')//'/'//trim(material_name_phase(p))//'/mechanics' - - call results_closeGroup(results_addGroup(group)) - - do o = 1, size(output_constituent(p)%label) - select case (output_constituent(p)%label(o)) - case('F') - selected_tensors = select_tensors(crystallite_partitionedF,p) - call results_writeDataset(group,selected_tensors,output_constituent(p)%label(o),& - 'deformation gradient','1') - case('F_e') - selected_tensors = select_tensors(crystallite_Fe,p) - call results_writeDataset(group,selected_tensors,output_constituent(p)%label(o),& - 'elastic deformation gradient','1') - case('F_p') - selected_tensors = select_tensors(crystallite_Fp,p) - call results_writeDataset(group,selected_tensors,output_constituent(p)%label(o),& - 'plastic deformation gradient','1') - case('F_i') - call results_writeDataset(group,constitutive_mech_Fi(p)%data,output_constituent(p)%label(o),& - 'inelastic deformation gradient','1') - case('L_p') - selected_tensors = select_tensors(crystallite_Lp,p) - call results_writeDataset(group,selected_tensors,output_constituent(p)%label(o),& - 'plastic velocity gradient','1/s') - case('L_i') - call results_writeDataset(group,constitutive_mech_Li(p)%data,output_constituent(p)%label(o),& - 'inelastic velocity gradient','1/s') - case('P') - selected_tensors = select_tensors(crystallite_P,p) - call results_writeDataset(group,selected_tensors,output_constituent(p)%label(o),& - 'First Piola-Kirchhoff stress','Pa') - case('S') - selected_tensors = select_tensors(crystallite_S,p) - call results_writeDataset(group,selected_tensors,output_constituent(p)%label(o),& - 'Second Piola-Kirchhoff stress','Pa') - case('O') - select case(lattice_structure(p)) - case(lattice_ISO_ID) - structureLabel = 'aP' - case(lattice_FCC_ID) - structureLabel = 'cF' - case(lattice_BCC_ID) - structureLabel = 'cI' - case(lattice_BCT_ID) - structureLabel = 'tI' - case(lattice_HEX_ID) - structureLabel = 'hP' - case(lattice_ORT_ID) - structureLabel = 'oP' - end select - selected_rotations = select_rotations(crystallite_orientation,p) - call results_writeDataset(group,selected_rotations,output_constituent(p)%label(o),& - 'crystal orientation as quaternion','q_0 ') - call results_addAttribute('Lattice',structureLabel,group//'/'//output_constituent(p)%label(o)) - end select - enddo - enddo - - contains - - !------------------------------------------------------------------------------------------------ - !> @brief select tensors for output - !------------------------------------------------------------------------------------------------ - function select_tensors(dataset,instance) - - integer, intent(in) :: instance - real(pReal), dimension(:,:,:,:,:), intent(in) :: dataset - real(pReal), allocatable, dimension(:,:,:) :: select_tensors - integer :: e,i,c,j - - allocate(select_tensors(3,3,count(material_phaseAt==instance)*discretization_nIPs)) - - j=0 - do e = 1, size(material_phaseAt,2) - do i = 1, discretization_nIPs - do c = 1, size(material_phaseAt,1) !ToDo: this needs to be changed for varying Ngrains - if (material_phaseAt(c,e) == instance) then - j = j + 1 - select_tensors(1:3,1:3,j) = dataset(1:3,1:3,c,i,e) - endif - enddo - enddo - enddo - - end function select_tensors - - -!-------------------------------------------------------------------------------------------------- -!> @brief select rotations for output -!-------------------------------------------------------------------------------------------------- - function select_rotations(dataset,instance) - - integer, intent(in) :: instance - type(rotation), dimension(:,:,:), intent(in) :: dataset - real(pReal), allocatable, dimension(:,:) :: select_rotations - integer :: e,i,c,j - - allocate(select_rotations(4,count(material_phaseAt==instance)*homogenization_maxNconstituents*discretization_nIPs)) - - j=0 - do e = 1, size(material_phaseAt,2) - do i = 1, discretization_nIPs - do c = 1, size(material_phaseAt,1) !ToDo: this needs to be changed for varying Ngrains - if (material_phaseAt(c,e) == instance) then - j = j + 1 - select_rotations(1:4,j) = dataset(c,i,e)%asQuaternion() - endif - enddo - enddo - enddo - - end function select_rotations - -end subroutine crystallite_results - - !-------------------------------------------------------------------------------------------------- !> @brief integrate stress, state with adaptive 1st order explicit Euler method !> using Fixed Point Iteration to adapt the stepsize !-------------------------------------------------------------------------------------------------- -subroutine integrateSourceState(g,i,e) +function integrateSourceState(dt,co,ip,el) result(broken) + real(pReal), intent(in) :: dt integer, intent(in) :: & - e, & !< element index in element loop - i, & !< integration point index in ip loop - g !< grain index in grain loop + el, & !< element index in element loop + ip, & !< integration point index in ip loop + co !< grain index in grain loop + integer :: & NiterationState, & !< number of iterations in state loop - p, & - c, & - s, & - size_pl + ph, & + me, & + so integer, dimension(maxval(phase_Nsources)) :: & size_so real(pReal) :: & zeta - real(pReal), dimension(max(constitutive_plasticity_maxSizeDotState,constitutive_source_maxSizeDotState)) :: & + real(pReal), dimension(constitutive_source_maxSizeDotState) :: & r ! state residuum real(pReal), dimension(constitutive_source_maxSizeDotState,2,maxval(phase_Nsources)) :: source_dotState logical :: & - broken + broken, converged_ - p = material_phaseAt(g,e) - c = material_phaseMemberAt(g,i,e) - broken = constitutive_collectDotState_source(crystallite_S(1:3,1:3,g,i,e), g,i,e,p,c) + ph = material_phaseAt(co,el) + me = material_phaseMemberAt(co,ip,el) + + converged_ = .true. + broken = constitutive_thermal_collectDotState(ph,me) + broken = broken .or. constitutive_damage_collectDotState(crystallite_S(1:3,1:3,co,ip,el), co,ip,el,ph,me) if(broken) return - do s = 1, phase_Nsources(p) - size_so(s) = sourceState(p)%p(s)%sizeDotState - sourceState(p)%p(s)%state(1:size_so(s),c) = sourceState(p)%p(s)%subState0(1:size_so(s),c) & - + sourceState(p)%p(s)%dotState (1:size_so(s),c) & - * crystallite_subdt(g,i,e) - source_dotState(1:size_so(s),2,s) = 0.0_pReal + do so = 1, phase_Nsources(ph) + size_so(so) = sourceState(ph)%p(so)%sizeDotState + sourceState(ph)%p(so)%state(1:size_so(so),me) = sourceState(ph)%p(so)%subState0(1:size_so(so),me) & + + sourceState(ph)%p(so)%dotState (1:size_so(so),me) * dt + source_dotState(1:size_so(so),2,so) = 0.0_pReal enddo iteration: do NiterationState = 1, num%nState - do s = 1, phase_Nsources(p) - if(nIterationState > 1) source_dotState(1:size_so(s),2,s) = source_dotState(1:size_so(s),1,s) - source_dotState(1:size_so(s),1,s) = sourceState(p)%p(s)%dotState(:,c) + do so = 1, phase_Nsources(ph) + if(nIterationState > 1) source_dotState(1:size_so(so),2,so) = source_dotState(1:size_so(so),1,so) + source_dotState(1:size_so(so),1,so) = sourceState(ph)%p(so)%dotState(:,me) enddo - broken = constitutive_collectDotState_source(crystallite_S(1:3,1:3,g,i,e), g,i,e,p,c) + broken = constitutive_thermal_collectDotState(ph,me) + broken = broken .or. constitutive_damage_collectDotState(crystallite_S(1:3,1:3,co,ip,el), co,ip,el,ph,me) if(broken) exit iteration - do s = 1, phase_Nsources(p) - zeta = damper(sourceState(p)%p(s)%dotState(:,c), & - source_dotState(1:size_so(s),1,s),& - source_dotState(1:size_so(s),2,s)) - sourceState(p)%p(s)%dotState(:,c) = sourceState(p)%p(s)%dotState(:,c) * zeta & - + source_dotState(1:size_so(s),1,s)* (1.0_pReal - zeta) - r(1:size_so(s)) = sourceState(p)%p(s)%state (1:size_so(s),c) & - - sourceState(p)%p(s)%subState0(1:size_so(s),c) & - - sourceState(p)%p(s)%dotState (1:size_so(s),c) * crystallite_subdt(g,i,e) - sourceState(p)%p(s)%state(1:size_so(s),c) = sourceState(p)%p(s)%state(1:size_so(s),c) & - - r(1:size_so(s)) - crystallite_converged(g,i,e) = & - crystallite_converged(g,i,e) .and. converged(r(1:size_so(s)), & - sourceState(p)%p(s)%state(1:size_so(s),c), & - sourceState(p)%p(s)%atol(1:size_so(s))) + do so = 1, phase_Nsources(ph) + zeta = damper(sourceState(ph)%p(so)%dotState(:,me), & + source_dotState(1:size_so(so),1,so),& + source_dotState(1:size_so(so),2,so)) + sourceState(ph)%p(so)%dotState(:,me) = sourceState(ph)%p(so)%dotState(:,me) * zeta & + + source_dotState(1:size_so(so),1,so)* (1.0_pReal - zeta) + r(1:size_so(so)) = sourceState(ph)%p(so)%state (1:size_so(so),me) & + - sourceState(ph)%p(so)%subState0(1:size_so(so),me) & + - sourceState(ph)%p(so)%dotState (1:size_so(so),me) * dt + sourceState(ph)%p(so)%state(1:size_so(so),me) = sourceState(ph)%p(so)%state(1:size_so(so),me) & + - r(1:size_so(so)) + converged_ = converged_ .and. converged(r(1:size_so(so)), & + sourceState(ph)%p(so)%state(1:size_so(so),me), & + sourceState(ph)%p(so)%atol(1:size_so(so))) enddo - if(crystallite_converged(g,i,e)) then - broken = constitutive_deltaState_source(crystallite_Fe(1:3,1:3,g,i,e),g,i,e,p,c) + if(converged_) then + broken = constitutive_damage_deltaState(crystallite_Fe(1:3,1:3,co,ip,el),co,ip,el,ph,me) exit iteration endif enddo iteration + broken = broken .or. .not. converged_ + contains @@ -1679,7 +1323,7 @@ subroutine integrateSourceState(g,i,e) end function damper -end subroutine integrateSourceState +end function integrateSourceState !-------------------------------------------------------------------------------------------------- @@ -1705,7 +1349,7 @@ end function converged !-------------------------------------------------------------------------------------------------- subroutine crystallite_restartWrite - integer :: i + integer :: ph integer(HID_T) :: fileHandle, groupHandle character(len=pStringLen) :: fileName, datasetName @@ -1714,26 +1358,27 @@ subroutine crystallite_restartWrite write(fileName,'(a,i0,a)') trim(getSolverJobName())//'_',worldrank,'.hdf5' fileHandle = HDF5_openFile(fileName,'a') - call HDF5_write(fileHandle,crystallite_partitionedF,'F') - call HDF5_write(fileHandle,crystallite_Fp, 'F_p') + call HDF5_write(fileHandle,crystallite_F,'F') call HDF5_write(fileHandle,crystallite_Lp, 'L_p') call HDF5_write(fileHandle,crystallite_S, 'S') groupHandle = HDF5_addGroup(fileHandle,'phase') - do i = 1,size(material_name_phase) - write(datasetName,'(i0,a)') i,'_omega' - call HDF5_write(groupHandle,plasticState(i)%state,datasetName) - write(datasetName,'(i0,a)') i,'_F_i' - call HDF5_write(groupHandle,constitutive_mech_Fi(i)%data,datasetName) - write(datasetName,'(i0,a)') i,'_L_i' - call HDF5_write(groupHandle,constitutive_mech_Li(i)%data,datasetName) + do ph = 1,size(material_name_phase) + write(datasetName,'(i0,a)') ph,'_omega' + call HDF5_write(groupHandle,plasticState(ph)%state,datasetName) + write(datasetName,'(i0,a)') ph,'_F_i' + call HDF5_write(groupHandle,constitutive_mech_Fi(ph)%data,datasetName) + write(datasetName,'(i0,a)') ph,'_L_i' + call HDF5_write(groupHandle,constitutive_mech_Li(ph)%data,datasetName) + write(datasetName,'(i0,a)') ph,'_F_p' + call HDF5_write(groupHandle,constitutive_mech_Fp(ph)%data,datasetName) enddo call HDF5_closeGroup(groupHandle) groupHandle = HDF5_addGroup(fileHandle,'homogenization') - do i = 1, size(material_name_homogenization) - write(datasetName,'(i0,a)') i,'_omega' - call HDF5_write(groupHandle,homogState(i)%state,datasetName) + do ph = 1, size(material_name_homogenization) + write(datasetName,'(i0,a)') ph,'_omega' + call HDF5_write(groupHandle,homogState(ph)%state,datasetName) enddo call HDF5_closeGroup(groupHandle) @@ -1748,7 +1393,7 @@ end subroutine crystallite_restartWrite !-------------------------------------------------------------------------------------------------- subroutine crystallite_restartRead - integer :: i + integer :: ph integer(HID_T) :: fileHandle, groupHandle character(len=pStringLen) :: fileName, datasetName @@ -1758,25 +1403,26 @@ subroutine crystallite_restartRead fileHandle = HDF5_openFile(fileName) call HDF5_read(fileHandle,crystallite_F0, 'F') - call HDF5_read(fileHandle,crystallite_Fp0,'F_p') call HDF5_read(fileHandle,crystallite_Lp0,'L_p') call HDF5_read(fileHandle,crystallite_S0, 'S') groupHandle = HDF5_openGroup(fileHandle,'phase') - do i = 1,size(material_name_phase) - write(datasetName,'(i0,a)') i,'_omega' - call HDF5_read(groupHandle,plasticState(i)%state0,datasetName) - write(datasetName,'(i0,a)') i,'_F_i' - call HDF5_read(groupHandle,constitutive_mech_Fi0(i)%data,datasetName) - write(datasetName,'(i0,a)') i,'_L_i' - call HDF5_read(groupHandle,constitutive_mech_Li0(i)%data,datasetName) + do ph = 1,size(material_name_phase) + write(datasetName,'(i0,a)') ph,'_omega' + call HDF5_read(groupHandle,plasticState(ph)%state0,datasetName) + write(datasetName,'(i0,a)') ph,'_F_i' + call HDF5_read(groupHandle,constitutive_mech_Fi0(ph)%data,datasetName) + write(datasetName,'(i0,a)') ph,'_L_i' + call HDF5_read(groupHandle,constitutive_mech_Li0(ph)%data,datasetName) + write(datasetName,'(i0,a)') ph,'_F_p' + call HDF5_read(groupHandle,constitutive_mech_Fp0(ph)%data,datasetName) enddo call HDF5_closeGroup(groupHandle) groupHandle = HDF5_openGroup(fileHandle,'homogenization') - do i = 1,size(material_name_homogenization) - write(datasetName,'(i0,a)') i,'_omega' - call HDF5_read(groupHandle,homogState(i)%state0,datasetName) + do ph = 1,size(material_name_homogenization) + write(datasetName,'(i0,a)') ph,'_omega' + call HDF5_read(groupHandle,homogState(ph)%state0,datasetName) enddo call HDF5_closeGroup(groupHandle) @@ -1785,29 +1431,4 @@ subroutine crystallite_restartRead end subroutine crystallite_restartRead -!-------------------------------------------------------------------------------------------------- -!> @brief Forward data after successful increment. -! ToDo: Any guessing for the current states possible? -!-------------------------------------------------------------------------------------------------- -subroutine crystallite_forward - - integer :: i, j - - crystallite_F0 = crystallite_partitionedF - crystallite_Fp0 = crystallite_Fp - crystallite_Lp0 = crystallite_Lp - crystallite_S0 = crystallite_S - - do i = 1, size(plasticState) - plasticState(i)%state0 = plasticState(i)%state - constitutive_mech_Fi0(i) = constitutive_mech_Fi(i) - constitutive_mech_Li0(i) = constitutive_mech_Li(i) - enddo - do i = 1,size(material_name_homogenization) - homogState (i)%state0 = homogState (i)%state - damageState (i)%state0 = damageState (i)%state - enddo - -end subroutine crystallite_forward - end module constitutive diff --git a/src/constitutive_damage.f90 b/src/constitutive_damage.f90 index a864ca1b8..3ce614666 100644 --- a/src/constitutive_damage.f90 +++ b/src/constitutive_damage.f90 @@ -1,5 +1,5 @@ !---------------------------------------------------------------------------------------------------- -!> @brief internal microstructure state for all damage sources and kinematics constitutive models +!> @brief internal microstructure state for all damage sources and kinematics constitutive models !---------------------------------------------------------------------------------------------------- submodule(constitutive) constitutive_damage @@ -8,7 +8,7 @@ submodule(constitutive) constitutive_damage module function source_damage_anisoBrittle_init(source_length) result(mySources) integer, intent(in) :: source_length logical, dimension(:,:), allocatable :: mySources - end function source_damage_anisoBrittle_init + end function source_damage_anisoBrittle_init module function source_damage_anisoDuctile_init(source_length) result(mySources) integer, intent(in) :: source_length @@ -23,7 +23,7 @@ submodule(constitutive) constitutive_damage module function source_damage_isoDuctile_init(source_length) result(mySources) integer, intent(in) :: source_length logical, dimension(:,:), allocatable :: mySources - end function source_damage_isoDuctile_init + end function source_damage_isoDuctile_init module function kinematics_cleavage_opening_init(kinematics_length) result(myKinematics) integer, intent(in) :: kinematics_length @@ -39,14 +39,14 @@ submodule(constitutive) constitutive_damage module subroutine source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) integer, intent(in) :: & phase, & !< phase ID of element - constituent !< position of element within its phase instance + constituent !< position of element within its phase instance real(pReal), intent(in) :: & - phi !< damage parameter + phi !< damage parameter real(pReal), intent(out) :: & localphiDot, & dLocalphiDot_dPhi end subroutine source_damage_anisoBrittle_getRateAndItsTangent - + module subroutine source_damage_anisoDuctile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) integer, intent(in) :: & phase, & !< phase ID of element @@ -129,7 +129,7 @@ module subroutine damage_init allocate(sourceState(ph)%p(phase_Nsources(ph))) enddo - allocate(phase_source(maxval(phase_Nsources),phases%length), source = SOURCE_undefined_ID) + allocate(phase_source(maxval(phase_Nsources),phases%length), source = SOURCE_undefined_ID) ! initialize source mechanisms if(maxval(phase_Nsources) /= 0) then @@ -141,19 +141,19 @@ module subroutine damage_init !-------------------------------------------------------------------------------------------------- ! initialize kinematic mechanisms - allocate(phase_Nkinematics(phases%length),source = 0) + allocate(phase_Nkinematics(phases%length),source = 0) do ph = 1,phases%length phase => phases%get(ph) kinematics => phase%get('kinematics',defaultVal=emptyList) phase_Nkinematics(ph) = kinematics%length enddo - - allocate(phase_kinematics(maxval(phase_Nkinematics),phases%length), source = KINEMATICS_undefined_ID) + + allocate(phase_kinematics(maxval(phase_Nkinematics),phases%length), source = KINEMATICS_undefined_ID) if(maxval(phase_Nkinematics) /= 0) then where(kinematics_cleavage_opening_init(maxval(phase_Nkinematics))) phase_kinematics = KINEMATICS_cleavage_opening_ID where(kinematics_slipplane_opening_init(maxval(phase_Nkinematics))) phase_kinematics = KINEMATICS_slipplane_opening_ID - endif + endif end subroutine damage_init @@ -167,7 +167,7 @@ module subroutine constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi ip, & !< integration point number el !< element number real(pReal), intent(in) :: & - phi !< damage parameter + phi !< damage parameter real(pReal), intent(inout) :: & phiDot, & dPhiDot_dPhi @@ -183,7 +183,7 @@ module subroutine constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi phiDot = 0.0_pReal dPhiDot_dPhi = 0.0_pReal - + do grain = 1, homogenization_Nconstituents(material_homogenizationAt(el)) phase = material_phaseAt(grain,el) constituent = material_phasememberAt(grain,ip,el) @@ -217,32 +217,35 @@ end subroutine constitutive_damage_getRateAndItsTangents !---------------------------------------------------------------------------------------------- !< @brief writes damage sources results to HDF5 output file !---------------------------------------------------------------------------------------------- -module subroutine damage_results +module subroutine damage_results(group,ph) - integer :: p,i - character(len=pStringLen) :: group + character(len=*), intent(in) :: group + integer, intent(in) :: ph - do p = 1, size(material_name_phase) + integer :: so - sourceLoop: do i = 1, phase_Nsources(p) - group = trim('current/phase')//'/'//trim(material_name_phase(p)) - group = trim(group)//'/sources' - call results_closeGroup(results_addGroup(group)) + sourceLoop: do so = 1, phase_Nsources(ph) - sourceType: select case (phase_source(i,p)) + if (phase_source(so,ph) /= SOURCE_UNDEFINED_ID) & + call results_closeGroup(results_addGroup(group//'sources/')) ! should be 'damage' - case (SOURCE_damage_anisoBrittle_ID) sourceType - call source_damage_anisoBrittle_results(p,group) - case (SOURCE_damage_anisoDuctile_ID) sourceType - call source_damage_anisoDuctile_results(p,group) - case (SOURCE_damage_isoBrittle_ID) sourceType - call source_damage_isoBrittle_results(p,group) - case (SOURCE_damage_isoDuctile_ID) sourceType - call source_damage_isoDuctile_results(p,group) - end select sourceType + sourceType: select case (phase_source(so,ph)) - enddo SourceLoop - enddo + case (SOURCE_damage_anisoBrittle_ID) sourceType + call source_damage_anisoBrittle_results(ph,group//'sources/') + + case (SOURCE_damage_anisoDuctile_ID) sourceType + call source_damage_anisoDuctile_results(ph,group//'sources/') + + case (SOURCE_damage_isoBrittle_ID) sourceType + call source_damage_isoBrittle_results(ph,group//'sources/') + + case (SOURCE_damage_isoDuctile_ID) sourceType + call source_damage_isoDuctile_results(ph,group//'sources/') + + end select sourceType + + enddo SourceLoop end subroutine damage_results diff --git a/src/constitutive_mech.f90 b/src/constitutive_mech.f90 index 07e48b537..c48c59ec9 100644 --- a/src/constitutive_mech.f90 +++ b/src/constitutive_mech.f90 @@ -3,6 +3,13 @@ !---------------------------------------------------------------------------------------------------- submodule(constitutive) constitutive_mech + enum, bind(c); enumerator :: & + ELASTICITY_UNDEFINED_ID, & + ELASTICITY_HOOKE_ID, & + STIFFNESS_DEGRADATION_UNDEFINED_ID, & + STIFFNESS_DEGRADATION_DAMAGE_ID + end enum + integer(kind(ELASTICITY_undefined_ID)), dimension(:), allocatable :: & phase_elasticity !< elasticity of each phase integer(kind(SOURCE_undefined_ID)), dimension(:,:), allocatable :: & @@ -133,7 +140,7 @@ submodule(constitutive) constitutive_mech el !< current element number end subroutine plastic_nonlocal_LpAndItsTangent - module subroutine plastic_isotropic_dotState(Mp,instance,of) + module subroutine plastic_isotropic_dotState(Mp,instance,of) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress integer, intent(in) :: & @@ -177,13 +184,12 @@ submodule(constitutive) constitutive_mech of end subroutine plastic_disloTungsten_dotState - module subroutine plastic_nonlocal_dotState(Mp, F, Fp, Temperature,timestep, & + module subroutine plastic_nonlocal_dotState(Mp, F, Temperature,timestep, & instance,of,ip,el) real(pReal), dimension(3,3), intent(in) :: & Mp !< MandelStress real(pReal), dimension(3,3,homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems), intent(in) :: & - F, & !< deformation gradient - Fp !< plastic deformation gradient + F !< deformation gradient real(pReal), intent(in) :: & Temperature, & !< temperature timestep !< substepped crystallite time increment @@ -209,10 +215,9 @@ submodule(constitutive) constitutive_mech of end subroutine plastic_dislotungsten_dependentState - module subroutine plastic_nonlocal_dependentState(F, Fp, instance, of, ip, el) + module subroutine plastic_nonlocal_dependentState(F, instance, of, ip, el) real(pReal), dimension(3,3), intent(in) :: & - F, & !< deformation gradient - Fp !< plastic deformation gradient + F !< deformation gradient integer, intent(in) :: & instance, & of, & @@ -220,7 +225,7 @@ submodule(constitutive) constitutive_mech el !< current element number end subroutine plastic_nonlocal_dependentState - module subroutine plastic_kinehardening_deltaState(Mp,instance,of) + module subroutine plastic_kinehardening_deltaState(Mp,instance,of) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress integer, intent(in) :: & @@ -268,9 +273,24 @@ submodule(constitutive) constitutive_mech character(len=*), intent(in) :: group end subroutine plastic_nonlocal_results + module function plastic_dislotwin_homogenizedC(co,ip,el) result(homogenizedC) + real(pReal), dimension(6,6) :: & + homogenizedC + integer, intent(in) :: & + co, & !< component-ID of integration point + ip, & !< integration point + el !< element + end function plastic_dislotwin_homogenizedC + end interface + type :: tOutput !< new requested output (per phase) + character(len=pStringLen), allocatable, dimension(:) :: & + label + end type tOutput + type(tOutput), allocatable, dimension(:) :: output_constituent + procedure(integrateStateFPI), pointer :: integrateState contains @@ -282,9 +302,10 @@ contains module subroutine mech_init integer :: & - p, & + ph, & stiffDegradationCtr class(tNode), pointer :: & + num_crystallite, & phases, & phase, & mech, & @@ -299,31 +320,37 @@ module subroutine mech_init allocate(phase_elasticity(phases%length), source = ELASTICITY_undefined_ID) allocate(phase_elasticityInstance(phases%length), source = 0) allocate(phase_NstiffnessDegradations(phases%length),source=0) + allocate(output_constituent(phases%length)) - do p = 1, phases%length - phase => phases%get(p) + do ph = 1, phases%length + phase => phases%get(ph) mech => phase%get('mechanics') +#if defined(__GFORTRAN__) + output_constituent(ph)%label = output_asStrings(mech) +#else + output_constituent(ph)%label = mech%get_asStrings('output',defaultVal=emptyStringArray) +#endif elastic => mech%get('elasticity') if(elastic%get_asString('type') == 'hooke') then - phase_elasticity(p) = ELASTICITY_HOOKE_ID + phase_elasticity(ph) = ELASTICITY_HOOKE_ID else call IO_error(200,ext_msg=elastic%get_asString('type')) endif stiffDegradation => mech%get('stiffness_degradation',defaultVal=emptyList) ! check for stiffness degradation mechanisms - phase_NstiffnessDegradations(p) = stiffDegradation%length + phase_NstiffnessDegradations(ph) = stiffDegradation%length enddo allocate(phase_stiffnessDegradation(maxval(phase_NstiffnessDegradations),phases%length), & source=STIFFNESS_DEGRADATION_undefined_ID) if(maxVal(phase_NstiffnessDegradations)/=0) then - do p = 1, phases%length - phase => phases%get(p) + do ph = 1, phases%length + phase => phases%get(ph) mech => phase%get('mechanics') stiffDegradation => mech%get('stiffness_degradation',defaultVal=emptyList) do stiffDegradationCtr = 1, stiffDegradation%length if(stiffDegradation%get_asString(stiffDegradationCtr) == 'damage') & - phase_stiffnessDegradation(stiffDegradationCtr,p) = STIFFNESS_DEGRADATION_damage_ID + phase_stiffnessDegradation(stiffDegradationCtr,ph) = STIFFNESS_DEGRADATION_damage_ID enddo enddo endif @@ -343,18 +370,42 @@ module subroutine mech_init where(plastic_dislotungsten_init()) phase_plasticity = PLASTICITY_DISLOTUNGSTEN_ID where(plastic_nonlocal_init()) phase_plasticity = PLASTICITY_NONLOCAL_ID - do p = 1, phases%length - phase_elasticityInstance(p) = count(phase_elasticity(1:p) == phase_elasticity(p)) - phase_plasticityInstance(p) = count(phase_plasticity(1:p) == phase_plasticity(p)) + do ph = 1, phases%length + phase_elasticityInstance(ph) = count(phase_elasticity(1:ph) == phase_elasticity(ph)) + phase_plasticityInstance(ph) = count(phase_plasticity(1:ph) == phase_plasticity(ph)) enddo + num_crystallite => config_numerics%get('crystallite',defaultVal=emptyDict) + + select case(num_crystallite%get_asString('integrator',defaultVal='FPI')) + + case('FPI') + integrateState => integrateStateFPI + + case('Euler') + integrateState => integrateStateEuler + + case('AdaptiveEuler') + integrateState => integrateStateAdaptiveEuler + + case('RK4') + integrateState => integrateStateRK4 + + case('RKCK45') + integrateState => integrateStateRKCK45 + + case default + call IO_error(301,ext_msg='integrator') + + end select + end subroutine mech_init !-------------------------------------------------------------------------------------------------- !> @brief checks if a plastic module is active or not !-------------------------------------------------------------------------------------------------- -module function plastic_active(plastic_label) result(active_plastic) +function plastic_active(plastic_label) result(active_plastic) character(len=*), intent(in) :: plastic_label !< type of plasticity model logical, dimension(:), allocatable :: active_plastic @@ -364,15 +415,15 @@ module function plastic_active(plastic_label) result(active_plastic) phase, & mech, & pl - integer :: p + integer :: ph phases => config_material%get('phase') allocate(active_plastic(phases%length), source = .false. ) - do p = 1, phases%length - phase => phases%get(p) + do ph = 1, phases%length + phase => phases%get(ph) mech => phase%get('mechanics') pl => mech%get('plasticity') - if(pl%get_asString('type') == plastic_label) active_plastic(p) = .true. + if(pl%get_asString('type') == plastic_label) active_plastic(ph) = .true. enddo end function plastic_active @@ -383,10 +434,10 @@ end function plastic_active !> the elastic and intermediate deformation gradients using Hooke's law !-------------------------------------------------------------------------------------------------- module subroutine constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, & - Fe, Fi, ipc, ip, el) + Fe, Fi, co, ip, el) integer, intent(in) :: & - ipc, & !< component-ID of integration point + co, & !< component-ID of integration point ip, & !< integration point el !< element real(pReal), intent(in), dimension(3,3) :: & @@ -406,10 +457,10 @@ module subroutine constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, & i, j ho = material_homogenizationAt(el) - C = math_66toSym3333(constitutive_homogenizedC(ipc,ip,el)) + C = math_66toSym3333(constitutive_homogenizedC(co,ip,el)) - DegradationLoop: do d = 1, phase_NstiffnessDegradations(material_phaseAt(ipc,el)) - degradationType: select case(phase_stiffnessDegradation(d,material_phaseAt(ipc,el))) + DegradationLoop: do d = 1, phase_NstiffnessDegradations(material_phaseAt(co,el)) + degradationType: select case(phase_stiffnessDegradation(d,material_phaseAt(co,el))) case (STIFFNESS_DEGRADATION_damage_ID) degradationType C = C * damage(ho)%p(material_homogenizationMemberAt(ip,el))**2 end select degradationType @@ -429,15 +480,14 @@ end subroutine constitutive_hooke_SandItsTangents !-------------------------------------------------------------------------------------------------- !> @brief calls microstructure function of the different plasticity constitutive models !-------------------------------------------------------------------------------------------------- -module subroutine constitutive_plastic_dependentState(F, Fp, ipc, ip, el) +module subroutine constitutive_plastic_dependentState(F, co, ip, el) integer, intent(in) :: & - ipc, & !< component-ID of integration point + co, & !< component-ID of integration point ip, & !< integration point el !< element real(pReal), intent(in), dimension(3,3) :: & - F, & !< elastic deformation gradient - Fp !< plastic deformation gradient + F !< deformation gradient integer :: & ho, & !< homogenization @@ -446,16 +496,16 @@ module subroutine constitutive_plastic_dependentState(F, Fp, ipc, ip, el) ho = material_homogenizationAt(el) tme = material_homogenizationMemberAt(ip,el) - of = material_phasememberAt(ipc,ip,el) - instance = phase_plasticityInstance(material_phaseAt(ipc,el)) + of = material_phasememberAt(co,ip,el) + instance = phase_plasticityInstance(material_phaseAt(co,el)) - plasticityType: select case (phase_plasticity(material_phaseAt(ipc,el))) + plasticityType: select case (phase_plasticity(material_phaseAt(co,el))) case (PLASTICITY_DISLOTWIN_ID) plasticityType call plastic_dislotwin_dependentState(temperature(ho)%p(tme),instance,of) case (PLASTICITY_DISLOTUNGSTEN_ID) plasticityType call plastic_dislotungsten_dependentState(instance,of) case (PLASTICITY_NONLOCAL_ID) plasticityType - call plastic_nonlocal_dependentState (F,Fp,instance,of,ip,el) + call plastic_nonlocal_dependentState (F,instance,of,ip,el) end select plasticityType end subroutine constitutive_plastic_dependentState @@ -467,9 +517,9 @@ end subroutine constitutive_plastic_dependentState ! Mp in, dLp_dMp out !-------------------------------------------------------------------------------------------------- module subroutine constitutive_plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & - S, Fi, ipc, ip, el) + S, Fi, co, ip, el) integer, intent(in) :: & - ipc, & !< component-ID of integration point + co, & !< component-ID of integration point ip, & !< integration point el !< element real(pReal), intent(in), dimension(3,3) :: & @@ -495,10 +545,10 @@ module subroutine constitutive_plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & tme = material_homogenizationMemberAt(ip,el) Mp = matmul(matmul(transpose(Fi),Fi),S) - of = material_phasememberAt(ipc,ip,el) - instance = phase_plasticityInstance(material_phaseAt(ipc,el)) + of = material_phasememberAt(co,ip,el) + instance = phase_plasticityInstance(material_phaseAt(co,el)) - plasticityType: select case (phase_plasticity(material_phaseAt(ipc,el))) + plasticityType: select case (phase_plasticity(material_phaseAt(co,el))) case (PLASTICITY_NONE_ID) plasticityType Lp = 0.0_pReal @@ -536,23 +586,16 @@ end subroutine constitutive_plastic_LpAndItsTangents !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the rate of change of microstructure !-------------------------------------------------------------------------------------------------- -module function constitutive_collectDotState(S, FArray, Fi, FpArray, subdt, ipc, ip, el,phase,of) result(broken) +function mech_collectDotState(subdt, co, ip, el,ph,of) result(broken) integer, intent(in) :: & - ipc, & !< component-ID of integration point + co, & !< component-ID of integration point ip, & !< integration point el, & !< element - phase, & + ph, & of real(pReal), intent(in) :: & subdt !< timestep - real(pReal), intent(in), dimension(3,3,homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems) :: & - FArray, & !< elastic deformation gradient - FpArray !< plastic deformation gradient - real(pReal), intent(in), dimension(3,3) :: & - Fi !< intermediate deformation gradient - real(pReal), intent(in), dimension(3,3) :: & - S !< 2nd Piola Kirchhoff stress (vector notation) real(pReal), dimension(3,3) :: & Mp integer :: & @@ -561,14 +604,14 @@ module function constitutive_collectDotState(S, FArray, Fi, FpArray, subdt, ipc, i, & !< counter in source loop instance logical :: broken - ho = material_homogenizationAt(el) tme = material_homogenizationMemberAt(ip,el) - instance = phase_plasticityInstance(phase) + instance = phase_plasticityInstance(ph) - Mp = matmul(matmul(transpose(Fi),Fi),S) + Mp = matmul(matmul(transpose(constitutive_mech_Fi(ph)%data(1:3,1:3,of)),& + constitutive_mech_Fi(ph)%data(1:3,1:3,of)),crystallite_S(1:3,1:3,co,ip,el)) - plasticityType: select case (phase_plasticity(phase)) + plasticityType: select case (phase_plasticity(ph)) case (PLASTICITY_ISOTROPIC_ID) plasticityType call plastic_isotropic_dotState(Mp,instance,of) @@ -586,26 +629,26 @@ module function constitutive_collectDotState(S, FArray, Fi, FpArray, subdt, ipc, call plastic_disloTungsten_dotState(Mp,temperature(ho)%p(tme),instance,of) case (PLASTICITY_NONLOCAL_ID) plasticityType - call plastic_nonlocal_dotState(Mp,FArray,FpArray,temperature(ho)%p(tme),subdt, & + call plastic_nonlocal_dotState(Mp,crystallite_partitionedF0,temperature(ho)%p(tme),subdt, & instance,of,ip,el) end select plasticityType - broken = any(IEEE_is_NaN(plasticState(phase)%dotState(:,of))) + broken = any(IEEE_is_NaN(plasticState(ph)%dotState(:,of))) -end function constitutive_collectDotState +end function mech_collectDotState !-------------------------------------------------------------------------------------------------- !> @brief for constitutive models having an instantaneous change of state !> will return false if delta state is not needed/supported by the constitutive model !-------------------------------------------------------------------------------------------------- -module function constitutive_deltaState(S, Fi, ipc, ip, el, phase, of) result(broken) +function constitutive_deltaState(S, Fi, co, ip, el, ph, of) result(broken) integer, intent(in) :: & - ipc, & !< component-ID of integration point + co, & !< component-ID of integration point ip, & !< integration point el, & !< element - phase, & + ph, & of real(pReal), intent(in), dimension(3,3) :: & S, & !< 2nd Piola Kirchhoff stress @@ -620,17 +663,17 @@ module function constitutive_deltaState(S, Fi, ipc, ip, el, phase, of) result(br broken Mp = matmul(matmul(transpose(Fi),Fi),S) - instance = phase_plasticityInstance(phase) + instance = phase_plasticityInstance(ph) - plasticityType: select case (phase_plasticity(phase)) + plasticityType: select case (phase_plasticity(ph)) case (PLASTICITY_KINEHARDENING_ID) plasticityType call plastic_kinehardening_deltaState(Mp,instance,of) - broken = any(IEEE_is_NaN(plasticState(phase)%deltaState(:,of))) + broken = any(IEEE_is_NaN(plasticState(ph)%deltaState(:,of))) case (PLASTICITY_NONLOCAL_ID) plasticityType call plastic_nonlocal_deltaState(Mp,instance,of,ip,el) - broken = any(IEEE_is_NaN(plasticState(phase)%deltaState(:,of))) + broken = any(IEEE_is_NaN(plasticState(ph)%deltaState(:,of))) case default broken = .false. @@ -638,73 +681,71 @@ module function constitutive_deltaState(S, Fi, ipc, ip, el, phase, of) result(br end select plasticityType if(.not. broken) then - select case(phase_plasticity(phase)) + select case(phase_plasticity(ph)) case (PLASTICITY_NONLOCAL_ID,PLASTICITY_KINEHARDENING_ID) - myOffset = plasticState(phase)%offsetDeltaState - mySize = plasticState(phase)%sizeDeltaState - plasticState(phase)%state(myOffset + 1:myOffset + mySize,of) = & - plasticState(phase)%state(myOffset + 1:myOffset + mySize,of) + plasticState(phase)%deltaState(1:mySize,of) + myOffset = plasticState(ph)%offsetDeltaState + mySize = plasticState(ph)%sizeDeltaState + plasticState(ph)%state(myOffset + 1:myOffset + mySize,of) = & + plasticState(ph)%state(myOffset + 1:myOffset + mySize,of) + plasticState(ph)%deltaState(1:mySize,of) end select endif end function constitutive_deltaState -!-------------------------------------------------------------------------------------------- -!> @brief writes plasticity constitutive results to HDF5 output file -!-------------------------------------------------------------------------------------------- -module subroutine plastic_results +module subroutine mech_results(group,ph) - integer :: p - character(len=:), allocatable :: group + character(len=*), intent(in) :: group + integer, intent(in) :: ph - plasticityLoop: do p=1,size(material_name_phase) - group = '/current/phase/'//trim(material_name_phase(p)) - call results_closeGroup(results_addGroup(group)) + if (phase_plasticity(ph) /= PLASTICITY_NONE_ID) & + call results_closeGroup(results_addGroup(group//'plastic/')) - group = trim(group)//'/plastic' + select case(phase_plasticity(ph)) - call results_closeGroup(results_addGroup(group)) - select case(phase_plasticity(p)) + case(PLASTICITY_ISOTROPIC_ID) + call plastic_isotropic_results(phase_plasticityInstance(ph),group//'plastic/') - case(PLASTICITY_ISOTROPIC_ID) - call plastic_isotropic_results(phase_plasticityInstance(p),group) + case(PLASTICITY_PHENOPOWERLAW_ID) + call plastic_phenopowerlaw_results(phase_plasticityInstance(ph),group//'plastic/') - case(PLASTICITY_PHENOPOWERLAW_ID) - call plastic_phenopowerlaw_results(phase_plasticityInstance(p),group) + case(PLASTICITY_KINEHARDENING_ID) + call plastic_kinehardening_results(phase_plasticityInstance(ph),group//'plastic/') - case(PLASTICITY_KINEHARDENING_ID) - call plastic_kinehardening_results(phase_plasticityInstance(p),group) + case(PLASTICITY_DISLOTWIN_ID) + call plastic_dislotwin_results(phase_plasticityInstance(ph),group//'plastic/') - case(PLASTICITY_DISLOTWIN_ID) - call plastic_dislotwin_results(phase_plasticityInstance(p),group) + case(PLASTICITY_DISLOTUNGSTEN_ID) + call plastic_dislotungsten_results(phase_plasticityInstance(ph),group//'plastic/') - case(PLASTICITY_DISLOTUNGSTEN_ID) - call plastic_dislotungsten_results(phase_plasticityInstance(p),group) + case(PLASTICITY_NONLOCAL_ID) + call plastic_nonlocal_results(phase_plasticityInstance(ph),group//'plastic/') - case(PLASTICITY_NONLOCAL_ID) - call plastic_nonlocal_results(phase_plasticityInstance(p),group) - end select + end select - enddo plasticityLoop + call crystallite_results(group,ph) -end subroutine plastic_results +end subroutine mech_results + + module subroutine mech_restart_read(fileHandle) + integer(HID_T), intent(in) :: fileHandle + end subroutine mech_restart_read !-------------------------------------------------------------------------------------------------- !> @brief calculation of stress (P) with time integration based on a residuum in Lp and !> intermediate acceleration of the Newton-Raphson correction !-------------------------------------------------------------------------------------------------- -function integrateStress(ipc,ip,el,timeFraction) result(broken) +function integrateStress(F,Delta_t,co,ip,el) result(broken) + real(pReal), dimension(3,3), intent(in) :: F + real(pReal), intent(in) :: Delta_t integer, intent(in):: el, & ! element index ip, & ! integration point index - ipc ! grain index - real(pReal), optional, intent(in) :: timeFraction ! fraction of timestep + co ! grain index - real(pReal), dimension(3,3):: F, & ! deformation gradient at end of timestep - Fp_new, & ! plastic deformation gradient at end of timestep + real(pReal), dimension(3,3):: Fp_new, & ! plastic deformation gradient at end of timestep invFp_new, & ! inverse of Fp_new invFp_current, & ! inverse of Fp_current Lpguess, & ! current guess for plastic velocity gradient @@ -742,7 +783,6 @@ function integrateStress(ipc,ip,el,timeFraction) result(broken) dLi_dS real(pReal) steplengthLp, & steplengthLi, & - dt, & ! time increment atol_Lp, & atol_Li, & devNull @@ -752,33 +792,25 @@ function integrateStress(ipc,ip,el,timeFraction) result(broken) o, & p, & m, & + ph, & + me, & jacoCounterLp, & jacoCounterLi ! counters to check for Jacobian update logical :: error,broken broken = .true. - if (present(timeFraction)) then - dt = crystallite_subdt(ipc,ip,el) * timeFraction - F = crystallite_subF0(1:3,1:3,ipc,ip,el) & - + (crystallite_subF(1:3,1:3,ipc,ip,el) - crystallite_subF0(1:3,1:3,ipc,ip,el)) * timeFraction - else - dt = crystallite_subdt(ipc,ip,el) - F = crystallite_subF(1:3,1:3,ipc,ip,el) - endif + call constitutive_plastic_dependentState(crystallite_F(1:3,1:3,co,ip,el),co,ip,el) - call constitutive_plastic_dependentState(crystallite_partitionedF(1:3,1:3,ipc,ip,el), & - crystallite_Fp(1:3,1:3,ipc,ip,el),ipc,ip,el) + ph = material_phaseAt(co,el) + me = material_phaseMemberAt(co,ip,el) - p = material_phaseAt(ipc,el) - m = material_phaseMemberAt(ipc,ip,el) + Lpguess = crystallite_Lp(1:3,1:3,co,ip,el) ! take as first guess + Liguess = constitutive_mech_Li(ph)%data(1:3,1:3,me) ! take as first guess - Lpguess = crystallite_Lp(1:3,1:3,ipc,ip,el) ! take as first guess - Liguess = constitutive_mech_Li(p)%data(1:3,1:3,m) ! take as first guess - - call math_invert33(invFp_current,devNull,error,crystallite_subFp0(1:3,1:3,ipc,ip,el)) + call math_invert33(invFp_current,devNull,error,crystallite_subFp0(1:3,1:3,co,ip,el)) if (error) return ! error - call math_invert33(invFi_current,devNull,error,crystallite_subFi0(1:3,1:3,ipc,ip,el)) + call math_invert33(invFi_current,devNull,error,crystallite_subFi0(1:3,1:3,co,ip,el)) if (error) return ! error A = matmul(F,invFp_current) ! intermediate tensor needed later to calculate dFe_dLp @@ -793,7 +825,7 @@ function integrateStress(ipc,ip,el,timeFraction) result(broken) NiterationStressLi = NiterationStressLi + 1 if (NiterationStressLi>num%nStress) return ! error - invFi_new = matmul(invFi_current,math_I3 - dt*Liguess) + invFi_new = matmul(invFi_current,math_I3 - Delta_t*Liguess) Fi_new = math_inv33(invFi_new) jacoCounterLp = 0 @@ -806,13 +838,13 @@ function integrateStress(ipc,ip,el,timeFraction) result(broken) NiterationStressLp = NiterationStressLp + 1 if (NiterationStressLp>num%nStress) return ! error - B = math_I3 - dt*Lpguess + B = math_I3 - Delta_t*Lpguess Fe = matmul(matmul(A,B), invFi_new) call constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, & - Fe, Fi_new, ipc, ip, el) + Fe, Fi_new, co, ip, el) call constitutive_plastic_LpAndItsTangents(Lp_constitutive, dLp_dS, dLp_dFi, & - S, Fi_new, ipc, ip, el) + S, Fi_new, co, ip, el) !* update current residuum and check for convergence of loop atol_Lp = max(num%rtol_crystalliteStress * max(norm2(Lpguess),norm2(Lp_constitutive)), & ! absolute tolerance from largest acceptable relative error @@ -838,7 +870,7 @@ function integrateStress(ipc,ip,el,timeFraction) result(broken) jacoCounterLp = jacoCounterLp + 1 do o=1,3; do p=1,3 - dFe_dLp(o,1:3,p,1:3) = - dt * A(o,p)*transpose(invFi_new) ! dFe_dLp(i,j,k,l) = -dt * A(i,k) invFi(l,j) + dFe_dLp(o,1:3,p,1:3) = - Delta_t * A(o,p)*transpose(invFi_new) ! dFe_dLp(i,j,k,l) = -Delta_t * A(i,k) invFi(l,j) enddo; enddo dRLp_dLp = math_eye(9) & - math_3333to99(math_mul3333xx3333(math_mul3333xx3333(dLp_dS,dS_dFe),dFe_dLp)) @@ -853,7 +885,7 @@ function integrateStress(ipc,ip,el,timeFraction) result(broken) enddo LpLoop call constitutive_LiAndItsTangents(Li_constitutive, dLi_dS, dLi_dFi, & - S, Fi_new, ipc, ip, el) + S, Fi_new, co, ip, el) !* update current residuum and check for convergence of loop atol_Li = max(num%rtol_crystalliteStress * max(norm2(Liguess),norm2(Li_constitutive)), & ! absolute tolerance from largest acceptable relative error @@ -879,8 +911,8 @@ function integrateStress(ipc,ip,el,timeFraction) result(broken) temp_33 = matmul(matmul(A,B),invFi_current) do o=1,3; do p=1,3 - dFe_dLi(1:3,o,1:3,p) = -dt*math_I3(o,p)*temp_33 ! dFe_dLp(i,j,k,l) = -dt * A(i,k) invFi(l,j) - dFi_dLi(1:3,o,1:3,p) = -dt*math_I3(o,p)*invFi_current + dFe_dLi(1:3,o,1:3,p) = -Delta_t*math_I3(o,p)*temp_33 ! dFe_dLp(i,j,k,l) = -Delta_t * A(i,k) invFi(l,j) + dFi_dLi(1:3,o,1:3,p) = -Delta_t*math_I3(o,p)*invFi_current enddo; enddo do o=1,3; do p=1,3 dFi_dLi(1:3,1:3,o,p) = matmul(matmul(Fi_new,dFi_dLi(1:3,1:3,o,p)),Fi_new) @@ -903,16 +935,13 @@ function integrateStress(ipc,ip,el,timeFraction) result(broken) call math_invert33(Fp_new,devNull,error,invFp_new) if (error) return ! error - p = material_phaseAt(ipc,el) - m = material_phaseMemberAt(ipc,ip,el) - - crystallite_P (1:3,1:3,ipc,ip,el) = matmul(matmul(F,invFp_new),matmul(S,transpose(invFp_new))) - crystallite_S (1:3,1:3,ipc,ip,el) = S - crystallite_Lp (1:3,1:3,ipc,ip,el) = Lpguess - constitutive_mech_Li(p)%data(1:3,1:3,m) = Liguess - crystallite_Fp (1:3,1:3,ipc,ip,el) = Fp_new / math_det33(Fp_new)**(1.0_pReal/3.0_pReal) ! regularize - constitutive_mech_Fi(p)%data(1:3,1:3,m) = Fi_new - crystallite_Fe (1:3,1:3,ipc,ip,el) = matmul(matmul(F,invFp_new),invFi_new) + crystallite_P (1:3,1:3,co,ip,el) = matmul(matmul(F,invFp_new),matmul(S,transpose(invFp_new))) + crystallite_S (1:3,1:3,co,ip,el) = S + crystallite_Lp (1:3,1:3,co,ip,el) = Lpguess + constitutive_mech_Li(ph)%data(1:3,1:3,me) = Liguess + constitutive_mech_Fp(ph)%data(1:3,1:3,me) = Fp_new / math_det33(Fp_new)**(1.0_pReal/3.0_pReal) ! regularize + constitutive_mech_Fi(ph)%data(1:3,1:3,me) = Fi_new + crystallite_Fe (1:3,1:3,co,ip,el) = matmul(matmul(F,invFp_new),invFi_new) broken = .false. end function integrateStress @@ -922,77 +951,64 @@ end function integrateStress !> @brief integrate stress, state with adaptive 1st order explicit Euler method !> using Fixed Point Iteration to adapt the stepsize !-------------------------------------------------------------------------------------------------- -module subroutine integrateStateFPI(g,i,e) +function integrateStateFPI(F_0,F,Delta_t,co,ip,el) result(broken) + real(pReal), intent(in),dimension(3,3) :: F_0,F + real(pReal), intent(in) :: Delta_t integer, intent(in) :: & - e, & !< element index in element loop - i, & !< integration point index in ip loop - g !< grain index in grain loop - integer :: & - NiterationState, & !< number of iterations in state loop - p, & - c, & - s, & - size_pl - integer, dimension(maxval(phase_Nsources)) :: & - size_so - real(pReal) :: & - zeta - real(pReal), dimension(max(constitutive_plasticity_maxSizeDotState,constitutive_source_maxSizeDotState)) :: & - r ! state residuum - real(pReal), dimension(constitutive_plasticity_maxSizeDotState,2) :: & - plastic_dotState - real(pReal), dimension(constitutive_source_maxSizeDotState,2,maxval(phase_Nsources)) :: source_dotState + el, & !< element index in element loop + ip, & !< integration point index in ip loop + co !< grain index in grain loop logical :: & broken - p = material_phaseAt(g,e) - c = material_phaseMemberAt(g,i,e) + integer :: & + NiterationState, & !< number of iterations in state loop + ph, & + me, & + sizeDotState + real(pReal) :: & + zeta + real(pReal), dimension(constitutive_plasticity_maxSizeDotState) :: & + r ! state residuum + real(pReal), dimension(constitutive_plasticity_maxSizeDotState,2) :: & + dotState - broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & - crystallite_partitionedF0, & - constitutive_mech_Fi(p)%data(1:3,1:3,c), & - crystallite_partitionedFp0, & - crystallite_subdt(g,i,e), g,i,e,p,c) + + ph = material_phaseAt(co,el) + me = material_phaseMemberAt(co,ip,el) + + broken = mech_collectDotState(Delta_t, co,ip,el,ph,me) if(broken) return - size_pl = plasticState(p)%sizeDotState - plasticState(p)%state(1:size_pl,c) = plasticState(p)%subState0(1:size_pl,c) & - + plasticState(p)%dotState (1:size_pl,c) & - * crystallite_subdt(g,i,e) - plastic_dotState(1:size_pl,2) = 0.0_pReal + sizeDotState = plasticState(ph)%sizeDotState + plasticState(ph)%state(1:sizeDotState,me) = plasticState(ph)%subState0(1:sizeDotState,me) & + + plasticState(ph)%dotState (1:sizeDotState,me) * Delta_t + dotState(1:sizeDotState,2) = 0.0_pReal iteration: do NiterationState = 1, num%nState - if(nIterationState > 1) plastic_dotState(1:size_pl,2) = plastic_dotState(1:size_pl,1) - plastic_dotState(1:size_pl,1) = plasticState(p)%dotState(:,c) + if(nIterationState > 1) dotState(1:sizeDotState,2) = dotState(1:sizeDotState,1) + dotState(1:sizeDotState,1) = plasticState(ph)%dotState(:,me) - broken = integrateStress(g,i,e) + broken = integrateStress(F,Delta_t,co,ip,el) if(broken) exit iteration - broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & - crystallite_partitionedF0, & - constitutive_mech_Fi(p)%data(1:3,1:3,c), & - crystallite_partitionedFp0, & - crystallite_subdt(g,i,e), g,i,e,p,c) + broken = mech_collectDotState(Delta_t, co,ip,el,ph,me) if(broken) exit iteration - zeta = damper(plasticState(p)%dotState(:,c),plastic_dotState(1:size_pl,1),& - plastic_dotState(1:size_pl,2)) - plasticState(p)%dotState(:,c) = plasticState(p)%dotState(:,c) * zeta & - + plastic_dotState(1:size_pl,1) * (1.0_pReal - zeta) - r(1:size_pl) = plasticState(p)%state (1:size_pl,c) & - - plasticState(p)%subState0(1:size_pl,c) & - - plasticState(p)%dotState (1:size_pl,c) * crystallite_subdt(g,i,e) - plasticState(p)%state(1:size_pl,c) = plasticState(p)%state(1:size_pl,c) & - - r(1:size_pl) - crystallite_converged(g,i,e) = converged(r(1:size_pl), & - plasticState(p)%state(1:size_pl,c), & - plasticState(p)%atol(1:size_pl)) - - if(crystallite_converged(g,i,e)) then - broken = constitutive_deltaState(crystallite_S(1:3,1:3,g,i,e), & - constitutive_mech_Fi(p)%data(1:3,1:3,c),g,i,e,p,c) + zeta = damper(plasticState(ph)%dotState(:,me),dotState(1:sizeDotState,1),& + dotState(1:sizeDotState,2)) + plasticState(ph)%dotState(:,me) = plasticState(ph)%dotState(:,me) * zeta & + + dotState(1:sizeDotState,1) * (1.0_pReal - zeta) + r(1:sizeDotState) = plasticState(ph)%state (1:sizeDotState,me) & + - plasticState(ph)%subState0(1:sizeDotState,me) & + - plasticState(ph)%dotState (1:sizeDotState,me) * Delta_t + plasticState(ph)%state(1:sizeDotState,me) = plasticState(ph)%state(1:sizeDotState,me) & + - r(1:sizeDotState) + if (converged(r(1:sizeDotState),plasticState(ph)%state(1:sizeDotState,me),plasticState(ph)%atol(1:sizeDotState))) then + broken = constitutive_deltaState(crystallite_S(1:3,1:3,co,ip,el), & + constitutive_mech_Fi(ph)%data(1:3,1:3,me),co,ip,el,ph,me) exit iteration endif @@ -1021,115 +1037,107 @@ module subroutine integrateStateFPI(g,i,e) end function damper -end subroutine integrateStateFPI +end function integrateStateFPI !-------------------------------------------------------------------------------------------------- !> @brief integrate state with 1st order explicit Euler method !-------------------------------------------------------------------------------------------------- -module subroutine integrateStateEuler(g,i,e) +function integrateStateEuler(F_0,F,Delta_t,co,ip,el) result(broken) + real(pReal), intent(in),dimension(3,3) :: F_0,F + real(pReal), intent(in) :: Delta_t integer, intent(in) :: & - e, & !< element index in element loop - i, & !< integration point index in ip loop - g !< grain index in grain loop - integer :: & - p, & - c, & - sizeDotState + el, & !< element index in element loop + ip, & !< integration point index in ip loop + co !< grain index in grain loop logical :: & broken - p = material_phaseAt(g,e) - c = material_phaseMemberAt(g,i,e) + integer :: & + ph, & + me, & + sizeDotState - broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & - crystallite_partitionedF0, & - constitutive_mech_Fi(p)%data(1:3,1:3,c), & - crystallite_partitionedFp0, & - crystallite_subdt(g,i,e), g,i,e,p,c) + + ph = material_phaseAt(co,el) + me = material_phaseMemberAt(co,ip,el) + + broken = mech_collectDotState(Delta_t, co,ip,el,ph,me) if(broken) return - sizeDotState = plasticState(p)%sizeDotState - plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) & - + plasticState(p)%dotState (1:sizeDotState,c) & - * crystallite_subdt(g,i,e) + sizeDotState = plasticState(ph)%sizeDotState + plasticState(ph)%state(1:sizeDotState,me) = plasticState(ph)%subState0(1:sizeDotState,me) & + + plasticState(ph)%dotState (1:sizeDotState,me) * Delta_t - broken = constitutive_deltaState(crystallite_S(1:3,1:3,g,i,e), & - constitutive_mech_Fi(p)%data(1:3,1:3,c),g,i,e,p,c) + broken = constitutive_deltaState(crystallite_S(1:3,1:3,co,ip,el), & + constitutive_mech_Fi(ph)%data(1:3,1:3,me),co,ip,el,ph,me) if(broken) return - broken = integrateStress(g,i,e) - crystallite_converged(g,i,e) = .not. broken + broken = integrateStress(F,Delta_t,co,ip,el) -end subroutine integrateStateEuler +end function integrateStateEuler !-------------------------------------------------------------------------------------------------- !> @brief integrate stress, state with 1st order Euler method with adaptive step size !-------------------------------------------------------------------------------------------------- -module subroutine integrateStateAdaptiveEuler(g,i,e) +function integrateStateAdaptiveEuler(F_0,F,Delta_t,co,ip,el) result(broken) + real(pReal), intent(in),dimension(3,3) :: F_0,F + real(pReal), intent(in) :: Delta_t integer, intent(in) :: & - e, & !< element index in element loop - i, & !< integration point index in ip loop - g !< grain index in grain loop - integer :: & - p, & - c, & - sizeDotState + el, & !< element index in element loop + ip, & !< integration point index in ip loop + co !< grain index in grain loop logical :: & broken + integer :: & + ph, & + me, & + sizeDotState real(pReal), dimension(constitutive_plasticity_maxSizeDotState) :: residuum_plastic - p = material_phaseAt(g,e) - c = material_phaseMemberAt(g,i,e) + ph = material_phaseAt(co,el) + me = material_phaseMemberAt(co,ip,el) - broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & - crystallite_partitionedF0, & - constitutive_mech_Fi(p)%data(1:3,1:3,c), & - crystallite_partitionedFp0, & - crystallite_subdt(g,i,e), g,i,e,p,c) + broken = mech_collectDotState(Delta_t, co,ip,el,ph,me) if(broken) return - sizeDotState = plasticState(p)%sizeDotState + sizeDotState = plasticState(ph)%sizeDotState - residuum_plastic(1:sizeDotState) = - plasticState(p)%dotstate(1:sizeDotState,c) * 0.5_pReal * crystallite_subdt(g,i,e) - plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) & - + plasticState(p)%dotstate(1:sizeDotState,c) * crystallite_subdt(g,i,e) + residuum_plastic(1:sizeDotState) = - plasticState(ph)%dotstate(1:sizeDotState,me) * 0.5_pReal * Delta_t + plasticState(ph)%state(1:sizeDotState,me) = plasticState(ph)%subState0(1:sizeDotState,me) & + + plasticState(ph)%dotstate(1:sizeDotState,me) * Delta_t - broken = constitutive_deltaState(crystallite_S(1:3,1:3,g,i,e), & - constitutive_mech_Fi(p)%data(1:3,1:3,c),g,i,e,p,c) + broken = constitutive_deltaState(crystallite_S(1:3,1:3,co,ip,el), & + constitutive_mech_Fi(ph)%data(1:3,1:3,me),co,ip,el,ph,me) if(broken) return - broken = integrateStress(g,i,e) + broken = integrateStress(F,Delta_t,co,ip,el) if(broken) return - broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & - crystallite_partitionedF0, & - constitutive_mech_Fi(p)%data(1:3,1:3,c), & - crystallite_partitionedFp0, & - crystallite_subdt(g,i,e), g,i,e,p,c) + broken = mech_collectDotState(Delta_t, co,ip,el,ph,me) if(broken) return + broken = .not. converged(residuum_plastic(1:sizeDotState) + 0.5_pReal * plasticState(ph)%dotState(:,me) * Delta_t, & + plasticState(ph)%state(1:sizeDotState,me), & + plasticState(ph)%atol(1:sizeDotState)) - sizeDotState = plasticState(p)%sizeDotState - crystallite_converged(g,i,e) = converged(residuum_plastic(1:sizeDotState) & - + 0.5_pReal * plasticState(p)%dotState(:,c) * crystallite_subdt(g,i,e), & - plasticState(p)%state(1:sizeDotState,c), & - plasticState(p)%atol(1:sizeDotState)) - -end subroutine integrateStateAdaptiveEuler +end function integrateStateAdaptiveEuler !--------------------------------------------------------------------------------------------------- !> @brief Integrate state (including stress integration) with the classic Runge Kutta method !--------------------------------------------------------------------------------------------------- -module subroutine integrateStateRK4(g,i,e) +function integrateStateRK4(F_0,F,Delta_t,co,ip,el) result(broken) - integer, intent(in) :: g,i,e + real(pReal), intent(in),dimension(3,3) :: F_0,F + real(pReal), intent(in) :: Delta_t + integer, intent(in) :: co,ip,el + logical :: broken real(pReal), dimension(3,3), parameter :: & A = reshape([& @@ -1142,17 +1150,21 @@ module subroutine integrateStateRK4(g,i,e) real(pReal), dimension(4), parameter :: & B = [1.0_pReal/6.0_pReal, 1.0_pReal/3.0_pReal, 1.0_pReal/3.0_pReal, 1.0_pReal/6.0_pReal] - call integrateStateRK(g,i,e,A,B,C) -end subroutine integrateStateRK4 + broken = integrateStateRK(F_0,F,Delta_t,co,ip,el,A,B,C) + +end function integrateStateRK4 !--------------------------------------------------------------------------------------------------- !> @brief Integrate state (including stress integration) with the Cash-Carp method !--------------------------------------------------------------------------------------------------- -module subroutine integrateStateRKCK45(g,i,e) +function integrateStateRKCK45(F_0,F,Delta_t,co,ip,el) result(broken) - integer, intent(in) :: g,i,e + real(pReal), intent(in),dimension(3,3) :: F_0,F + real(pReal), intent(in) :: Delta_t + integer, intent(in) :: co,ip,el + logical :: broken real(pReal), dimension(5,5), parameter :: & A = reshape([& @@ -1172,100 +1184,428 @@ module subroutine integrateStateRKCK45(g,i,e) [2825.0_pReal/27648.0_pReal, .0_pReal, 18575.0_pReal/48384.0_pReal,& 13525.0_pReal/55296.0_pReal, 277.0_pReal/14336.0_pReal, 1._pReal/4._pReal] - call integrateStateRK(g,i,e,A,B,C,DB) -end subroutine integrateStateRKCK45 + broken = integrateStateRK(F_0,F,Delta_t,co,ip,el,A,B,C,DB) + +end function integrateStateRKCK45 !-------------------------------------------------------------------------------------------------- !> @brief Integrate state (including stress integration) with an explicit Runge-Kutta method or an !! embedded explicit Runge-Kutta method !-------------------------------------------------------------------------------------------------- -subroutine integrateStateRK(g,i,e,A,B,CC,DB) - +function integrateStateRK(F_0,F,Delta_t,co,ip,el,A,B,C,DB) result(broken) + real(pReal), intent(in),dimension(3,3) :: F_0,F + real(pReal), intent(in) :: Delta_t real(pReal), dimension(:,:), intent(in) :: A - real(pReal), dimension(:), intent(in) :: B, CC + real(pReal), dimension(:), intent(in) :: B, C real(pReal), dimension(:), intent(in), optional :: DB - integer, intent(in) :: & - e, & !< element index in element loop - i, & !< integration point index in ip loop - g !< grain index in grain loop + el, & !< element index in element loop + ip, & !< integration point index in ip loop + co !< grain index in grain loop + logical :: broken + integer :: & stage, & ! stage index in integration stage loop n, & - p, & - c, & + ph, & + me, & sizeDotState - logical :: & - broken - real(pReal), dimension(constitutive_plasticity_maxSizeDotState,size(B)) :: plastic_RKdotState + real(pReal), dimension(constitutive_plasticity_maxSizeDotState,size(B)) :: plastic_RKdotState - p = material_phaseAt(g,e) - c = material_phaseMemberAt(g,i,e) - broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & - crystallite_partitionedF0, & - constitutive_mech_Fi(p)%data(1:3,1:3,c), & - crystallite_partitionedFp0, & - crystallite_subdt(g,i,e), g,i,e,p,c) + ph = material_phaseAt(co,el) + me = material_phaseMemberAt(co,ip,el) + + broken = mech_collectDotState(Delta_t,co,ip,el,ph,me) if(broken) return - do stage = 1,size(A,1) - sizeDotState = plasticState(p)%sizeDotState - plastic_RKdotState(1:sizeDotState,stage) = plasticState(p)%dotState(:,c) - plasticState(p)%dotState(:,c) = A(1,stage) * plastic_RKdotState(1:sizeDotState,1) + sizeDotState = plasticState(ph)%sizeDotState + + do stage = 1, size(A,1) + + plastic_RKdotState(1:sizeDotState,stage) = plasticState(ph)%dotState(:,me) + plasticState(ph)%dotState(:,me) = A(1,stage) * plastic_RKdotState(1:sizeDotState,1) do n = 2, stage - sizeDotState = plasticState(p)%sizeDotState - plasticState(p)%dotState(:,c) = plasticState(p)%dotState(:,c) & - + A(n,stage) * plastic_RKdotState(1:sizeDotState,n) + plasticState(ph)%dotState(:,me) = plasticState(ph)%dotState(:,me) & + + A(n,stage) * plastic_RKdotState(1:sizeDotState,n) enddo - sizeDotState = plasticState(p)%sizeDotState - plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) & - + plasticState(p)%dotState (1:sizeDotState,c) & - * crystallite_subdt(g,i,e) + plasticState(ph)%state(1:sizeDotState,me) = plasticState(ph)%subState0(1:sizeDotState,me) & + + plasticState(ph)%dotState (1:sizeDotState,me) * Delta_t - broken = integrateStress(g,i,e,CC(stage)) + broken = integrateStress(F_0 + (F - F_0) * Delta_t * C(stage),Delta_t * C(stage),co,ip,el) if(broken) exit - broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & - crystallite_partitionedF0, & - constitutive_mech_Fi(p)%data(1:3,1:3,c), & - crystallite_partitionedFp0, & - crystallite_subdt(g,i,e)*CC(stage), g,i,e,p,c) + broken = mech_collectDotState(Delta_t*C(stage),co,ip,el,ph,me) if(broken) exit enddo if(broken) return - sizeDotState = plasticState(p)%sizeDotState - plastic_RKdotState(1:sizeDotState,size(B)) = plasticState (p)%dotState(:,c) - plasticState(p)%dotState(:,c) = matmul(plastic_RKdotState(1:sizeDotState,1:size(B)),B) - plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) & - + plasticState(p)%dotState (1:sizeDotState,c) & - * crystallite_subdt(g,i,e) + plastic_RKdotState(1:sizeDotState,size(B)) = plasticState (ph)%dotState(:,me) + plasticState(ph)%dotState(:,me) = matmul(plastic_RKdotState(1:sizeDotState,1:size(B)),B) + plasticState(ph)%state(1:sizeDotState,me) = plasticState(ph)%subState0(1:sizeDotState,me) & + + plasticState(ph)%dotState (1:sizeDotState,me) * Delta_t + if(present(DB)) & - broken = .not. converged( matmul(plastic_RKdotState(1:sizeDotState,1:size(DB)),DB) & - * crystallite_subdt(g,i,e), & - plasticState(p)%state(1:sizeDotState,c), & - plasticState(p)%atol(1:sizeDotState)) + broken = .not. converged(matmul(plastic_RKdotState(1:sizeDotState,1:size(DB)),DB) * Delta_t, & + plasticState(ph)%state(1:sizeDotState,me), & + plasticState(ph)%atol(1:sizeDotState)) if(broken) return - broken = constitutive_deltaState(crystallite_S(1:3,1:3,g,i,e), & - constitutive_mech_Fi(p)%data(1:3,1:3,c),g,i,e,p,c) + broken = constitutive_deltaState(crystallite_S(1:3,1:3,co,ip,el), & + constitutive_mech_Fi(ph)%data(1:3,1:3,me),co,ip,el,ph,me) if(broken) return - broken = integrateStress(g,i,e) - crystallite_converged(g,i,e) = .not. broken + broken = integrateStress(F,Delta_t,co,ip,el) + +end function integrateStateRK -end subroutine integrateStateRK +!-------------------------------------------------------------------------------------------------- +!> @brief writes crystallite results to HDF5 output file +!-------------------------------------------------------------------------------------------------- +subroutine crystallite_results(group,ph) + character(len=*), intent(in) :: group + integer, intent(in) :: ph + + integer :: ou + real(pReal), allocatable, dimension(:,:,:) :: selected_tensors + real(pReal), allocatable, dimension(:,:) :: selected_rotations + character(len=:), allocatable :: structureLabel + + + call results_closeGroup(results_addGroup(group//'/mechanics/')) + + do ou = 1, size(output_constituent(ph)%label) + + select case (output_constituent(ph)%label(ou)) + case('F') + selected_tensors = select_tensors(crystallite_F,ph) + call results_writeDataset(group//'/mechanics/',selected_tensors,output_constituent(ph)%label(ou),& + 'deformation gradient','1') + case('F_e') + selected_tensors = select_tensors(crystallite_Fe,ph) + call results_writeDataset(group//'/mechanics/',selected_tensors,output_constituent(ph)%label(ou),& + 'elastic deformation gradient','1') + case('F_p') + call results_writeDataset(group//'/mechanics/',constitutive_mech_Fp(ph)%data,output_constituent(ph)%label(ou),& + 'plastic deformation gradient','1') + case('F_i') + call results_writeDataset(group//'/mechanics/',constitutive_mech_Fi(ph)%data,output_constituent(ph)%label(ou),& + 'inelastic deformation gradient','1') + case('L_p') + selected_tensors = select_tensors(crystallite_Lp,ph) + call results_writeDataset(group//'/mechanics/',selected_tensors,output_constituent(ph)%label(ou),& + 'plastic velocity gradient','1/s') + case('L_i') + call results_writeDataset(group//'/mechanics/',constitutive_mech_Li(ph)%data,output_constituent(ph)%label(ou),& + 'inelastic velocity gradient','1/s') + case('P') + selected_tensors = select_tensors(crystallite_P,ph) + call results_writeDataset(group//'/mechanics/',selected_tensors,output_constituent(ph)%label(ou),& + 'First Piola-Kirchhoff stress','Pa') + case('S') + selected_tensors = select_tensors(crystallite_S,ph) + call results_writeDataset(group//'/mechanics/',selected_tensors,output_constituent(ph)%label(ou),& + 'Second Piola-Kirchhoff stress','Pa') + case('O') + select case(lattice_structure(ph)) + case(lattice_ISO_ID) + structureLabel = 'aP' + case(lattice_FCC_ID) + structureLabel = 'cF' + case(lattice_BCC_ID) + structureLabel = 'cI' + case(lattice_BCT_ID) + structureLabel = 'tI' + case(lattice_HEX_ID) + structureLabel = 'hP' + case(lattice_ORT_ID) + structureLabel = 'oP' + end select + selected_rotations = select_rotations(crystallite_orientation,ph) + call results_writeDataset(group//'/mechanics/',selected_rotations,output_constituent(ph)%label(ou),& + 'crystal orientation as quaternion','q_0 (q_1 q_2 q_3)') + call results_addAttribute('Lattice',structureLabel,group//'/mechanics/'//output_constituent(ph)%label(ou)) + end select + enddo + + + contains + + !------------------------------------------------------------------------------------------------ + !> @brief select tensors for output + !------------------------------------------------------------------------------------------------ + function select_tensors(dataset,ph) + + integer, intent(in) :: ph + real(pReal), dimension(:,:,:,:,:), intent(in) :: dataset + real(pReal), allocatable, dimension(:,:,:) :: select_tensors + integer :: el,ip,co,j + + allocate(select_tensors(3,3,count(material_phaseAt==ph)*discretization_nIPs)) + + j=0 + do el = 1, size(material_phaseAt,2) + do ip = 1, discretization_nIPs + do co = 1, size(material_phaseAt,1) !ToDo: this needs to be changed for varying Ngrains + if (material_phaseAt(co,el) == ph) then + j = j + 1 + select_tensors(1:3,1:3,j) = dataset(1:3,1:3,co,ip,el) + endif + enddo + enddo + enddo + + end function select_tensors + + +!-------------------------------------------------------------------------------------------------- +!> @brief select rotations for output +!-------------------------------------------------------------------------------------------------- + function select_rotations(dataset,ph) + + integer, intent(in) :: ph + type(rotation), dimension(:,:,:), intent(in) :: dataset + real(pReal), allocatable, dimension(:,:) :: select_rotations + integer :: el,ip,co,j + + allocate(select_rotations(4,count(material_phaseAt==ph)*homogenization_maxNconstituents*discretization_nIPs)) + + j=0 + do el = 1, size(material_phaseAt,2) + do ip = 1, discretization_nIPs + do co = 1, size(material_phaseAt,1) !ToDo: this needs to be changed for varying Ngrains + if (material_phaseAt(co,el) == ph) then + j = j + 1 + select_rotations(1:4,j) = dataset(co,ip,el)%asQuaternion() + endif + enddo + enddo + enddo + + end function select_rotations + +end subroutine crystallite_results + + +!-------------------------------------------------------------------------------------------------- +!> @brief Backup data for homog cutback. +!-------------------------------------------------------------------------------------------------- +module subroutine mech_initializeRestorationPoints(ph,me) + + integer, intent(in) :: ph, me + + + constitutive_mech_partitionedFi0(ph)%data(1:3,1:3,me) = constitutive_mech_Fi0(ph)%data(1:3,1:3,me) + constitutive_mech_partitionedFp0(ph)%data(1:3,1:3,me) = constitutive_mech_Fp0(ph)%data(1:3,1:3,me) + constitutive_mech_partitionedLi0(ph)%data(1:3,1:3,me) = constitutive_mech_Li0(ph)%data(1:3,1:3,me) + plasticState(ph)%partitionedState0(:,me) = plasticState(ph)%state0(:,me) + +end subroutine mech_initializeRestorationPoints + + +!-------------------------------------------------------------------------------------------------- +!> @brief Wind homog inc forward. +!-------------------------------------------------------------------------------------------------- +module subroutine constitutive_mech_windForward(ph,me) + + integer, intent(in) :: ph, me + + + constitutive_mech_partitionedFp0(ph)%data(1:3,1:3,me) = constitutive_mech_Fp(ph)%data(1:3,1:3,me) + constitutive_mech_partitionedFi0(ph)%data(1:3,1:3,me) = constitutive_mech_Fi(ph)%data(1:3,1:3,me) + constitutive_mech_partitionedLi0(ph)%data(1:3,1:3,me) = constitutive_mech_Li(ph)%data(1:3,1:3,me) + + plasticState(ph)%partitionedState0(:,me) = plasticState(ph)%state(:,me) + +end subroutine constitutive_mech_windForward + + +!-------------------------------------------------------------------------------------------------- +!> @brief Forward data after successful increment. +! ToDo: Any guessing for the current states possible? +!-------------------------------------------------------------------------------------------------- +module subroutine constitutive_mech_forward() + + integer :: ph + + + do ph = 1, size(plasticState) + plasticState(ph)%state0 = plasticState(ph)%state + constitutive_mech_Fi0(ph) = constitutive_mech_Fi(ph) + constitutive_mech_Fp0(ph) = constitutive_mech_Fp(ph) + constitutive_mech_Li0(ph) = constitutive_mech_Li(ph) + enddo + +end subroutine constitutive_mech_forward + + + +!-------------------------------------------------------------------------------------------------- +!> @brief returns the homogenize elasticity matrix +!> ToDo: homogenizedC66 would be more consistent +!-------------------------------------------------------------------------------------------------- +module function constitutive_homogenizedC(co,ip,el) result(C) + + real(pReal), dimension(6,6) :: C + integer, intent(in) :: & + co, & !< component-ID of integration point + ip, & !< integration point + el !< element + + plasticityType: select case (phase_plasticity(material_phaseAt(co,el))) + case (PLASTICITY_DISLOTWIN_ID) plasticityType + C = plastic_dislotwin_homogenizedC(co,ip,el) + case default plasticityType + C = lattice_C66(1:6,1:6,material_phaseAt(co,el)) + end select plasticityType + +end function constitutive_homogenizedC + + +!-------------------------------------------------------------------------------------------------- +!> @brief calculate stress (P) +!-------------------------------------------------------------------------------------------------- +module function crystallite_stress(dt,co,ip,el) result(converged_) + + real(pReal), intent(in) :: dt + integer, intent(in) :: & + co, & + ip, & + el + logical :: converged_ + + real(pReal) :: & + formerSubStep + integer :: & + NiterationCrystallite, & ! number of iterations in crystallite loop + so, ph, me + logical :: todo + real(pReal) :: subFrac,subStep + real(pReal), dimension(3,3) :: & + subLp0, & !< plastic velocity grad at start of crystallite inc + subLi0, & !< intermediate velocity grad at start of crystallite inc + subF0, & + subF + + + ph = material_phaseAt(co,el) + me = material_phaseMemberAt(co,ip,el) + subLi0 = constitutive_mech_partitionedLi0(ph)%data(1:3,1:3,me) + subLp0 = crystallite_partitionedLp0(1:3,1:3,co,ip,el) + + plasticState(ph)%subState0(:,me) = plasticState(ph)%partitionedState0(:,me) + do so = 1, phase_Nsources(ph) + sourceState(ph)%p(so)%subState0(:,me) = sourceState(ph)%p(so)%partitionedState0(:,me) + enddo + crystallite_subFp0(1:3,1:3,co,ip,el) = constitutive_mech_partitionedFp0(ph)%data(1:3,1:3,me) + crystallite_subFi0(1:3,1:3,co,ip,el) = constitutive_mech_partitionedFi0(ph)%data(1:3,1:3,me) + subF0 = crystallite_partitionedF0(1:3,1:3,co,ip,el) + subFrac = 0.0_pReal + subStep = 1.0_pReal/num%subStepSizeCryst + todo = .true. + converged_ = .false. ! pretend failed step of 1/subStepSizeCryst + + todo = .true. + NiterationCrystallite = 0 + cutbackLooping: do while (todo) + NiterationCrystallite = NiterationCrystallite + 1 + +!-------------------------------------------------------------------------------------------------- +! wind forward + if (converged_) then + formerSubStep = subStep + subFrac = subFrac + subStep + subStep = min(1.0_pReal - subFrac, num%stepIncreaseCryst * subStep) + + todo = subStep > 0.0_pReal ! still time left to integrate on? + + if (todo) then + subF0 = subF + subLp0 = crystallite_Lp (1:3,1:3,co,ip,el) + subLi0 = constitutive_mech_Li(ph)%data(1:3,1:3,me) + crystallite_subFp0(1:3,1:3,co,ip,el) = constitutive_mech_Fp(ph)%data(1:3,1:3,me) + crystallite_subFi0(1:3,1:3,co,ip,el) = constitutive_mech_Fi(ph)%data(1:3,1:3,me) + plasticState(ph)%subState0(:,me) = plasticState(ph)%state(:,me) + do so = 1, phase_Nsources(ph) + sourceState(ph)%p(so)%subState0(:,me) = sourceState(ph)%p(so)%state(:,me) + enddo + endif +!-------------------------------------------------------------------------------------------------- +! cut back (reduced time and restore) + else + subStep = num%subStepSizeCryst * subStep + constitutive_mech_Fp(ph)%data(1:3,1:3,me) = crystallite_subFp0(1:3,1:3,co,ip,el) + constitutive_mech_Fi(ph)%data(1:3,1:3,me) = crystallite_subFi0(1:3,1:3,co,ip,el) + crystallite_S (1:3,1:3,co,ip,el) = crystallite_S0 (1:3,1:3,co,ip,el) + if (subStep < 1.0_pReal) then ! actual (not initial) cutback + crystallite_Lp (1:3,1:3,co,ip,el) = subLp0 + constitutive_mech_Li(ph)%data(1:3,1:3,me) = subLi0 + endif + plasticState(ph)%state(:,me) = plasticState(ph)%subState0(:,me) + do so = 1, phase_Nsources(ph) + sourceState(ph)%p(so)%state(:,me) = sourceState(ph)%p(so)%subState0(:,me) + enddo + + todo = subStep > num%subStepMinCryst ! still on track or already done (beyond repair) + endif + +!-------------------------------------------------------------------------------------------------- +! prepare for integration + if (todo) then + subF = subF0 & + + subStep * (crystallite_F(1:3,1:3,co,ip,el) - crystallite_partitionedF0(1:3,1:3,co,ip,el)) + crystallite_Fe(1:3,1:3,co,ip,el) = matmul(subF,math_inv33(matmul(constitutive_mech_Fi(ph)%data(1:3,1:3,me), & + constitutive_mech_Fp(ph)%data(1:3,1:3,me)))) + crystallite_subdt(co,ip,el) = subStep * dt + converged_ = .not. integrateState(subF0,subF,subStep * dt,co,ip,el) + converged_ = converged_ .and. .not. integrateSourceState(subStep * dt,co,ip,el) + endif + + enddo cutbackLooping + +end function crystallite_stress + + +!-------------------------------------------------------------------------------------------------- +!> @brief Restore data after homog cutback. +!-------------------------------------------------------------------------------------------------- +module subroutine mech_restore(ip,el,includeL) + + integer, intent(in) :: & + ip, & !< integration point number + el !< element number + logical, intent(in) :: & + includeL !< protect agains fake cutback + integer :: & + co, p, m !< constituent number + + do co = 1,homogenization_Nconstituents(material_homogenizationAt(el)) + p = material_phaseAt(co,el) + m = material_phaseMemberAt(co,ip,el) + if (includeL) then + crystallite_Lp(1:3,1:3,co,ip,el) = crystallite_partitionedLp0(1:3,1:3,co,ip,el) + constitutive_mech_Li(p)%data(1:3,1:3,m) = constitutive_mech_partitionedLi0(p)%data(1:3,1:3,m) + endif ! maybe protecting everything from overwriting makes more sense + + constitutive_mech_Fp(p)%data(1:3,1:3,m) = constitutive_mech_partitionedFp0(p)%data(1:3,1:3,m) + constitutive_mech_Fi(p)%data(1:3,1:3,m) = constitutive_mech_partitionedFi0(p)%data(1:3,1:3,m) + crystallite_S (1:3,1:3,co,ip,el) = crystallite_partitionedS0 (1:3,1:3,co,ip,el) + + plasticState (material_phaseAt(co,el))%state( :,material_phasememberAt(co,ip,el)) = & + plasticState (material_phaseAt(co,el))%partitionedState0(:,material_phasememberAt(co,ip,el)) + enddo + +end subroutine mech_restore end submodule constitutive_mech diff --git a/src/constitutive_plastic_dislotwin.f90 b/src/constitutive_plastic_dislotwin.f90 index 4234a55b8..0474427fe 100644 --- a/src/constitutive_plastic_dislotwin.f90 +++ b/src/constitutive_plastic_dislotwin.f90 @@ -485,12 +485,12 @@ end function plastic_dislotwin_init !-------------------------------------------------------------------------------------------------- !> @brief Return the homogenized elasticity matrix. !-------------------------------------------------------------------------------------------------- -module function plastic_dislotwin_homogenizedC(ipc,ip,el) result(homogenizedC) +module function plastic_dislotwin_homogenizedC(co,ip,el) result(homogenizedC) real(pReal), dimension(6,6) :: & homogenizedC integer, intent(in) :: & - ipc, & !< component-ID of integration point + co, & !< component-ID of integration point ip, & !< integration point el !< element @@ -498,9 +498,9 @@ module function plastic_dislotwin_homogenizedC(ipc,ip,el) result(homogenizedC) of real(pReal) :: f_unrotated - of = material_phasememberAt(ipc,ip,el) - associate(prm => param(phase_plasticityInstance(material_phaseAt(ipc,el))),& - stt => state(phase_plasticityInstance(material_phaseAT(ipc,el)))) + of = material_phasememberAt(co,ip,el) + associate(prm => param(phase_plasticityInstance(material_phaseAt(co,el))),& + stt => state(phase_plasticityInstance(material_phaseAT(co,el)))) f_unrotated = 1.0_pReal & - sum(stt%f_tw(1:prm%sum_N_tw,of)) & diff --git a/src/constitutive_plastic_nonlocal.f90 b/src/constitutive_plastic_nonlocal.f90 index 65f74b6cc..0d7875291 100644 --- a/src/constitutive_plastic_nonlocal.f90 +++ b/src/constitutive_plastic_nonlocal.f90 @@ -552,11 +552,10 @@ end function plastic_nonlocal_init !-------------------------------------------------------------------------------------------------- !> @brief calculates quantities characterizing the microstructure !-------------------------------------------------------------------------------------------------- -module subroutine plastic_nonlocal_dependentState(F, Fp, instance, of, ip, el) +module subroutine plastic_nonlocal_dependentState(F, instance, of, ip, el) real(pReal), dimension(3,3), intent(in) :: & - F, & - Fp + F integer, intent(in) :: & instance, & of, & @@ -564,6 +563,8 @@ module subroutine plastic_nonlocal_dependentState(F, Fp, instance, of, ip, el) el integer :: & + ph, & + me, & no, & !< neighbor offset neighbor_el, & ! element number of neighboring material point neighbor_ip, & ! integration point of neighboring material point @@ -643,8 +644,10 @@ module subroutine plastic_nonlocal_dependentState(F, Fp, instance, of, ip, el) rho0 = getRho0(instance,of,ip,el) if (.not. phase_localPlasticity(material_phaseAt(1,el)) .and. prm%shortRangeStressCorrection) then - invFp = math_inv33(Fp) - invFe = matmul(Fp,math_inv33(F)) + ph = material_phaseAt(1,el) + me = material_phaseMemberAt(1,ip,el) + invFp = math_inv33(constitutive_mech_Fp(ph)%data(1:3,1:3,me)) + invFe = matmul(constitutive_mech_Fp(ph)%data(1:3,1:3,me),math_inv33(F)) rho_edg_delta = rho0(:,mob_edg_pos) - rho0(:,mob_edg_neg) rho_scr_delta = rho0(:,mob_scr_pos) - rho0(:,mob_scr_neg) @@ -973,14 +976,13 @@ end subroutine plastic_nonlocal_deltaState !--------------------------------------------------------------------------------------------------- !> @brief calculates the rate of change of microstructure !--------------------------------------------------------------------------------------------------- -module subroutine plastic_nonlocal_dotState(Mp, F, Fp, Temperature,timestep, & +module subroutine plastic_nonlocal_dotState(Mp, F, Temperature,timestep, & instance,of,ip,el) real(pReal), dimension(3,3), intent(in) :: & Mp !< MandelStress real(pReal), dimension(3,3,homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems), intent(in) :: & - F, & !< elastic deformation gradient - Fp !< plastic deformation gradient + F !< Deformation gradient real(pReal), intent(in) :: & Temperature, & !< temperature timestep !< substepped crystallite time increment @@ -1147,7 +1149,7 @@ module subroutine plastic_nonlocal_dotState(Mp, F, Fp, Temperature,timestep, & - rhoDip(s,1) / timestep - rhoDotAthermalAnnihilation(s,9) & - rhoDotSingle2DipoleGlide(s,9)) ! make sure that we do not annihilate more dipoles than we have - rhoDot = rhoDotFlux(F,Fp,timestep, instance,of,ip,el) & + rhoDot = rhoDotFlux(F,timestep, instance,of,ip,el) & + rhoDotMultiplication & + rhoDotSingle2DipoleGlide & + rhoDotAthermalAnnihilation & @@ -1176,11 +1178,10 @@ end subroutine plastic_nonlocal_dotState !--------------------------------------------------------------------------------------------------- !> @brief calculates the rate of change of microstructure !--------------------------------------------------------------------------------------------------- -function rhoDotFlux(F,Fp,timestep, instance,of,ip,el) +function rhoDotFlux(F,timestep, instance,of,ip,el) real(pReal), dimension(3,3,homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems), intent(in) :: & - F, & !< elastic deformation gradient - Fp !< plastic deformation gradient + F !< Deformation gradient real(pReal), intent(in) :: & timestep !< substepped crystallite time increment integer, intent(in) :: & @@ -1293,7 +1294,7 @@ function rhoDotFlux(F,Fp,timestep, instance,of,ip,el) m(1:3,:,4) = prm%slip_transverse my_F = F(1:3,1:3,1,ip,el) - my_Fe = matmul(my_F, math_inv33(Fp(1:3,1:3,1,ip,el))) + my_Fe = matmul(my_F, math_inv33(constitutive_mech_Fp(ph)%data(1:3,1:3,of))) neighbors: do n = 1,nIPneighbors @@ -1311,7 +1312,7 @@ function rhoDotFlux(F,Fp,timestep, instance,of,ip,el) if (neighbor_n > 0) then ! if neighbor exists, average deformation gradient neighbor_instance = phase_plasticityInstance(material_phaseAt(1,neighbor_el)) neighbor_F = F(1:3,1:3,1,neighbor_ip,neighbor_el) - neighbor_Fe = matmul(neighbor_F, math_inv33(Fp(1:3,1:3,1,neighbor_ip,neighbor_el))) + neighbor_Fe = matmul(neighbor_F, math_inv33(constitutive_mech_Fp(np)%data(1:3,1:3,no))) Favg = 0.5_pReal * (my_F + neighbor_F) else ! if no neighbor, take my value as average Favg = my_F diff --git a/src/damage_nonlocal.f90 b/src/damage_nonlocal.f90 index ac4d8636a..3db63cab2 100644 --- a/src/damage_nonlocal.f90 +++ b/src/damage_nonlocal.f90 @@ -148,12 +148,12 @@ real(pReal) function damage_nonlocal_getMobility(ip,el) ip, & !< integration point number el !< element number integer :: & - ipc + co damage_nonlocal_getMobility = 0.0_pReal - do ipc = 1, homogenization_Nconstituents(material_homogenizationAt(el)) - damage_nonlocal_getMobility = damage_nonlocal_getMobility + lattice_M(material_phaseAt(ipc,el)) + do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) + damage_nonlocal_getMobility = damage_nonlocal_getMobility + lattice_M(material_phaseAt(co,el)) enddo damage_nonlocal_getMobility = damage_nonlocal_getMobility/& diff --git a/src/grid/discretization_grid.f90 b/src/grid/discretization_grid.f90 index 1b3700c14..48ad5b7e1 100644 --- a/src/grid/discretization_grid.f90 +++ b/src/grid/discretization_grid.f90 @@ -19,7 +19,6 @@ module discretization_grid use results use discretization use geometry_plastic_nonlocal - use FEsolving implicit none private @@ -117,9 +116,6 @@ subroutine discretization_grid_init(restart) (grid(1)+1) * (grid(2)+1) * grid3,& ! ...unless not last process worldrank+1==worldsize)) - FEsolving_execElem = [1,product(myGrid)] ! parallel loop bounds set to comprise all elements - FEsolving_execIP = [1,1] ! parallel loop bounds set to comprise the only IP - !-------------------------------------------------------------------------------------------------- ! store geometry information for post processing if(.not. restart) then diff --git a/src/grid/grid_mech_FEM.f90 b/src/grid/grid_mech_FEM.f90 index cdf806b35..003f568c6 100644 --- a/src/grid/grid_mech_FEM.f90 +++ b/src/grid/grid_mech_FEM.f90 @@ -18,7 +18,6 @@ module grid_mech_FEM use math use rotations use spectral_utilities - use FEsolving use config use homogenization use discretization diff --git a/src/grid/grid_mech_spectral_basic.f90 b/src/grid/grid_mech_spectral_basic.f90 index ebaaf3b55..9bc36165f 100644 --- a/src/grid/grid_mech_spectral_basic.f90 +++ b/src/grid/grid_mech_spectral_basic.f90 @@ -18,7 +18,6 @@ module grid_mech_spectral_basic use math use rotations use spectral_utilities - use FEsolving use config use homogenization use discretization_grid diff --git a/src/grid/grid_mech_spectral_polarisation.f90 b/src/grid/grid_mech_spectral_polarisation.f90 index 9f2a17c97..7160c1adc 100644 --- a/src/grid/grid_mech_spectral_polarisation.f90 +++ b/src/grid/grid_mech_spectral_polarisation.f90 @@ -18,7 +18,6 @@ module grid_mech_spectral_polarisation use math use rotations use spectral_utilities - use FEsolving use config use homogenization use discretization_grid diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index c0c84233d..e8bae223a 100644 --- a/src/grid/spectral_utilities.f90 +++ b/src/grid/spectral_utilities.f90 @@ -810,9 +810,9 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,& print'(/,a)', ' ... evaluating constitutive response ......................................' flush(IO_STDOUT) - homogenization_F = reshape(F,[3,3,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 + call materialpoint_stressAndItsTangent(timeinc,[1,1],[1,product(grid(1:2))*grid3]) ! calculate P field P = reshape(homogenization_P, [3,3,grid(1),grid(2),grid3]) P_av = sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt ! average of P diff --git a/src/homogenization.f90 b/src/homogenization.f90 index cbab8e468..52553b57b 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -11,7 +11,6 @@ module homogenization use math use material use constitutive - use FEsolving use discretization use thermal_isothermal use thermal_conduction @@ -48,20 +47,6 @@ module homogenization type(tNumerics) :: num - type :: tDebugOptions - logical :: & - basic, & - extensive, & - selective - integer :: & - element, & - ip, & - grain - end type tDebugOptions - - type(tDebugOptions) :: debugHomog - - !-------------------------------------------------------------------------------------------------- interface @@ -85,33 +70,27 @@ module homogenization end subroutine mech_homogenize module subroutine mech_results(group_base,h) - character(len=*), intent(in) :: group_base integer, intent(in) :: h - end subroutine mech_results -! -------- ToDo --------------------------------------------------------- - module function mech_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) - logical, dimension(2) :: mech_RGC_updateState - real(pReal), dimension(:,:,:), intent(in) :: & - P,& !< partitioned stresses - F,& !< partitioned deformation gradients - F0 !< partitioned initial deformation gradients - real(pReal), dimension(:,:,:,:,:), intent(in) :: dPdF !< partitioned stiffnesses - real(pReal), dimension(3,3), intent(in) :: avgF !< average F - real(pReal), intent(in) :: dt !< time increment - integer, intent(in) :: & - ip, & !< integration point number - el !< element number - end function mech_RGC_updateState + module function mech_updateState(subdt,subF,ip,el) result(doneAndHappy) + real(pReal), intent(in) :: & + subdt !< current time step + real(pReal), intent(in), dimension(3,3) :: & + subF + integer, intent(in) :: & + ip, & !< integration point + el !< element number + logical, dimension(2) :: doneAndHappy + end function mech_updateState end interface -! ----------------------------------------------------------------------- public :: & homogenization_init, & materialpoint_stressAndItsTangent, & + homogenization_forward, & homogenization_results contains @@ -124,24 +103,10 @@ subroutine homogenization_init class (tNode) , pointer :: & num_homog, & - num_homogGeneric, & - debug_homogenization + num_homogGeneric print'(/,a)', ' <<<+- homogenization init -+>>>'; flush(IO_STDOUT) - debug_homogenization => config_debug%get('homogenization', defaultVal=emptyList) - debugHomog%basic = debug_homogenization%contains('basic') - debugHomog%extensive = debug_homogenization%contains('extensive') - debugHomog%selective = debug_homogenization%contains('selective') - debugHomog%element = config_debug%get_asInt('element',defaultVal = 1) - debugHomog%ip = config_debug%get_asInt('integrationpoint',defaultVal = 1) - debugHomog%grain = config_debug%get_asInt('grain',defaultVal = 1) - - if (debugHomog%grain < 1 & - .or. debugHomog%grain > homogenization_Nconstituents(material_homogenizationAt(debugHomog%element))) & - call IO_error(602,ext_msg='constituent', el=debugHomog%element, g=debugHomog%grain) - - num_homog => config_numerics%get('homogenization',defaultVal=emptyDict) num_homogGeneric => num_homog%get('generic',defaultVal=emptyDict) @@ -171,178 +136,130 @@ end subroutine homogenization_init !-------------------------------------------------------------------------------------------------- !> @brief parallelized calculation of stress and corresponding tangent at material points !-------------------------------------------------------------------------------------------------- -subroutine materialpoint_stressAndItsTangent(dt) +subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execElem) real(pReal), intent(in) :: dt !< time increment + integer, dimension(2), intent(in) :: FEsolving_execElem, FEsolving_execIP integer :: & - NiterationHomog, & NiterationMPstate, & - i, & !< integration point number - e, & !< element number - myNgrains - real(pReal), dimension(discretization_nIPs,discretization_Nelems) :: & + ip, & !< integration point number + el, & !< element number + myNgrains, co, ce, ho, me + real(pReal) :: & subFrac, & subStep - logical, dimension(discretization_nIPs,discretization_Nelems) :: & - requested, & + logical :: & converged - logical, dimension(2,discretization_nIPs,discretization_Nelems) :: & + logical, dimension(2) :: & doneAndHappy - integer :: m + !$OMP PARALLEL DO PRIVATE(ce,me,ho,myNgrains,NiterationMPstate,subFrac,converged,subStep,doneAndHappy) + do el = FEsolving_execElem(1),FEsolving_execElem(2) + ho = material_homogenizationAt(el) + myNgrains = homogenization_Nconstituents(ho) + do ip = FEsolving_execIP(1),FEsolving_execIP(2) + me = material_homogenizationMemberAt(ip,el) !-------------------------------------------------------------------------------------------------- ! initialize restoration points - do e = FEsolving_execElem(1),FEsolving_execElem(2) - do i = FEsolving_execIP(1),FEsolving_execIP(2); + call constitutive_initializeRestorationPoints(ip,el) - call crystallite_initializeRestorationPoints(i,e) + subFrac = 0.0_pReal + converged = .false. ! pretend failed step ... + subStep = 1.0_pReal/num%subStepSizeHomog ! ... larger then the requested calculation - subFrac(i,e) = 0.0_pReal - converged(i,e) = .false. ! pretend failed step ... - subStep(i,e) = 1.0_pReal/num%subStepSizeHomog ! ... larger then the requested calculation - requested(i,e) = .true. ! everybody requires calculation + if (homogState(ho)%sizeState > 0) & + homogState(ho)%subState0(:,me) = homogState(ho)%State0(:,me) + if (damageState(ho)%sizeState > 0) & + damageState(ho)%subState0(:,me) = damageState(ho)%State0(:,me) - if (homogState(material_homogenizationAt(e))%sizeState > 0) & - homogState(material_homogenizationAt(e))%subState0(:,material_homogenizationMemberAt(i,e)) = & - homogState(material_homogenizationAt(e))%State0( :,material_homogenizationMemberAt(i,e)) + cutBackLooping: do while (.not. terminallyIll .and. subStep > num%subStepMinHomog) - if (damageState(material_homogenizationAt(e))%sizeState > 0) & - damageState(material_homogenizationAt(e))%subState0(:,material_homogenizationMemberAt(i,e)) = & - damageState(material_homogenizationAt(e))%State0( :,material_homogenizationMemberAt(i,e)) - enddo - enddo + if (converged) then + subFrac = subFrac + subStep + subStep = min(1.0_pReal-subFrac,num%stepIncreaseHomog*subStep) ! introduce flexibility for step increase/acceleration - NiterationHomog = 0 - - cutBackLooping: do while (.not. terminallyIll .and. & - any(subStep(FEsolving_execIP(1):FEsolving_execIP(2),& - FEsolving_execElem(1):FEsolving_execElem(2)) > num%subStepMinHomog)) - - !$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) - - if (converged(i,e)) then - subFrac(i,e) = subFrac(i,e) + subStep(i,e) - subStep(i,e) = min(1.0_pReal-subFrac(i,e),num%stepIncreaseHomog*subStep(i,e)) ! introduce flexibility for step increase/acceleration - - steppingNeeded: if (subStep(i,e) > num%subStepMinHomog) then + steppingNeeded: if (subStep > num%subStepMinHomog) then ! wind forward grain starting point - call crystallite_windForward(i,e) + call constitutive_windForward(ip,el) - if(homogState(material_homogenizationAt(e))%sizeState > 0) & - homogState(material_homogenizationAt(e))%subState0(:,material_homogenizationMemberAt(i,e)) = & - homogState(material_homogenizationAt(e))%State (:,material_homogenizationMemberAt(i,e)) - if(damageState(material_homogenizationAt(e))%sizeState > 0) & - damageState(material_homogenizationAt(e))%subState0(:,material_homogenizationMemberAt(i,e)) = & - damageState(material_homogenizationAt(e))%State (:,material_homogenizationMemberAt(i,e)) + if(homogState(ho)%sizeState > 0) & + homogState(ho)%subState0(:,me) = homogState(ho)%State(:,me) + if(damageState(ho)%sizeState > 0) & + damageState(ho)%subState0(:,me) = damageState(ho)%State(:,me) endif steppingNeeded - - else - if ( (myNgrains == 1 .and. subStep(i,e) <= 1.0 ) .or. & ! single grain already tried internal subStepping in crystallite - num%subStepSizeHomog * subStep(i,e) <= num%subStepMinHomog ) then ! would require too small subStep + elseif ( (myNgrains == 1 .and. subStep <= 1.0 ) .or. & ! single grain already tried internal subStepping in crystallite + num%subStepSizeHomog * subStep <= num%subStepMinHomog ) then ! would require too small subStep ! cutback makes no sense - if (.not. terminallyIll) then ! so first signals terminally ill... - print*, ' Integration point ', i,' at element ', e, ' terminally ill' - endif - terminallyIll = .true. ! ...and kills all others - else ! cutback makes sense - subStep(i,e) = num%subStepSizeHomog * subStep(i,e) ! crystallite had severe trouble, so do a significant cutback + if (.not. terminallyIll) & ! so first signals terminally ill... + print*, ' Integration point ', ip,' at element ', el, ' terminally ill' + terminallyIll = .true. ! ...and kills all others + else ! cutback makes sense + subStep = num%subStepSizeHomog * subStep ! crystallite had severe trouble, so do a significant cutback - call crystallite_restore(i,e,subStep(i,e) < 1.0_pReal) - call constitutive_restore(i,e) + call constitutive_restore(ip,el,subStep < 1.0_pReal) - if(homogState(material_homogenizationAt(e))%sizeState > 0) & - homogState(material_homogenizationAt(e))%State( :,material_homogenizationMemberAt(i,e)) = & - homogState(material_homogenizationAt(e))%subState0(:,material_homogenizationMemberAt(i,e)) - if(damageState(material_homogenizationAt(e))%sizeState > 0) & - damageState(material_homogenizationAt(e))%State( :,material_homogenizationMemberAt(i,e)) = & - damageState(material_homogenizationAt(e))%subState0(:,material_homogenizationMemberAt(i,e)) - endif + if(homogState(ho)%sizeState > 0) & + homogState(ho)%State(:,me) = homogState(ho)%subState0(:,me) + if(damageState(ho)%sizeState > 0) & + damageState(ho)%State(:,me) = damageState(ho)%subState0(:,me) endif - if (subStep(i,e) > num%subStepMinHomog) then - requested(i,e) = .true. - doneAndHappy(1:2,i,e) = [.false.,.true.] - endif - enddo IpLooping1 - enddo elementLooping1 - !$OMP END PARALLEL DO + if (subStep > num%subStepMinHomog) doneAndHappy = [.false.,.true.] - NiterationMPstate = 0 - - convergenceLooping: do while (.not. terminallyIll .and. & - any( requested(:,FEsolving_execELem(1):FEsolving_execElem(2)) & - .and. .not. doneAndHappy(1,:,FEsolving_execELem(1):FEsolving_execElem(2)) & - ) .and. & - NiterationMPstate < num%nMPstate) - NiterationMPstate = NiterationMPstate + 1 + NiterationMPstate = 0 + convergenceLooping: do while (.not. terminallyIll & + .and. .not. doneAndHappy(1) & + .and. NiterationMPstate < num%nMPstate) + NiterationMPstate = NiterationMPstate + 1 !-------------------------------------------------------------------------------------------------- ! deformation partitioning - !$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 - 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 - crystallite_requested(1:myNgrains,i,e) = .true. ! request calculation for constituents - else - crystallite_requested(1:myNgrains,i,e) = .false. ! calculation for constituents not required anymore - endif - enddo IpLooping2 - enddo elementLooping2 - !$OMP END PARALLEL DO -!-------------------------------------------------------------------------------------------------- -! crystallite integration - converged = crystallite_stress() !ToDo: MD not sure if that is the best logic + if (.not. doneAndHappy(1)) then + ce = (el-1)*discretization_nIPs + ip + call mech_partition(homogenization_F0(1:3,1:3,ce) & + + (homogenization_F(1:3,1:3,ce)-homogenization_F0(1:3,1:3,ce))& + *(subStep+subFrac), & + ip,el) + converged = .true. + do co = 1, myNgrains + converged = converged .and. crystallite_stress(dt*subStep,co,ip,el) + enddo -!-------------------------------------------------------------------------------------------------- -! state update - !$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.] + if (.not. converged) then + doneAndHappy = [.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,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 + ce = (el-1)*discretization_nIPs + ip + doneAndHappy = mech_updateState(dt*subStep, & + homogenization_F0(1:3,1:3,ce) & + + (homogenization_F(1:3,1:3,ce)-homogenization_F0(1:3,1:3,ce)) & + *(subStep+subFrac), & + ip,el) + converged = all(doneAndHappy) endif endif - enddo IpLooping3 - enddo elementLooping3 - !$OMP END PARALLEL DO - enddo convergenceLooping - - NiterationHomog = NiterationHomog + 1 - - enddo cutBackLooping + enddo convergenceLooping + enddo cutBackLooping + enddo + enddo + !$OMP END PARALLEL DO if (.not. terminallyIll ) then - call crystallite_orientations() ! calculate crystal orientations - !$OMP PARALLEL DO - elementLooping4: do e = FEsolving_execElem(1),FEsolving_execElem(2) - IpLooping4: do i = FEsolving_execIP(1),FEsolving_execIP(2) - call mech_homogenize(i,e) - enddo IpLooping4 - enddo elementLooping4 + !$OMP PARALLEL DO PRIVATE(ho,myNgrains) + elementLooping3: do el = FEsolving_execElem(1),FEsolving_execElem(2) + ho = material_homogenizationAt(el) + myNgrains = homogenization_Nconstituents(ho) + IpLooping3: do ip = FEsolving_execIP(1),FEsolving_execIP(2) + do co = 1, myNgrains + call crystallite_orientations(co,ip,el) + enddo + call mech_homogenize(ip,el) + enddo IpLooping3 + enddo elementLooping3 !$OMP END PARALLEL DO else print'(/,a,/)', ' << HOMOG >> Material Point terminally ill' @@ -351,77 +268,56 @@ subroutine materialpoint_stressAndItsTangent(dt) end subroutine materialpoint_stressAndItsTangent -!-------------------------------------------------------------------------------------------------- -!> @brief update the internal state of the homogenization scheme and tell whether "done" and -!> "happy" with result -!-------------------------------------------------------------------------------------------------- -function updateState(subdt,subF,ip,el) - - real(pReal), intent(in) :: & - subdt !< current time step - real(pReal), intent(in), dimension(3,3) :: & - subF - integer, intent(in) :: & - ip, & !< integration point - el !< element number - integer :: c - logical, dimension(2) :: updateState - real(pReal) :: dPdFs(3,3,3,3,homogenization_Nconstituents(material_homogenizationAt(el))) - - updateState = .true. - chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el))) - case (HOMOGENIZATION_RGC_ID) chosenHomogenization - do c=1,homogenization_Nconstituents(material_homogenizationAt(el)) - dPdFs(:,:,:,:,c) = crystallite_stressTangent(c,ip,el) - enddo - updateState = & - updateState .and. & - mech_RGC_updateState(crystallite_P(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), & - crystallite_partitionedF(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), & - crystallite_partitionedF0(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el),& - subF,& - subdt, & - dPdFs, & - ip, & - el) - end select chosenHomogenization - -end function updateState - - !-------------------------------------------------------------------------------------------------- !> @brief writes homogenization results to HDF5 output file !-------------------------------------------------------------------------------------------------- subroutine homogenization_results - use material, only: & - material_homogenization_type => homogenization_type - integer :: p + integer :: ho character(len=:), allocatable :: group_base,group - do p=1,size(material_name_homogenization) - group_base = 'current/homogenization/'//trim(material_name_homogenization(p)) + call results_closeGroup(results_addGroup('current/homogenization/')) + + do ho=1,size(material_name_homogenization) + group_base = 'current/homogenization/'//trim(material_name_homogenization(ho)) call results_closeGroup(results_addGroup(group_base)) - call mech_results(group_base,p) + call mech_results(group_base,ho) group = trim(group_base)//'/damage' call results_closeGroup(results_addGroup(group)) - select case(damage_type(p)) + select case(damage_type(ho)) case(DAMAGE_NONLOCAL_ID) - call damage_nonlocal_results(p,group) + call damage_nonlocal_results(ho,group) end select group = trim(group_base)//'/thermal' call results_closeGroup(results_addGroup(group)) - select case(thermal_type(p)) + select case(thermal_type(ho)) case(THERMAL_CONDUCTION_ID) - call thermal_conduction_results(p,group) + call thermal_conduction_results(ho,group) end select enddo end subroutine homogenization_results + +!-------------------------------------------------------------------------------------------------- +!> @brief Forward data after successful increment. +! ToDo: Any guessing for the current states possible? +!-------------------------------------------------------------------------------------------------- +subroutine homogenization_forward + + integer :: ho + + + do ho = 1, size(material_name_homogenization) + homogState (ho)%state0 = homogState (ho)%state + damageState(ho)%state0 = damageState(ho)%state + enddo + +end subroutine homogenization_forward + end module homogenization diff --git a/src/homogenization_mech.f90 b/src/homogenization_mech.f90 index b0641be07..641e960fd 100644 --- a/src/homogenization_mech.f90 +++ b/src/homogenization_mech.f90 @@ -4,6 +4,7 @@ !-------------------------------------------------------------------------------------------------- submodule(homogenization) homogenization_mech + interface module subroutine mech_none_init @@ -51,6 +52,21 @@ submodule(homogenization) homogenization_mech end subroutine mech_RGC_averageStressAndItsTangent + module function mech_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) result(doneAndHappy) + logical, dimension(2) :: doneAndHappy + real(pReal), dimension(:,:,:), intent(in) :: & + P,& !< partitioned stresses + F,& !< partitioned deformation gradients + F0 !< partitioned initial deformation gradients + real(pReal), dimension(:,:,:,:,:), intent(in) :: dPdF !< partitioned stiffnesses + real(pReal), dimension(3,3), intent(in) :: avgF !< average F + real(pReal), intent(in) :: dt !< time increment + integer, intent(in) :: & + ip, & !< integration point number + el !< element number + end function mech_RGC_updateState + + module subroutine mech_RGC_results(instance,group) integer, intent(in) :: instance !< homogenization instance character(len=*), intent(in) :: group !< group name in HDF5 file @@ -100,16 +116,16 @@ module subroutine mech_partition(subF,ip,el) chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el))) case (HOMOGENIZATION_NONE_ID) chosenHomogenization - crystallite_partitionedF(1:3,1:3,1,ip,el) = subF + crystallite_F(1:3,1:3,1,ip,el) = subF case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization call mech_isostrain_partitionDeformation(& - crystallite_partitionedF(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), & + crystallite_F(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), & subF) case (HOMOGENIZATION_RGC_ID) chosenHomogenization call mech_RGC_partitionDeformation(& - crystallite_partitionedF(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), & + crystallite_F(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), & subF,& ip, & el) @@ -127,35 +143,35 @@ module subroutine mech_homogenize(ip,el) integer, intent(in) :: & ip, & !< integration point el !< element number - integer :: c,m + integer :: co,ce real(pReal) :: dPdFs(3,3,3,3,homogenization_Nconstituents(material_homogenizationAt(el))) - m = (el-1)* discretization_nIPs + ip + ce = (el-1)* discretization_nIPs + ip chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el))) case (HOMOGENIZATION_NONE_ID) chosenHomogenization - 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) + homogenization_P(1:3,1:3,ce) = crystallite_P(1:3,1:3,1,ip,el) + homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = 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) + do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) + dPdFs(:,:,:,:,co) = crystallite_stressTangent(co,ip,el) enddo call mech_isostrain_averageStressAndItsTangent(& - homogenization_P(1:3,1:3,m), & - homogenization_dPdF(1:3,1:3,1:3,1:3,m),& + homogenization_P(1:3,1:3,ce), & + homogenization_dPdF(1:3,1:3,1:3,1:3,ce),& crystallite_P(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), & dPdFs, & homogenization_typeInstance(material_homogenizationAt(el))) case (HOMOGENIZATION_RGC_ID) chosenHomogenization - do c = 1, homogenization_Nconstituents(material_homogenizationAt(el)) - dPdFs(:,:,:,:,c) = crystallite_stressTangent(c,ip,el) + do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) + dPdFs(:,:,:,:,co) = crystallite_stressTangent(co,ip,el) enddo call mech_RGC_averageStressAndItsTangent(& - homogenization_P(1:3,1:3,m), & - homogenization_dPdF(1:3,1:3,1:3,1:3,m),& + homogenization_P(1:3,1:3,ce), & + homogenization_dPdF(1:3,1:3,1:3,1:3,ce),& crystallite_P(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), & dPdFs, & homogenization_typeInstance(material_homogenizationAt(el))) @@ -165,6 +181,45 @@ module subroutine mech_homogenize(ip,el) end subroutine mech_homogenize +!-------------------------------------------------------------------------------------------------- +!> @brief update the internal state of the homogenization scheme and tell whether "done" and +!> "happy" with result +!-------------------------------------------------------------------------------------------------- +module function mech_updateState(subdt,subF,ip,el) result(doneAndHappy) + + real(pReal), intent(in) :: & + subdt !< current time step + real(pReal), intent(in), dimension(3,3) :: & + subF + integer, intent(in) :: & + ip, & !< integration point + el !< element number + logical, dimension(2) :: doneAndHappy + + integer :: co + real(pReal) :: dPdFs(3,3,3,3,homogenization_Nconstituents(material_homogenizationAt(el))) + + + if (homogenization_type(material_homogenizationAt(el)) == HOMOGENIZATION_RGC_ID) then + do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) + dPdFs(:,:,:,:,co) = crystallite_stressTangent(co,ip,el) + enddo + doneAndHappy = & + mech_RGC_updateState(crystallite_P(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), & + crystallite_F(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), & + crystallite_partitionedF0(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el),& + subF,& + subdt, & + dPdFs, & + ip, & + el) + else + doneAndHappy = .true. + endif + +end function mech_updateState + + !-------------------------------------------------------------------------------------------------- !> @brief Write results to file. !-------------------------------------------------------------------------------------------------- diff --git a/src/homogenization_mech_RGC.f90 b/src/homogenization_mech_RGC.f90 index 0a9d0ac92..04ec73845 100644 --- a/src/homogenization_mech_RGC.f90 +++ b/src/homogenization_mech_RGC.f90 @@ -8,6 +8,7 @@ !-------------------------------------------------------------------------------------------------- submodule(homogenization:homogenization_mech) homogenization_mech_RGC use rotations + use lattice type :: tParameters integer, dimension(:), allocatable :: & @@ -18,8 +19,6 @@ submodule(homogenization:homogenization_mech) homogenization_mech_RGC real(pReal), dimension(:), allocatable :: & D_alpha, & a_g - integer :: & - of_debug = 0 character(len=pStringLen), allocatable, dimension(:) :: & output end type tParameters @@ -151,12 +150,6 @@ module subroutine mech_RGC_init(num_homogMech) st0 => state0(homogenization_typeInstance(h)), & dst => dependentState(homogenization_typeInstance(h))) -#ifdef DEBUG - if (h==material_homogenizationAt(debugHomog%element)) then - prm%of_debug = material_homogenizationMemberAt(debugHomog%ip,debugHomog%element) - endif -#endif - #if defined (__GFORTRAN__) prm%output = output_asStrings(homogMech) #else @@ -239,17 +232,6 @@ module subroutine mech_RGC_partitionDeformation(F,avgF,instance,of) F(i,j,iGrain) = F(i,j,iGrain) + aVect(i)*nVect(j) ! calculating deformation relaxations due to interface relaxation enddo F(1:3,1:3,iGrain) = F(1:3,1:3,iGrain) + avgF ! resulting relaxed deformation gradient - -#ifdef DEBUG - if (debugHomog%extensive) then - print'(a,i3)',' Deformation gradient of grain: ',iGrain - do i = 1,3 - print'(1x,3(e15.8,1x))',(F(i,j,iGrain), j = 1,3) - enddo - print*,' ' - flush(IO_STDOUT) - endif -#endif enddo end associate @@ -261,7 +243,18 @@ end subroutine mech_RGC_partitionDeformation !> @brief update the internal state of the homogenization scheme and tell whether "done" and ! "happy" with result !-------------------------------------------------------------------------------------------------- -module procedure mech_RGC_updateState +module function mech_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) result(doneAndHappy) + logical, dimension(2) :: doneAndHappy + real(pReal), dimension(:,:,:), intent(in) :: & + P,& !< partitioned stresses + F,& !< partitioned deformation gradients + F0 !< partitioned initial deformation gradients + real(pReal), dimension(:,:,:,:,:), intent(in) :: dPdF !< partitioned stiffnesses + real(pReal), dimension(3,3), intent(in) :: avgF !< average F + real(pReal), intent(in) :: dt !< time increment + integer, intent(in) :: & + ip, & !< integration point number + el !< element number integer, dimension(4) :: intFaceN,intFaceP,faceID integer, dimension(3) :: nGDim,iGr3N,iGr3P @@ -273,13 +266,9 @@ module procedure mech_RGC_updateState logical :: error real(pReal), dimension(:,:), allocatable :: tract,jmatrix,jnverse,smatrix,pmatrix,rmatrix real(pReal), dimension(:), allocatable :: resid,relax,p_relax,p_resid,drelax -#ifdef DEBUG - integer, dimension(3) :: stresLoc - integer, dimension(2) :: residLoc -#endif zeroTimeStep: if(dEq0(dt)) then - mech_RGC_updateState = .true. ! pretend everything is fine and return + doneAndHappy = .true. ! pretend everything is fine and return return endif zeroTimeStep @@ -303,16 +292,6 @@ module procedure mech_RGC_updateState relax = stt%relaxationVector(:,of) drelax = stt%relaxationVector(:,of) - st0%relaxationVector(:,of) -#ifdef DEBUG - if (debugHomog%extensive) then - print*, 'Obtained state: ' - do i = 1,size(stt%relaxationVector(:,of)) - print'(1x,2(e15.8,1x))', stt%relaxationVector(i,of) - enddo - print*,' ' - endif -#endif - !-------------------------------------------------------------------------------------------------- ! computing interface mismatch and stress penalty tensor for all interfaces of all grains call stressPenalty(R,NN,avgF,F,ip,el,instance,of) @@ -353,13 +332,6 @@ module procedure mech_RGC_updateState enddo enddo -#ifdef DEBUG - if (debugHomog%extensive) then - print'(a,i3)',' Traction at interface: ',iNum - print'(1x,3(e15.8,1x))',(tract(iNum,j), j = 1,3) - print*,' ' - endif -#endif enddo !-------------------------------------------------------------------------------------------------- @@ -367,29 +339,12 @@ module procedure mech_RGC_updateState stresMax = maxval(abs(P)) ! get the maximum of first Piola-Kirchhoff (material) stress residMax = maxval(abs(tract)) ! get the maximum of the residual -#ifdef DEBUG - if (debugHomog%extensive .and. prm%of_debug == of) then - stresLoc = maxloc(abs(P)) - residLoc = maxloc(abs(tract)) - print'(a,i2,1x,i4)',' RGC residual check ... ',ip,el - print'(a,e15.8,a,i3,a,i2,i2)', ' Max stress: ',stresMax, & - '@ grain ',stresLoc(3),' in component ',stresLoc(1),stresLoc(2) - print'(a,e15.8,a,i3,a,i2)',' Max residual: ',residMax, & - ' @ iface ',residLoc(1),' in direction ',residLoc(2) - flush(IO_STDOUT) - endif -#endif - - mech_RGC_updateState = .false. + doneAndHappy = .false. !-------------------------------------------------------------------------------------------------- ! If convergence reached => done and happy if (residMax < num%rtol*stresMax .or. residMax < num%atol) then - mech_RGC_updateState = .true. -#ifdef DEBUG - if (debugHomog%extensive .and. prm%of_debug == of) & - print*, '... done and happy'; flush(IO_STDOUT) -#endif + doneAndHappy = .true. !-------------------------------------------------------------------------------------------------- ! compute/update the state for postResult, i.e., all energy densities computed by time-integration @@ -406,41 +361,14 @@ module procedure mech_RGC_updateState dst%relaxationRate_avg(of) = sum(abs(drelax))/dt/real(3*nIntFaceTot,pReal) dst%relaxationRate_max(of) = maxval(abs(drelax))/dt -#ifdef DEBUG - if (debugHomog%extensive .and. prm%of_debug == of) then - print'(a,e15.8)', ' Constitutive work: ',stt%work(of) - print'(a,3(1x,e15.8))', ' Magnitude mismatch: ',dst%mismatch(1,of), & - dst%mismatch(2,of), & - dst%mismatch(3,of) - print'(a,e15.8)', ' Penalty energy: ', stt%penaltyEnergy(of) - print'(a,e15.8,/)', ' Volume discrepancy: ', dst%volumeDiscrepancy(of) - print'(a,e15.8)', ' Maximum relaxation rate: ', dst%relaxationRate_max(of) - print'(a,e15.8,/)', ' Average relaxation rate: ', dst%relaxationRate_avg(of) - flush(IO_STDOUT) - endif -#endif - return !-------------------------------------------------------------------------------------------------- ! if residual blows-up => done but unhappy elseif (residMax > num%relMax*stresMax .or. residMax > num%absMax) then ! try to restart when residual blows up exceeding maximum bound - mech_RGC_updateState = [.true.,.false.] ! with direct cut-back - -#ifdef DEBUG - if (debugHomog%extensive .and. prm%of_debug == of) & - print'(a,/)', ' ... broken'; flush(IO_STDOUT) -#endif - + doneAndHappy = [.true.,.false.] ! with direct cut-back return - - else ! proceed with computing the Jacobian and state update -#ifdef DEBUG - if (debugHomog%extensive .and. prm%of_debug == of) & - print'(a,/)', ' ... not yet done'; flush(IO_STDOUT) -#endif - - endif + endif !--------------------------------------------------------------------------------------------------- ! construct the global Jacobian matrix for updating the global relaxation vector array when @@ -492,17 +420,6 @@ module procedure mech_RGC_updateState enddo enddo -#ifdef DEBUG - if (debugHomog%extensive) then - print*, 'Jacobian matrix of stress' - do i = 1,3*nIntFaceTot - print'(1x,100(e11.4,1x))',(smatrix(i,j), j = 1,3*nIntFaceTot) - enddo - print*,' ' - flush(IO_STDOUT) - endif -#endif - !-------------------------------------------------------------------------------------------------- ! ... of the stress penalty tangent (mismatch penalty and volume penalty, computed using numerical ! perturbation method) "pmatrix" @@ -552,16 +469,6 @@ module procedure mech_RGC_updateState pmatrix(:,ipert) = p_resid/num%pPert enddo -#ifdef DEBUG - if (debugHomog%extensive) then - print*, 'Jacobian matrix of penalty' - do i = 1,3*nIntFaceTot - print'(1x,100(e11.4,1x))',(pmatrix(i,j), j = 1,3*nIntFaceTot) - enddo - print*,' ' - flush(IO_STDOUT) - endif -#endif !-------------------------------------------------------------------------------------------------- ! ... of the numerical viscosity traction "rmatrix" @@ -571,48 +478,16 @@ module procedure mech_RGC_updateState (abs(drelax(i))/(num%refRelaxRate*dt))**(num%viscPower - 1.0_pReal) ! only in the main diagonal term enddo -#ifdef DEBUG - if (debugHomog%extensive) then - print*, 'Jacobian matrix of penalty' - do i = 1,3*nIntFaceTot - print'(1x,100(e11.4,1x))',(rmatrix(i,j), j = 1,3*nIntFaceTot) - enddo - print*,' ' - flush(IO_STDOUT) - endif -#endif !-------------------------------------------------------------------------------------------------- ! The overall Jacobian matrix summarizing contributions of smatrix, pmatrix, rmatrix allocate(jmatrix(3*nIntFaceTot,3*nIntFaceTot)); jmatrix = smatrix + pmatrix + rmatrix -#ifdef DEBUG - if (debugHomog%extensive) then - print*, 'Jacobian matrix (total)' - do i = 1,3*nIntFaceTot - print'(1x,100(e11.4,1x))',(jmatrix(i,j), j = 1,3*nIntFaceTot) - enddo - print*,' ' - flush(IO_STDOUT) - endif -#endif - !-------------------------------------------------------------------------------------------------- ! computing the update of the state variable (relaxation vectors) using the Jacobian matrix allocate(jnverse(3*nIntFaceTot,3*nIntFaceTot),source=0.0_pReal) call math_invert(jnverse,error,jmatrix) -#ifdef DEBUG - if (debugHomog%extensive) then - print*, 'Jacobian inverse' - do i = 1,3*nIntFaceTot - print'(1x,100(e11.4,1x))',(jnverse(i,j), j = 1,3*nIntFaceTot) - enddo - print*,' ' - flush(IO_STDOUT) - endif -#endif - !-------------------------------------------------------------------------------------------------- ! calculate the state update (global relaxation vectors) for the next Newton-Raphson iteration drelax = 0.0_pReal @@ -621,7 +496,7 @@ module procedure mech_RGC_updateState enddo; enddo stt%relaxationVector(:,of) = relax + drelax ! Updateing the state variable for the next iteration if (any(abs(drelax) > num%maxdRelax)) then ! Forcing cutback when the incremental change of relaxation vector becomes too large - mech_RGC_updateState = [.true.,.false.] + doneAndHappy = [.true.,.false.] !$OMP CRITICAL (write2out) print'(a,i3,a,i3,a)',' RGC_updateState: ip ',ip,' | el ',el,' enforces cutback' print'(a,e15.8)',' due to large relaxation change = ',maxval(abs(drelax)) @@ -629,17 +504,6 @@ module procedure mech_RGC_updateState !$OMP END CRITICAL (write2out) endif -#ifdef DEBUG - if (debugHomog%extensive) then - print*, 'Returned state: ' - do i = 1,size(stt%relaxationVector(:,of)) - print'(1x,2(e15.8,1x))', stt%relaxationVector(i,of) - enddo - print*,' ' - flush(IO_STDOUT) - endif -#endif - end associate contains @@ -661,8 +525,10 @@ module procedure mech_RGC_updateState real(pReal), dimension (3) :: nVect,surfCorr real(pReal), dimension (2) :: Gmoduli integer :: iGrain,iGNghb,iFace,i,j,k,l - real(pReal) :: muGrain,muGNghb,nDefNorm,bgGrain,bgGNghb - real(pReal), parameter :: nDefToler = 1.0e-10_pReal + real(pReal) :: muGrain,muGNghb,nDefNorm + real(pReal), parameter :: & + nDefToler = 1.0e-10_pReal, & + b = 2.5e-10_pReal ! Length of Burgers vector nGDim = param(instance)%N_constituents rPen = 0.0_pReal @@ -676,19 +542,11 @@ module procedure mech_RGC_updateState associate(prm => param(instance)) -#ifdef DEBUG - if (debugHomog%extensive .and. prm%of_debug == of) then - print'(a,2(1x,i3))', ' Correction factor: ',ip,el - print*, surfCorr - endif -#endif !----------------------------------------------------------------------------------------------- ! computing the mismatch and penalty stress tensor of all grains grainLoop: do iGrain = 1,product(prm%N_constituents) - Gmoduli = equivalentModuli(iGrain,ip,el) - muGrain = Gmoduli(1) ! collecting the equivalent shear modulus of grain - bgGrain = Gmoduli(2) ! and the lengthh of Burgers vector + muGrain = equivalentMu(iGrain,ip,el) iGrain3 = grain1to3(iGrain,prm%N_constituents) ! get the grain ID in local 3-dimensional index (x,y,z)-position interfaceLoop: do iFace = 1,6 @@ -700,9 +558,7 @@ module procedure mech_RGC_updateState where(iGNghb3 < 1) iGNghb3 = nGDim where(iGNghb3 >nGDim) iGNghb3 = 1 iGNghb = grain3to1(iGNghb3,prm%N_constituents) ! get the ID of the neighboring grain - Gmoduli = equivalentModuli(iGNghb,ip,el) ! collect the shear modulus and Burgers vector of the neighbor - muGNghb = Gmoduli(1) - bgGNghb = Gmoduli(2) + muGNghb = equivalentMu(iGNghb,ip,el) gDef = 0.5_pReal*(fDef(1:3,1:3,iGNghb) - fDef(1:3,1:3,iGrain)) ! difference/jump in deformation gradeint across the neighbor !------------------------------------------------------------------------------------------- @@ -717,30 +573,19 @@ module procedure mech_RGC_updateState enddo; enddo nDefNorm = max(nDefToler,sqrt(nDefNorm)) ! approximation to zero mismatch if mismatch is zero (singularity) nMis(abs(intFace(1)),iGrain) = nMis(abs(intFace(1)),iGrain) + nDefNorm ! total amount of mismatch experienced by the grain (at all six interfaces) -#ifdef DEBUG - if (debugHomog%extensive .and. prm%of_debug == of) then - print'(a,i2,a,i3)',' Mismatch to face: ',intFace(1),' neighbor grain: ',iGNghb - print*, transpose(nDef) - print'(a,e11.4)', ' with magnitude: ',nDefNorm - endif -#endif + !------------------------------------------------------------------------------------------- ! compute the stress penalty of all interfaces do i = 1,3; do j = 1,3; do k = 1,3; do l = 1,3 - rPen(i,j,iGrain) = rPen(i,j,iGrain) + 0.5_pReal*(muGrain*bgGrain + muGNghb*bgGNghb)*prm%xi_alpha & + rPen(i,j,iGrain) = rPen(i,j,iGrain) + 0.5_pReal*(muGrain*b + muGNghb*b)*prm%xi_alpha & *surfCorr(abs(intFace(1)))/prm%D_alpha(abs(intFace(1))) & *cosh(prm%c_alpha*nDefNorm) & *0.5_pReal*nVect(l)*nDef(i,k)/nDefNorm*math_LeviCivita(k,l,j) & *tanh(nDefNorm/num%xSmoo) enddo; enddo;enddo; enddo enddo interfaceLoop -#ifdef DEBUG - if (debugHomog%extensive .and. prm%of_debug == of) then - print'(a,i2)', ' Penalty of grain: ',iGrain - print*, transpose(rPen(1:3,1:3,iGrain)) - endif -#endif + enddo grainLoop @@ -783,13 +628,6 @@ module procedure mech_RGC_updateState vPen(:,:,i) = -1.0_pReal/real(nGrain,pReal)*num%volDiscrMod*num%volDiscrPow/num%maxVolDiscr* & sign((abs(vDiscrep)/num%maxVolDiscr)**(num%volDiscrPow - 1.0),vDiscrep)* & gVol(i)*transpose(math_inv33(fDef(:,:,i))) - -#ifdef DEBUG - if (debugHomog%extensive .and. param(instance)%of_debug == of) then - print'(a,i2)',' Volume penalty of grain: ',i - print*, transpose(vPen(:,:,i)) - endif -#endif enddo end subroutine volumePenalty @@ -827,44 +665,26 @@ module procedure mech_RGC_updateState end function surfaceCorrection - !-------------------------------------------------------------------------------------------------- + !------------------------------------------------------------------------------------------------- !> @brief compute the equivalent shear and bulk moduli from the elasticity tensor - !-------------------------------------------------------------------------------------------------- - function equivalentModuli(grainID,ip,el) - - real(pReal), dimension(2) :: equivalentModuli + !------------------------------------------------------------------------------------------------- + real(pReal) function equivalentMu(grainID,ip,el) integer, intent(in) :: & grainID,& ip, & !< integration point number el !< element number - real(pReal), dimension(6,6) :: elasTens - real(pReal) :: & - cEquiv_11, & - cEquiv_12, & - cEquiv_44 - - elasTens = constitutive_homogenizedC(grainID,ip,el) - - !---------------------------------------------------------------------------------------------- - ! compute the equivalent shear modulus after Turterltaub and Suiker, JMPS (2005) - cEquiv_11 = (elasTens(1,1) + elasTens(2,2) + elasTens(3,3))/3.0_pReal - cEquiv_12 = (elasTens(1,2) + elasTens(2,3) + elasTens(3,1) + & - elasTens(1,3) + elasTens(2,1) + elasTens(3,2))/6.0_pReal - cEquiv_44 = (elasTens(4,4) + elasTens(5,5) + elasTens(6,6))/3.0_pReal - equivalentModuli(1) = 0.2_pReal*(cEquiv_11 - cEquiv_12) + 0.6_pReal*cEquiv_44 - - !---------------------------------------------------------------------------------------------- - ! obtain the length of Burgers vector (could be model dependend) - equivalentModuli(2) = 2.5e-10_pReal - - end function equivalentModuli - !-------------------------------------------------------------------------------------------------- + equivalentMu = lattice_equivalent_mu(constitutive_homogenizedC(grainID,ip,el),'voigt') + + end function equivalentMu + + + !------------------------------------------------------------------------------------------------- !> @brief calculating the grain deformation gradient (the same with ! homogenization_RGC_partitionDeformation, but used only for perturbation scheme) - !-------------------------------------------------------------------------------------------------- + !------------------------------------------------------------------------------------------------- subroutine grainDeformation(F, avgF, instance, of) real(pReal), dimension(:,:,:), intent(out) :: F !< partitioned F per grain @@ -879,7 +699,7 @@ module procedure mech_RGC_updateState integer, dimension(3) :: iGrain3 integer :: iGrain,iFace,i,j - !------------------------------------------------------------------------------------------------- + !----------------------------------------------------------------------------------------------- ! compute the deformation gradient of individual grains due to relaxations associate(prm => param(instance)) @@ -901,7 +721,7 @@ module procedure mech_RGC_updateState end subroutine grainDeformation -end procedure mech_RGC_updateState +end function mech_RGC_updateState !-------------------------------------------------------------------------------------------------- diff --git a/src/kinematics_cleavage_opening.f90 b/src/kinematics_cleavage_opening.f90 index 6fe8ed7f6..a29a290f8 100644 --- a/src/kinematics_cleavage_opening.f90 +++ b/src/kinematics_cleavage_opening.f90 @@ -99,10 +99,10 @@ end function kinematics_cleavage_opening_init !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the velocity gradient !-------------------------------------------------------------------------------------------------- -module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, ip, el) +module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, co, ip, el) integer, intent(in) :: & - ipc, & !< grain number + co, & !< grain number ip, & !< integration point number el !< element number real(pReal), intent(in), dimension(3,3) :: & @@ -124,7 +124,7 @@ module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, Ld = 0.0_pReal dLd_dTstar = 0.0_pReal - associate(prm => param(kinematics_cleavage_opening_instance(material_phaseAt(ipc,el)))) + associate(prm => param(kinematics_cleavage_opening_instance(material_phaseAt(co,el)))) do i = 1,prm%sum_N_cl traction_crit = prm%g_crit(i)* damage(homog)%p(damageOffset)**2.0_pReal diff --git a/src/kinematics_slipplane_opening.f90 b/src/kinematics_slipplane_opening.f90 index b7adb6807..84edab122 100644 --- a/src/kinematics_slipplane_opening.f90 +++ b/src/kinematics_slipplane_opening.f90 @@ -117,10 +117,10 @@ end function kinematics_slipplane_opening_init !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the velocity gradient !-------------------------------------------------------------------------------------------------- -module subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, ip, el) +module subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, co, ip, el) integer, intent(in) :: & - ipc, & !< grain number + co, & !< grain number ip, & !< integration point number el !< element number real(pReal), intent(in), dimension(3,3) :: & @@ -138,7 +138,7 @@ module subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S traction_d, traction_t, traction_n, traction_crit, & udotd, dudotd_dt, udott, dudott_dt, udotn, dudotn_dt - phase = material_phaseAt(ipc,el) + phase = material_phaseAt(co,el) instance = kinematics_slipplane_opening_instance(phase) homog = material_homogenizationAt(el) damageOffset = material_homogenizationMemberAt(ip,el) diff --git a/src/kinematics_thermal_expansion.f90 b/src/kinematics_thermal_expansion.f90 index 5265d6172..6d4a39632 100644 --- a/src/kinematics_thermal_expansion.f90 +++ b/src/kinematics_thermal_expansion.f90 @@ -84,10 +84,10 @@ end function kinematics_thermal_expansion_init !-------------------------------------------------------------------------------------------------- !> @brief constitutive equation for calculating the velocity gradient !-------------------------------------------------------------------------------------------------- -module subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar, ipc, ip, el) +module subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar, co, ip, el) integer, intent(in) :: & - ipc, & !< grain number + co, & !< grain number ip, & !< integration point number el !< element number real(pReal), intent(out), dimension(3,3) :: & @@ -101,7 +101,7 @@ module subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar, i real(pReal) :: & T, TDot - phase = material_phaseAt(ipc,el) + phase = material_phaseAt(co,el) homog = material_homogenizationAt(el) T = temperature(homog)%p(material_homogenizationMemberAt(ip,el)) TDot = temperatureRate(homog)%p(material_homogenizationMemberAt(ip,el)) diff --git a/src/lattice.f90 b/src/lattice.f90 index 676232efe..6af135e4e 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -421,6 +421,8 @@ module lattice lattice_BCT_ID, & lattice_HEX_ID, & lattice_ORT_ID, & + lattice_equivalent_nu, & + lattice_equivalent_mu, & lattice_applyLatticeSymmetry33, & lattice_SchmidMatrix_slip, & lattice_SchmidMatrix_twin, & @@ -508,8 +510,8 @@ subroutine lattice_init lattice_C66(1:6,1:6,p) = applyLatticeSymmetryC66(lattice_C66(1:6,1:6,p),phase%get_asString('lattice')) - lattice_mu(p) = equivalent_mu(lattice_C66(1:6,1:6,p),'voigt') - lattice_nu(p) = equivalent_nu(lattice_C66(1:6,1:6,p),'voigt') + lattice_nu(p) = lattice_equivalent_nu(lattice_C66(1:6,1:6,p),'voigt') + lattice_mu(p) = lattice_equivalent_mu(lattice_C66(1:6,1:6,p),'voigt') lattice_C66(1:6,1:6,p) = math_sym3333to66(math_Voigt66to3333(lattice_C66(1:6,1:6,p))) ! Literature data is in Voigt notation do i = 1, 6 @@ -2188,15 +2190,16 @@ end function getlabels !> @brief Equivalent Poisson's ratio (ν) !> @details https://doi.org/10.1143/JPSJ.20.635 !-------------------------------------------------------------------------------------------------- -function equivalent_nu(C,assumption) result(nu) +function lattice_equivalent_nu(C,assumption) result(nu) real(pReal), dimension(6,6), intent(in) :: C !< Stiffness tensor (Voigt notation) character(len=*), intent(in) :: assumption !< Assumption ('Voigt' = isostrain, 'Reuss' = isostress) - real(pReal) :: K, mu, nu + logical :: error real(pReal), dimension(6,6) :: S + if (IO_lc(assumption) == 'voigt') then K = (C(1,1)+C(2,2)+C(3,3) +2.0_pReal*(C(1,2)+C(2,3)+C(1,3))) & / 9.0_pReal @@ -2210,25 +2213,26 @@ function equivalent_nu(C,assumption) result(nu) K = 0.0_pReal endif - mu = equivalent_mu(C,assumption) + mu = lattice_equivalent_mu(C,assumption) nu = (1.5_pReal*K -mu)/(3.0_pReal*K+mu) -end function equivalent_nu +end function lattice_equivalent_nu !-------------------------------------------------------------------------------------------------- !> @brief Equivalent shear modulus (μ) !> @details https://doi.org/10.1143/JPSJ.20.635 !-------------------------------------------------------------------------------------------------- -function equivalent_mu(C,assumption) result(mu) +function lattice_equivalent_mu(C,assumption) result(mu) real(pReal), dimension(6,6), intent(in) :: C !< Stiffness tensor (Voigt notation) character(len=*), intent(in) :: assumption !< Assumption ('Voigt' = isostrain, 'Reuss' = isostress) - real(pReal) :: mu + logical :: error real(pReal), dimension(6,6) :: S + if (IO_lc(assumption) == 'voigt') then mu = (1.0_pReal*(C(1,1)+C(2,2)+C(3,3)) -1.0_pReal*(C(1,2)+C(2,3)+C(1,3)) +3.0_pReal*(C(4,4)+C(5,5)+C(6,6))) & / 15.0_pReal @@ -2242,7 +2246,7 @@ function equivalent_mu(C,assumption) result(mu) mu = 0.0_pReal endif -end function equivalent_mu +end function lattice_equivalent_mu !-------------------------------------------------------------------------------------------------- @@ -2266,14 +2270,14 @@ subroutine selfTest call random_number(C) C(1,1) = C(1,1) + 1.0_pReal C = applyLatticeSymmetryC66(C,'aP') - if(dNeq(C(6,6),equivalent_mu(C,'voigt'),1.0e-12_pReal)) error stop 'equivalent_mu/voigt' - if(dNeq(C(6,6),equivalent_mu(C,'voigt'),1.0e-12_pReal)) error stop 'equivalent_mu/reuss' + if(dNeq(C(6,6),lattice_equivalent_mu(C,'voigt'),1.0e-12_pReal)) error stop 'equivalent_mu/voigt' + if(dNeq(C(6,6),lattice_equivalent_mu(C,'voigt'),1.0e-12_pReal)) error stop 'equivalent_mu/reuss' lambda = C(1,2) - if(dNeq(lambda*0.5_pReal/(lambda+equivalent_mu(C,'voigt')),equivalent_nu(C,'voigt'),1.0e-12_pReal)) & - error stop 'equivalent_nu/voigt' - if(dNeq(lambda*0.5_pReal/(lambda+equivalent_mu(C,'reuss')),equivalent_nu(C,'reuss'),1.0e-12_pReal)) & - error stop 'equivalent_nu/reuss' + if(dNeq(lambda*0.5_pReal/(lambda+lattice_equivalent_mu(C,'voigt')), & + lattice_equivalent_nu(C,'voigt'),1.0e-12_pReal)) error stop 'equivalent_nu/voigt' + if(dNeq(lambda*0.5_pReal/(lambda+lattice_equivalent_mu(C,'reuss')), & + lattice_equivalent_nu(C,'reuss'),1.0e-12_pReal)) error stop 'equivalent_nu/reuss' end subroutine selfTest diff --git a/src/marc/discretization_marc.f90 b/src/marc/discretization_marc.f90 index ca0b54b73..675e57bd3 100644 --- a/src/marc/discretization_marc.f90 +++ b/src/marc/discretization_marc.f90 @@ -12,7 +12,6 @@ module discretization_marc use DAMASK_interface use IO use config - use FEsolving use element use discretization use geometry_plastic_nonlocal @@ -89,9 +88,6 @@ subroutine discretization_marc_init if (debug_e < 1 .or. debug_e > nElems) call IO_error(602,ext_msg='element') if (debug_i < 1 .or. debug_i > elem%nIPs) call IO_error(602,ext_msg='IP') - FEsolving_execElem = [1,nElems] - FEsolving_execIP = [1,elem%nIPs] - allocate(cellNodeDefinition(elem%nNodes-1)) allocate(connectivity_cell(elem%NcellNodesPerCell,elem%nIPs,nElems)) call buildCells(connectivity_cell,cellNodeDefinition,& diff --git a/src/material.f90 b/src/material.f90 index 1f2437ad3..581182d22 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -17,29 +17,6 @@ module material private enum, bind(c); enumerator :: & - ELASTICITY_UNDEFINED_ID, & - ELASTICITY_HOOKE_ID, & - PLASTICITY_UNDEFINED_ID, & - PLASTICITY_NONE_ID, & - PLASTICITY_ISOTROPIC_ID, & - PLASTICITY_PHENOPOWERLAW_ID, & - PLASTICITY_KINEHARDENING_ID, & - PLASTICITY_DISLOTWIN_ID, & - PLASTICITY_DISLOTUNGSTEN_ID, & - PLASTICITY_NONLOCAL_ID, & - SOURCE_UNDEFINED_ID ,& - SOURCE_THERMAL_DISSIPATION_ID, & - SOURCE_THERMAL_EXTERNALHEAT_ID, & - SOURCE_DAMAGE_ISOBRITTLE_ID, & - SOURCE_DAMAGE_ISODUCTILE_ID, & - SOURCE_DAMAGE_ANISOBRITTLE_ID, & - SOURCE_DAMAGE_ANISODUCTILE_ID, & - KINEMATICS_UNDEFINED_ID ,& - KINEMATICS_CLEAVAGE_OPENING_ID, & - KINEMATICS_SLIPPLANE_OPENING_ID, & - KINEMATICS_THERMAL_EXPANSION_ID, & - STIFFNESS_DEGRADATION_UNDEFINED_ID, & - STIFFNESS_DEGRADATION_DAMAGE_ID, & THERMAL_ISOTHERMAL_ID, & THERMAL_CONDUCTION_ID, & DAMAGE_NONE_ID, & @@ -96,29 +73,6 @@ module material public :: & material_init, & - ELASTICITY_UNDEFINED_ID, & - ELASTICITY_HOOKE_ID, & - PLASTICITY_UNDEFINED_ID, & - PLASTICITY_NONE_ID, & - PLASTICITY_ISOTROPIC_ID, & - PLASTICITY_PHENOPOWERLAW_ID, & - PLASTICITY_KINEHARDENING_ID, & - PLASTICITY_DISLOTWIN_ID, & - PLASTICITY_DISLOTUNGSTEN_ID, & - PLASTICITY_NONLOCAL_ID, & - SOURCE_UNDEFINED_ID ,& - SOURCE_THERMAL_DISSIPATION_ID, & - SOURCE_THERMAL_EXTERNALHEAT_ID, & - SOURCE_DAMAGE_ISOBRITTLE_ID, & - SOURCE_DAMAGE_ISODUCTILE_ID, & - SOURCE_DAMAGE_ANISOBRITTLE_ID, & - SOURCE_DAMAGE_ANISODUCTILE_ID, & - KINEMATICS_UNDEFINED_ID ,& - KINEMATICS_CLEAVAGE_OPENING_ID, & - KINEMATICS_SLIPPLANE_OPENING_ID, & - KINEMATICS_THERMAL_EXPANSION_ID, & - STIFFNESS_DEGRADATION_UNDEFINED_ID, & - STIFFNESS_DEGRADATION_DAMAGE_ID, & THERMAL_ISOTHERMAL_ID, & THERMAL_CONDUCTION_ID, & DAMAGE_NONE_ID, & diff --git a/src/math.f90 b/src/math.f90 index 8005b5406..6b89a9923 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -279,9 +279,12 @@ real(pReal) pure function math_LeviCivita(i,j,k) integer, intent(in) :: i,j,k - if (all([i,j,k] == [1,2,3]) .or. all([i,j,k] == [2,3,1]) .or. all([i,j,k] == [3,1,2])) then + integer :: o + + + if (any([(all(cshift([i,j,k],o) == [1,2,3]),o=0,2)])) then math_LeviCivita = +1.0_pReal - elseif (all([i,j,k] == [3,2,1]) .or. all([i,j,k] == [2,1,3]) .or. all([i,j,k] == [1,3,2])) then + elseif (any([(all(cshift([i,j,k],o) == [3,2,1]),o=0,2)])) then math_LeviCivita = -1.0_pReal else math_LeviCivita = 0.0_pReal diff --git a/src/mesh/DAMASK_mesh.f90 b/src/mesh/DAMASK_mesh.f90 index 1e353892c..7369520c1 100644 --- a/src/mesh/DAMASK_mesh.f90 +++ b/src/mesh/DAMASK_mesh.f90 @@ -15,7 +15,6 @@ program DAMASK_mesh use IO use math use CPFEM2 - use FEsolving use config use discretization_mesh use FEM_Utilities diff --git a/src/mesh/FEM_utilities.f90 b/src/mesh/FEM_utilities.f90 index cb81f1f0c..2f3633e11 100644 --- a/src/mesh/FEM_utilities.f90 +++ b/src/mesh/FEM_utilities.f90 @@ -160,7 +160,7 @@ subroutine utilities_constitutiveResponse(timeinc,P_av,forwardData) print'(/,a)', ' ... evaluating constitutive response ......................................' - call materialpoint_stressAndItsTangent(timeinc) ! calculate P field + call materialpoint_stressAndItsTangent(timeinc,[1,mesh_maxNips],[1,mesh_NcpElems]) ! calculate P field cutBack = .false. ! reset cutBack status diff --git a/src/mesh/discretization_mesh.f90 b/src/mesh/discretization_mesh.f90 index 7dbd05e46..21c5feace 100644 --- a/src/mesh/discretization_mesh.f90 +++ b/src/mesh/discretization_mesh.f90 @@ -18,7 +18,6 @@ module discretization_mesh use config use discretization use results - use FEsolving use FEM_quadrature use YAML_types use prec @@ -30,7 +29,7 @@ module discretization_mesh mesh_Nboundaries, & mesh_NcpElemsGlobal - integer :: & + integer, public, protected :: & mesh_NcpElems !< total number of CP elements in mesh !!!! BEGIN DEPRECATED !!!!! @@ -174,9 +173,6 @@ subroutine discretization_mesh_init(restart) if (debug_element < 1 .or. debug_element > mesh_NcpElems) call IO_error(602,ext_msg='element') if (debug_ip < 1 .or. debug_ip > mesh_maxNips) call IO_error(602,ext_msg='IP') - FEsolving_execElem = [1,mesh_NcpElems] ! parallel loop bounds set to comprise all DAMASK elements - FEsolving_execIP = [1,mesh_maxNips] - allocate(mesh_node0(3,mesh_Nnodes),source=0.0_pReal) call discretization_init(materialAt,& diff --git a/src/prec.f90 b/src/prec.f90 index 95b1116cd..4d73462c4 100644 --- a/src/prec.f90 +++ b/src/prec.f90 @@ -108,8 +108,10 @@ logical elemental pure function dEq(a,b,tol) real(pReal), intent(in) :: a,b real(pReal), intent(in), optional :: tol + real(pReal) :: eps + if (present(tol)) then eps = tol else @@ -132,11 +134,8 @@ logical elemental pure function dNeq(a,b,tol) real(pReal), intent(in) :: a,b real(pReal), intent(in), optional :: tol - if (present(tol)) then - dNeq = .not. dEq(a,b,tol) - else - dNeq = .not. dEq(a,b) - endif + + dNeq = .not. dEq(a,b,tol) end function dNeq @@ -151,8 +150,10 @@ logical elemental pure function dEq0(a,tol) real(pReal), intent(in) :: a real(pReal), intent(in), optional :: tol + real(pReal) :: eps + if (present(tol)) then eps = tol else @@ -175,11 +176,8 @@ logical elemental pure function dNeq0(a,tol) real(pReal), intent(in) :: a real(pReal), intent(in), optional :: tol - if (present(tol)) then - dNeq0 = .not. dEq0(a,tol) - else - dNeq0 = .not. dEq0(a) - endif + + dNeq0 = .not. dEq0(a,tol) end function dNeq0 @@ -195,8 +193,10 @@ logical elemental pure function cEq(a,b,tol) complex(pReal), intent(in) :: a,b real(pReal), intent(in), optional :: tol + real(pReal) :: eps + if (present(tol)) then eps = tol else @@ -220,11 +220,8 @@ logical elemental pure function cNeq(a,b,tol) complex(pReal), intent(in) :: a,b real(pReal), intent(in), optional :: tol - if (present(tol)) then - cNeq = .not. cEq(a,b,tol) - else - cNeq = .not. cEq(a,b) - endif + + cNeq = .not. cEq(a,b,tol) end function cNeq @@ -238,6 +235,7 @@ pure function prec_bytesToC_FLOAT(bytes) real(C_FLOAT), dimension(size(bytes,kind=pI64)/(storage_size(0._C_FLOAT,pI64)/8_pI64)) :: & prec_bytesToC_FLOAT + prec_bytesToC_FLOAT = transfer(bytes,prec_bytesToC_FLOAT,size(prec_bytesToC_FLOAT)) end function prec_bytesToC_FLOAT @@ -252,6 +250,7 @@ pure function prec_bytesToC_DOUBLE(bytes) real(C_DOUBLE), dimension(size(bytes,kind=pI64)/(storage_size(0._C_DOUBLE,pI64)/8_pI64)) :: & prec_bytesToC_DOUBLE + prec_bytesToC_DOUBLE = transfer(bytes,prec_bytesToC_DOUBLE,size(prec_bytesToC_DOUBLE)) end function prec_bytesToC_DOUBLE @@ -266,6 +265,7 @@ pure function prec_bytesToC_INT32_T(bytes) integer(C_INT32_T), dimension(size(bytes,kind=pI64)/(storage_size(0_C_INT32_T,pI64)/8_pI64)) :: & prec_bytesToC_INT32_T + prec_bytesToC_INT32_T = transfer(bytes,prec_bytesToC_INT32_T,size(prec_bytesToC_INT32_T)) end function prec_bytesToC_INT32_T @@ -280,6 +280,7 @@ pure function prec_bytesToC_INT64_T(bytes) integer(C_INT64_T), dimension(size(bytes,kind=pI64)/(storage_size(0_C_INT64_T,pI64)/8_pI64)) :: & prec_bytesToC_INT64_T + prec_bytesToC_INT64_T = transfer(bytes,prec_bytesToC_INT64_T,size(prec_bytesToC_INT64_T)) end function prec_bytesToC_INT64_T @@ -295,6 +296,7 @@ subroutine selfTest integer(pInt), dimension(1) :: i real(pReal), dimension(2) :: r + realloc_lhs_test = [1,2] if (any(realloc_lhs_test/=[1,2])) error stop 'LHS allocation' diff --git a/src/quaternions.f90 b/src/quaternions.f90 deleted file mode 100644 index c5c43e3c1..000000000 --- a/src/quaternions.f90 +++ /dev/null @@ -1,534 +0,0 @@ -!--------------------------------------------------------------------------------------------------- -!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH -!> @author Philip Eisenlohr, Michigan State University -!> @brief general quaternion math, not limited to unit quaternions -!> @details w is the real part, (x, y, z) are the imaginary parts. -!> @details https://en.wikipedia.org/wiki/Quaternion -!--------------------------------------------------------------------------------------------------- -module quaternions - use prec - - implicit none - private - - real(pReal), parameter, public :: P = -1.0_pReal !< parameter for orientation conversion. - - type, public :: quaternion - real(pReal), private :: w = 0.0_pReal - real(pReal), private :: x = 0.0_pReal - real(pReal), private :: y = 0.0_pReal - real(pReal), private :: z = 0.0_pReal - - - contains - procedure, private :: add__ - procedure, private :: pos__ - generic, public :: operator(+) => add__,pos__ - - procedure, private :: sub__ - procedure, private :: neg__ - generic, public :: operator(-) => sub__,neg__ - - procedure, private :: mul_quat__ - procedure, private :: mul_scal__ - generic, public :: operator(*) => mul_quat__, mul_scal__ - - procedure, private :: div_quat__ - procedure, private :: div_scal__ - generic, public :: operator(/) => div_quat__, div_scal__ - - procedure, private :: eq__ - generic, public :: operator(==) => eq__ - - procedure, private :: neq__ - generic, public :: operator(/=) => neq__ - - procedure, private :: pow_quat__ - procedure, private :: pow_scal__ - generic, public :: operator(**) => pow_quat__, pow_scal__ - - procedure, public :: abs => abs__ - procedure, public :: conjg => conjg__ - procedure, public :: real => real__ - procedure, public :: aimag => aimag__ - - procedure, public :: homomorphed - procedure, public :: asArray - procedure, public :: inverse - - end type - - interface assignment (=) - module procedure assign_quat__ - module procedure assign_vec__ - end interface assignment (=) - - interface quaternion - module procedure init__ - end interface quaternion - - interface abs - procedure abs__ - end interface abs - - interface dot_product - procedure dot_product__ - end interface dot_product - - interface conjg - module procedure conjg__ - end interface conjg - - interface exp - module procedure exp__ - end interface exp - - interface log - module procedure log__ - end interface log - - interface real - module procedure real__ - end interface real - - interface aimag - module procedure aimag__ - end interface aimag - - public :: & - quaternions_init, & - assignment(=), & - conjg, aimag, & - log, exp, & - abs, dot_product, & - inverse, & - real - -contains - - -!-------------------------------------------------------------------------------------------------- -!> @brief Do self test. -!-------------------------------------------------------------------------------------------------- -subroutine quaternions_init - - print'(/,a)', ' <<<+- quaternions init -+>>>'; flush(6) - - call selfTest - -end subroutine quaternions_init - - -!--------------------------------------------------------------------------------------------------- -!> @brief construct a quaternion from a 4-vector -!--------------------------------------------------------------------------------------------------- -type(quaternion) pure function init__(array) - - real(pReal), intent(in), dimension(4) :: array - - init__%w = array(1) - init__%x = array(2) - init__%y = array(3) - init__%z = array(4) - -end function init__ - - -!--------------------------------------------------------------------------------------------------- -!> @brief assign a quaternion -!--------------------------------------------------------------------------------------------------- -elemental pure subroutine assign_quat__(self,other) - - type(quaternion), intent(out) :: self - type(quaternion), intent(in) :: other - - self = [other%w,other%x,other%y,other%z] - -end subroutine assign_quat__ - - -!--------------------------------------------------------------------------------------------------- -!> @brief assign a 4-vector -!--------------------------------------------------------------------------------------------------- -pure subroutine assign_vec__(self,other) - - type(quaternion), intent(out) :: self - real(pReal), intent(in), dimension(4) :: other - - self%w = other(1) - self%x = other(2) - self%y = other(3) - self%z = other(4) - -end subroutine assign_vec__ - - -!--------------------------------------------------------------------------------------------------- -!> @brief add a quaternion -!--------------------------------------------------------------------------------------------------- -type(quaternion) elemental pure function add__(self,other) - - class(quaternion), intent(in) :: self,other - - add__ = [ self%w, self%x, self%y ,self%z] & - + [other%w, other%x, other%y,other%z] - -end function add__ - - -!--------------------------------------------------------------------------------------------------- -!> @brief return (unary positive operator) -!--------------------------------------------------------------------------------------------------- -type(quaternion) elemental pure function pos__(self) - - class(quaternion), intent(in) :: self - - pos__ = self * (+1.0_pReal) - -end function pos__ - - -!--------------------------------------------------------------------------------------------------- -!> @brief subtract a quaternion -!--------------------------------------------------------------------------------------------------- -type(quaternion) elemental pure function sub__(self,other) - - class(quaternion), intent(in) :: self,other - - sub__ = [ self%w, self%x, self%y ,self%z] & - - [other%w, other%x, other%y,other%z] - -end function sub__ - - -!--------------------------------------------------------------------------------------------------- -!> @brief negate (unary negative operator) -!--------------------------------------------------------------------------------------------------- -type(quaternion) elemental pure function neg__(self) - - class(quaternion), intent(in) :: self - - neg__ = self * (-1.0_pReal) - -end function neg__ - - -!--------------------------------------------------------------------------------------------------- -!> @brief multiply with a quaternion -!--------------------------------------------------------------------------------------------------- -type(quaternion) elemental pure function mul_quat__(self,other) - - class(quaternion), intent(in) :: self, other - - mul_quat__%w = self%w*other%w - self%x*other%x - self%y*other%y - self%z*other%z - mul_quat__%x = self%w*other%x + self%x*other%w + P * (self%y*other%z - self%z*other%y) - mul_quat__%y = self%w*other%y + self%y*other%w + P * (self%z*other%x - self%x*other%z) - mul_quat__%z = self%w*other%z + self%z*other%w + P * (self%x*other%y - self%y*other%x) - -end function mul_quat__ - - -!--------------------------------------------------------------------------------------------------- -!> @brief multiply with a scalar -!--------------------------------------------------------------------------------------------------- -type(quaternion) elemental pure function mul_scal__(self,scal) - - class(quaternion), intent(in) :: self - real(pReal), intent(in) :: scal - - mul_scal__ = [self%w,self%x,self%y,self%z]*scal - -end function mul_scal__ - - -!--------------------------------------------------------------------------------------------------- -!> @brief divide by a quaternion -!--------------------------------------------------------------------------------------------------- -type(quaternion) elemental pure function div_quat__(self,other) - - class(quaternion), intent(in) :: self, other - - div_quat__ = self * (conjg(other)/(abs(other)**2.0_pReal)) - -end function div_quat__ - - -!--------------------------------------------------------------------------------------------------- -!> @brief divide by a scalar -!--------------------------------------------------------------------------------------------------- -type(quaternion) elemental pure function div_scal__(self,scal) - - class(quaternion), intent(in) :: self - real(pReal), intent(in) :: scal - - div_scal__ = [self%w,self%x,self%y,self%z]/scal - -end function div_scal__ - - -!--------------------------------------------------------------------------------------------------- -!> @brief test equality -!--------------------------------------------------------------------------------------------------- -logical elemental pure function eq__(self,other) - - class(quaternion), intent(in) :: self,other - - eq__ = all(dEq([ self%w, self%x, self%y, self%z], & - [other%w,other%x,other%y,other%z])) - -end function eq__ - - -!--------------------------------------------------------------------------------------------------- -!> @brief test inequality -!--------------------------------------------------------------------------------------------------- -logical elemental pure function neq__(self,other) - - class(quaternion), intent(in) :: self,other - - neq__ = .not. self%eq__(other) - -end function neq__ - - -!--------------------------------------------------------------------------------------------------- -!> @brief raise to the power of a quaternion -!--------------------------------------------------------------------------------------------------- -type(quaternion) elemental pure function pow_quat__(self,expon) - - class(quaternion), intent(in) :: self - type(quaternion), intent(in) :: expon - - pow_quat__ = exp(log(self)*expon) - -end function pow_quat__ - - -!--------------------------------------------------------------------------------------------------- -!> @brief raise to the power of a scalar -!--------------------------------------------------------------------------------------------------- -type(quaternion) elemental pure function pow_scal__(self,expon) - - class(quaternion), intent(in) :: self - real(pReal), intent(in) :: expon - - pow_scal__ = exp(log(self)*expon) - -end function pow_scal__ - - -!--------------------------------------------------------------------------------------------------- -!> @brief take exponential -!--------------------------------------------------------------------------------------------------- -type(quaternion) elemental pure function exp__(a) - - class(quaternion), intent(in) :: a - real(pReal) :: absImag - - absImag = norm2(aimag(a)) - - exp__ = merge(exp(a%w) * [ cos(absImag), & - a%x/absImag * sin(absImag), & - a%y/absImag * sin(absImag), & - a%z/absImag * sin(absImag)], & - IEEE_value(1.0_pReal,IEEE_SIGNALING_NAN), & - dNeq0(absImag)) - -end function exp__ - - -!--------------------------------------------------------------------------------------------------- -!> @brief take logarithm -!--------------------------------------------------------------------------------------------------- -type(quaternion) elemental pure function log__(a) - - class(quaternion), intent(in) :: a - real(pReal) :: absImag - - absImag = norm2(aimag(a)) - - log__ = merge([log(abs(a)), & - a%x/absImag * acos(a%w/abs(a)), & - a%y/absImag * acos(a%w/abs(a)), & - a%z/absImag * acos(a%w/abs(a))], & - IEEE_value(1.0_pReal,IEEE_SIGNALING_NAN), & - dNeq0(absImag)) - -end function log__ - - -!--------------------------------------------------------------------------------------------------- -!> @brief return norm -!--------------------------------------------------------------------------------------------------- -real(pReal) elemental pure function abs__(self) - - class(quaternion), intent(in) :: self - - abs__ = norm2([self%w,self%x,self%y,self%z]) - -end function abs__ - - -!--------------------------------------------------------------------------------------------------- -!> @brief calculate dot product -!--------------------------------------------------------------------------------------------------- -real(pReal) elemental pure function dot_product__(a,b) - - class(quaternion), intent(in) :: a,b - - dot_product__ = a%w*b%w + a%x*b%x + a%y*b%y + a%z*b%z - -end function dot_product__ - - -!--------------------------------------------------------------------------------------------------- -!> @brief take conjugate complex -!--------------------------------------------------------------------------------------------------- -type(quaternion) elemental pure function conjg__(self) - - class(quaternion), intent(in) :: self - - conjg__ = [self%w,-self%x,-self%y,-self%z] - -end function conjg__ - - -!--------------------------------------------------------------------------------------------------- -!> @brief homomorph -!--------------------------------------------------------------------------------------------------- -type(quaternion) elemental pure function homomorphed(self) - - class(quaternion), intent(in) :: self - - homomorphed = - self - -end function homomorphed - - -!--------------------------------------------------------------------------------------------------- -!> @brief return as plain array -!--------------------------------------------------------------------------------------------------- -pure function asArray(self) - - real(pReal), dimension(4) :: asArray - class(quaternion), intent(in) :: self - - asArray = [self%w,self%x,self%y,self%z] - -end function asArray - - -!--------------------------------------------------------------------------------------------------- -!> @brief real part (scalar) -!--------------------------------------------------------------------------------------------------- -pure function real__(self) - - real(pReal) :: real__ - class(quaternion), intent(in) :: self - - real__ = self%w - -end function real__ - - -!--------------------------------------------------------------------------------------------------- -!> @brief imaginary part (3-vector) -!--------------------------------------------------------------------------------------------------- -pure function aimag__(self) - - real(pReal), dimension(3) :: aimag__ - class(quaternion), intent(in) :: self - - aimag__ = [self%x,self%y,self%z] - -end function aimag__ - - -!--------------------------------------------------------------------------------------------------- -!> @brief inverse -!--------------------------------------------------------------------------------------------------- -type(quaternion) elemental pure function inverse(self) - - class(quaternion), intent(in) :: self - - inverse = conjg(self)/abs(self)**2.0_pReal - -end function inverse - - -!-------------------------------------------------------------------------------------------------- -!> @brief check correctness of some quaternions functions -!-------------------------------------------------------------------------------------------------- -subroutine selfTest - - real(pReal), dimension(4) :: qu - type(quaternion) :: q, q_2 - - if(dNeq(abs(P),1.0_pReal)) error stop 'P not in {-1,+1}' - - call random_number(qu) - qu = (qu-0.5_pReal) * 2.0_pReal - q = quaternion(qu) - - q_2= qu - if(any(dNeq(q%asArray(),q_2%asArray()))) error stop 'assign_vec__' - - q_2 = q + q - if(any(dNeq(q_2%asArray(),2.0_pReal*qu))) error stop 'add__' - - q_2 = q - q - if(any(dNeq0(q_2%asArray()))) error stop 'sub__' - - q_2 = q * 5.0_pReal - if(any(dNeq(q_2%asArray(),5.0_pReal*qu))) error stop 'mul__' - - q_2 = q / 0.5_pReal - if(any(dNeq(q_2%asArray(),2.0_pReal*qu))) error stop 'div__' - - q_2 = q * 0.3_pReal - if(dNeq0(abs(q)) .and. q_2 == q) error stop 'eq__' - - q_2 = q - if(q_2 /= q) error stop 'neq__' - - if(dNeq(abs(q),norm2(qu))) error stop 'abs__' - if(dNeq(abs(q)**2.0_pReal, real(q*q%conjg()),1.0e-14_pReal)) & - error stop 'abs__/*conjg' - - if(any(dNeq(q%asArray(),qu))) error stop 'eq__' - if(dNeq(q%real(), qu(1))) error stop 'real()' - if(any(dNeq(q%aimag(), qu(2:4)))) error stop 'aimag()' - - q_2 = q%homomorphed() - if(q /= q_2* (-1.0_pReal)) error stop 'homomorphed' - if(dNeq(q_2%real(), qu(1)* (-1.0_pReal))) error stop 'homomorphed/real' - if(any(dNeq(q_2%aimag(),qu(2:4)*(-1.0_pReal)))) error stop 'homomorphed/aimag' - - q_2 = conjg(q) - if(dNeq(abs(q),abs(q_2))) error stop 'conjg/abs' - if(q /= conjg(q_2)) error stop 'conjg/involution' - if(dNeq(q_2%real(), q%real())) error stop 'conjg/real' - if(any(dNeq(q_2%aimag(),q%aimag()*(-1.0_pReal)))) error stop 'conjg/aimag' - - if(abs(q) > 0.0_pReal) then - q_2 = q * q%inverse() - if( dNeq(real(q_2), 1.0_pReal,1.0e-15_pReal)) error stop 'inverse/real' - if(any(dNeq0(aimag(q_2), 1.0e-15_pReal))) error stop 'inverse/aimag' - - q_2 = q/abs(q) - q_2 = conjg(q_2) - inverse(q_2) - if(any(dNeq0(q_2%asArray(),1.0e-15_pReal))) error stop 'inverse/conjg' - endif - if(dNeq(dot_product(qu,qu),dot_product(q,q))) error stop 'dot_product' - -#if !(defined(__GFORTRAN__) && __GNUC__ < 9) - if (norm2(aimag(q)) > 0.0_pReal) then - if (dNeq0(abs(q-exp(log(q))),1.0e-13_pReal)) error stop 'exp/log' - if (dNeq0(abs(q-log(exp(q))),1.0e-13_pReal)) error stop 'log/exp' - endif -#endif - -end subroutine selfTest - - -end module quaternions diff --git a/src/results.f90 b/src/results.f90 index ea9fd62d4..6363e3efc 100644 --- a/src/results.f90 +++ b/src/results.f90 @@ -111,8 +111,6 @@ subroutine results_addIncrement(inc,time) call results_closeGroup(results_addGroup(trim('inc'//trim(adjustl(incChar))))) call results_setLink(trim('inc'//trim(adjustl(incChar))),'current') call results_addAttribute('time/s',time,trim('inc'//trim(adjustl(incChar)))) - call results_closeGroup(results_addGroup('current/phase')) - call results_closeGroup(results_addGroup('current/homogenization')) end subroutine results_addIncrement diff --git a/src/rotations.f90 b/src/rotations.f90 index 888e73762..57dd16b53 100644 --- a/src/rotations.f90 +++ b/src/rotations.f90 @@ -47,16 +47,16 @@ !--------------------------------------------------------------------------------------------------- module rotations - use prec use IO use math - use quaternions implicit none private + real(pReal), parameter :: P = -1.0_pReal !< parameter for orientation conversion. + type, public :: rotation - type(quaternion) :: q + real(pReal), dimension(4) :: q contains procedure, public :: asQuaternion procedure, public :: asEulers @@ -103,7 +103,6 @@ contains !-------------------------------------------------------------------------------------------------- subroutine rotations_init - call quaternions_init print'(/,a)', ' <<<+- rotations init -+>>>'; flush(IO_STDOUT) print*, 'Rowenhorst et al., Modelling and Simulation in Materials Science and Engineering 23:083501, 2015' @@ -122,7 +121,7 @@ pure function asQuaternion(self) class(rotation), intent(in) :: self real(pReal), dimension(4) :: asQuaternion - asQuaternion = self%q%asArray() + asQuaternion = self%q end function asQuaternion !--------------------------------------------------------------------------------------------------- @@ -131,7 +130,7 @@ pure function asEulers(self) class(rotation), intent(in) :: self real(pReal), dimension(3) :: asEulers - asEulers = qu2eu(self%q%asArray()) + asEulers = qu2eu(self%q) end function asEulers !--------------------------------------------------------------------------------------------------- @@ -140,7 +139,7 @@ pure function asAxisAngle(self) class(rotation), intent(in) :: self real(pReal), dimension(4) :: asAxisAngle - asAxisAngle = qu2ax(self%q%asArray()) + asAxisAngle = qu2ax(self%q) end function asAxisAngle !--------------------------------------------------------------------------------------------------- @@ -149,7 +148,7 @@ pure function asMatrix(self) class(rotation), intent(in) :: self real(pReal), dimension(3,3) :: asMatrix - asMatrix = qu2om(self%q%asArray()) + asMatrix = qu2om(self%q) end function asMatrix !--------------------------------------------------------------------------------------------------- @@ -158,7 +157,7 @@ pure function asRodrigues(self) class(rotation), intent(in) :: self real(pReal), dimension(4) :: asRodrigues - asRodrigues = qu2ro(self%q%asArray()) + asRodrigues = qu2ro(self%q) end function asRodrigues !--------------------------------------------------------------------------------------------------- @@ -167,7 +166,7 @@ pure function asHomochoric(self) class(rotation), intent(in) :: self real(pReal), dimension(3) :: asHomochoric - asHomochoric = qu2ho(self%q%asArray()) + asHomochoric = qu2ho(self%q) end function asHomochoric @@ -259,7 +258,7 @@ pure elemental function rotRot__(self,R) result(rRot) type(rotation) :: rRot class(rotation), intent(in) :: self,R - rRot = rotation(self%q*R%q) + rRot = rotation(multiply_quaternion(self%q,R%q)) call rRot%standardize() end function rotRot__ @@ -272,14 +271,14 @@ pure elemental subroutine standardize(self) class(rotation), intent(inout) :: self - if (real(self%q) < 0.0_pReal) self%q = self%q%homomorphed() + if (self%q(1) < 0.0_pReal) self%q = - self%q end subroutine standardize !--------------------------------------------------------------------------------------------------- !> @author Marc De Graef, Carnegie Mellon University -!> @brief rotate a vector passively (default) or actively +!> @brief Rotate a vector passively (default) or actively. !--------------------------------------------------------------------------------------------------- pure function rotVector(self,v,active) result(vRot) @@ -288,9 +287,8 @@ pure function rotVector(self,v,active) result(vRot) real(pReal), intent(in), dimension(3) :: v logical, intent(in), optional :: active - real(pReal), dimension(3) :: v_normed - type(quaternion) :: q - logical :: passive + real(pReal), dimension(4) :: v_normed, q + logical :: passive if (present(active)) then passive = .not. active @@ -301,13 +299,13 @@ pure function rotVector(self,v,active) result(vRot) if (dEq0(norm2(v))) then vRot = v else - v_normed = v/norm2(v) + v_normed = [0.0_pReal,v]/norm2(v) if (passive) then - q = self%q * (quaternion([0.0_pReal, v_normed(1), v_normed(2), v_normed(3)]) * conjg(self%q) ) + q = multiply_quaternion(self%q, multiply_quaternion(v_normed, conjugate_quaternion(self%q))) else - q = conjg(self%q) * (quaternion([0.0_pReal, v_normed(1), v_normed(2), v_normed(3)]) * self%q ) + q = multiply_quaternion(conjugate_quaternion(self%q), multiply_quaternion(v_normed, self%q)) endif - vRot = q%aimag()*norm2(v) + vRot = q(2:4)*norm2(v) endif end function rotVector @@ -315,8 +313,8 @@ end function rotVector !--------------------------------------------------------------------------------------------------- !> @author Marc De Graef, Carnegie Mellon University -!> @brief rotate a rank-2 tensor passively (default) or actively -!> @details: rotation is based on rotation matrix +!> @brief Rotate a rank-2 tensor passively (default) or actively. +!> @details: Rotation is based on rotation matrix !--------------------------------------------------------------------------------------------------- pure function rotTensor2(self,T,active) result(tRot) @@ -403,7 +401,7 @@ pure elemental function misorientation(self,other) type(rotation) :: misorientation class(rotation), intent(in) :: self, other - misorientation%q = other%q * conjg(self%q) + misorientation%q = multiply_quaternion(other%q, conjugate_quaternion(self%q)) end function misorientation @@ -1338,7 +1336,7 @@ end function cu2ho !-------------------------------------------------------------------------- !> @author Marc De Graef, Carnegie Mellon University !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH -!> @brief determine to which pyramid a point in a cubic grid belongs +!> @brief Determine to which pyramid a point in a cubic grid belongs. !-------------------------------------------------------------------------- pure function GetPyramidOrder(xyz) @@ -1362,7 +1360,39 @@ end function GetPyramidOrder !-------------------------------------------------------------------------------------------------- -!> @brief check correctness of some rotations functions +!> @brief Multiply two quaternions. +!-------------------------------------------------------------------------------------------------- +pure function multiply_quaternion(qu1,qu2) + + real(pReal), dimension(4), intent(in) :: qu1, qu2 + real(pReal), dimension(4) :: multiply_quaternion + + + multiply_quaternion(1) = qu1(1)*qu2(1) - qu1(2)*qu2(2) - qu1(3)*qu2(3) - qu1(4)*qu2(4) + multiply_quaternion(2) = qu1(1)*qu2(2) + qu1(2)*qu2(1) + P * (qu1(3)*qu2(4) - qu1(4)*qu2(3)) + multiply_quaternion(3) = qu1(1)*qu2(3) + qu1(3)*qu2(1) + P * (qu1(4)*qu2(2) - qu1(2)*qu2(4)) + multiply_quaternion(4) = qu1(1)*qu2(4) + qu1(4)*qu2(1) + P * (qu1(2)*qu2(3) - qu1(3)*qu2(2)) + +end function multiply_quaternion + + +!-------------------------------------------------------------------------------------------------- +!> @brief Calculate conjugate complex of a quaternion. +!-------------------------------------------------------------------------------------------------- +pure function conjugate_quaternion(qu) + + real(pReal), dimension(4), intent(in) :: qu + real(pReal), dimension(4) :: conjugate_quaternion + + + conjugate_quaternion = [qu(1), -qu(2), -qu(3), -qu(4)] + + +end function conjugate_quaternion + + +!-------------------------------------------------------------------------------------------------- +!> @brief Check correctness of some rotations functions. !-------------------------------------------------------------------------------------------------- subroutine selfTest @@ -1374,7 +1404,8 @@ subroutine selfTest real :: A,B integer :: i - do i=1,10 + + do i = 1, 10 #if defined(__GFORTRAN__) && __GNUC__<9 if(i<7) cycle diff --git a/src/source_damage_anisoBrittle.f90 b/src/source_damage_anisoBrittle.f90 index 55d5546fc..0f923ceba 100644 --- a/src/source_damage_anisoBrittle.f90 +++ b/src/source_damage_anisoBrittle.f90 @@ -120,10 +120,10 @@ end function source_damage_anisoBrittle_init !-------------------------------------------------------------------------------------------------- !> @brief calculates derived quantities from state !-------------------------------------------------------------------------------------------------- -module subroutine source_damage_anisoBrittle_dotState(S, ipc, ip, el) +module subroutine source_damage_anisoBrittle_dotState(S, co, ip, el) integer, intent(in) :: & - ipc, & !< component-ID of integration point + co, & !< component-ID of integration point ip, & !< integration point el !< element real(pReal), intent(in), dimension(3,3) :: & @@ -139,8 +139,8 @@ module subroutine source_damage_anisoBrittle_dotState(S, ipc, ip, el) real(pReal) :: & traction_d, traction_t, traction_n, traction_crit - phase = material_phaseAt(ipc,el) - constituent = material_phasememberAt(ipc,ip,el) + phase = material_phaseAt(co,el) + constituent = material_phasememberAt(co,ip,el) sourceOffset = source_damage_anisoBrittle_offset(phase) homog = material_homogenizationAt(el) damageOffset = material_homogenizationMemberAt(ip,el) diff --git a/src/source_damage_anisoDuctile.f90 b/src/source_damage_anisoDuctile.f90 index 912fe1387..6f71fc145 100644 --- a/src/source_damage_anisoDuctile.f90 +++ b/src/source_damage_anisoDuctile.f90 @@ -107,10 +107,10 @@ end function source_damage_anisoDuctile_init !-------------------------------------------------------------------------------------------------- !> @brief calculates derived quantities from state !-------------------------------------------------------------------------------------------------- -module subroutine source_damage_anisoDuctile_dotState(ipc, ip, el) +module subroutine source_damage_anisoDuctile_dotState(co, ip, el) integer, intent(in) :: & - ipc, & !< component-ID of integration point + co, & !< component-ID of integration point ip, & !< integration point el !< element @@ -121,8 +121,8 @@ module subroutine source_damage_anisoDuctile_dotState(ipc, ip, el) damageOffset, & homog - phase = material_phaseAt(ipc,el) - constituent = material_phasememberAt(ipc,ip,el) + phase = material_phaseAt(co,el) + constituent = material_phasememberAt(co,ip,el) sourceOffset = source_damage_anisoDuctile_offset(phase) homog = material_homogenizationAt(el) damageOffset = material_homogenizationMemberAt(ip,el) diff --git a/src/source_damage_isoBrittle.f90 b/src/source_damage_isoBrittle.f90 index 7fcf17ee0..8c768b08d 100644 --- a/src/source_damage_isoBrittle.f90 +++ b/src/source_damage_isoBrittle.f90 @@ -94,10 +94,10 @@ end function source_damage_isoBrittle_init !-------------------------------------------------------------------------------------------------- !> @brief calculates derived quantities from state !-------------------------------------------------------------------------------------------------- -module subroutine source_damage_isoBrittle_deltaState(C, Fe, ipc, ip, el) +module subroutine source_damage_isoBrittle_deltaState(C, Fe, co, ip, el) integer, intent(in) :: & - ipc, & !< component-ID of integration point + co, & !< component-ID of integration point ip, & !< integration point el !< element real(pReal), intent(in), dimension(3,3) :: & @@ -114,8 +114,8 @@ module subroutine source_damage_isoBrittle_deltaState(C, Fe, ipc, ip, el) real(pReal) :: & strainenergy - phase = material_phaseAt(ipc,el) !< phase ID at ipc,ip,el - constituent = material_phasememberAt(ipc,ip,el) !< state array offset for phase ID at ipc,ip,el + phase = material_phaseAt(co,el) !< phase ID at co,ip,el + constituent = material_phasememberAt(co,ip,el) !< state array offset for phase ID at co,ip,el sourceOffset = source_damage_isoBrittle_offset(phase) strain = 0.5_pReal*math_sym33to6(matmul(transpose(Fe),Fe)-math_I3) diff --git a/src/source_damage_isoDuctile.f90 b/src/source_damage_isoDuctile.f90 index b66e220d9..86222bbf9 100644 --- a/src/source_damage_isoDuctile.f90 +++ b/src/source_damage_isoDuctile.f90 @@ -98,10 +98,10 @@ end function source_damage_isoDuctile_init !-------------------------------------------------------------------------------------------------- !> @brief calculates derived quantities from state !-------------------------------------------------------------------------------------------------- -module subroutine source_damage_isoDuctile_dotState(ipc, ip, el) +module subroutine source_damage_isoDuctile_dotState(co, ip, el) integer, intent(in) :: & - ipc, & !< component-ID of integration point + co, & !< component-ID of integration point ip, & !< integration point el !< element @@ -112,8 +112,8 @@ module subroutine source_damage_isoDuctile_dotState(ipc, ip, el) damageOffset, & homog - phase = material_phaseAt(ipc,el) - constituent = material_phasememberAt(ipc,ip,el) + phase = material_phaseAt(co,el) + constituent = material_phasememberAt(co,ip,el) sourceOffset = source_damage_isoDuctile_offset(phase) homog = material_homogenizationAt(el) damageOffset = material_homogenizationMemberAt(ip,el)