From 6d31f77debdb14de9bc2e8d4c7768d9f59a9a32e Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 22 May 2021 17:21:07 +0200 Subject: [PATCH] natural data structure --- src/discretization.f90 | 30 +++++++++++------------ src/material.f90 | 21 +++++++++++++++- src/phase.f90 | 28 ++++++++++++++++++--- src/phase_mechanical.f90 | 11 ++------- src/phase_mechanical_plastic_nonlocal.f90 | 4 +-- 5 files changed, 63 insertions(+), 31 deletions(-) diff --git a/src/discretization.f90 b/src/discretization.f90 index 22abbd644..cb74056cc 100644 --- a/src/discretization.f90 +++ b/src/discretization.f90 @@ -9,20 +9,20 @@ module discretization implicit none private - + integer, public, protected :: & discretization_nIPs, & discretization_Nelems - - integer, public, protected, dimension(:), allocatable :: & - discretization_materialAt !ToDo: discretization_materialID - real(pReal), public, protected, dimension(:,:), allocatable :: & + integer, public, protected, dimension(:), allocatable :: & + discretization_materialAt !ToDo: discretization_ID_material + + real(pReal), public, protected, dimension(:,:), allocatable :: & discretization_IPcoords0, & discretization_IPcoords, & discretization_NodeCoords0, & discretization_NodeCoords - + integer :: & discretization_sharedNodesBegin @@ -33,7 +33,7 @@ module discretization discretization_setNodeCoords contains - + !-------------------------------------------------------------------------------------------------- !> @brief stores the relevant information in globally accesible variables !-------------------------------------------------------------------------------------------------- @@ -54,20 +54,20 @@ subroutine discretization_init(materialAt,& discretization_Nelems = size(materialAt,1) discretization_nIPs = size(IPcoords0,2)/discretization_Nelems - discretization_materialAt = materialAt + discretization_materialAt = materialAt discretization_IPcoords0 = IPcoords0 discretization_IPcoords = IPcoords0 discretization_NodeCoords0 = NodeCoords0 discretization_NodeCoords = NodeCoords0 - + if(present(sharedNodesBegin)) then discretization_sharedNodesBegin = sharedNodesBegin else discretization_sharedNodesBegin = size(discretization_NodeCoords0,2) endif - + end subroutine discretization_init @@ -77,13 +77,13 @@ end subroutine discretization_init subroutine discretization_results real(pReal), dimension(:,:), allocatable :: u - + call results_closeGroup(results_addGroup('current/geometry')) - + u = discretization_NodeCoords (1:3,:discretization_sharedNodesBegin) & - discretization_NodeCoords0(1:3,:discretization_sharedNodesBegin) call results_writeDataset('current/geometry',u,'u_n','displacements of the nodes','m') - + u = discretization_IPcoords & - discretization_IPcoords0 call results_writeDataset('current/geometry',u,'u_p','displacements of the materialpoints (cell centers)','m') @@ -97,7 +97,7 @@ end subroutine discretization_results subroutine discretization_setIPcoords(IPcoords) real(pReal), dimension(:,:), intent(in) :: IPcoords - + discretization_IPcoords = IPcoords end subroutine discretization_setIPcoords @@ -109,7 +109,7 @@ end subroutine discretization_setIPcoords subroutine discretization_setNodeCoords(NodeCoords) real(pReal), dimension(:,:), intent(in) :: NodeCoords - + discretization_NodeCoords = NodeCoords end subroutine discretization_setNodeCoords diff --git a/src/material.f90 b/src/material.f90 index e072d9a1c..c62442fcd 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -16,6 +16,12 @@ module material implicit none private + type :: tRotationContainer + type(Rotation), dimension(:), allocatable :: data + end type + + type(tRotationContainer), dimension(:), allocatable :: material_orientation0_2 + integer, dimension(:), allocatable, public, protected :: & homogenization_Nconstituents !< number of grains in each homogenization integer, public, protected :: & @@ -39,6 +45,8 @@ module material material_phaseMemberAt !< position of the element within its phase instance public :: & + tRotationContainer, & + material_orientation0_2, & material_init contains @@ -88,7 +96,7 @@ subroutine parse() real(pReal) :: & frac integer :: & - el, ip, co, & + el, ip, co, ma, & h, ce materials => config_material%get('material') @@ -153,6 +161,17 @@ subroutine parse() enddo + allocate(material_orientation0_2(materials%length)) + + do ma = 1, materials%length + material => materials%get(ma) + constituents => material%get('constituents') + allocate(material_orientation0_2(ma)%data(constituents%length)) + do co = 1, constituents%length + call material_orientation0_2(ma)%data(co)%fromQuaternion(constituent%get_as1dFloat('O',requiredSize=4)) + enddo + enddo + end subroutine parse diff --git a/src/phase.f90 b/src/phase.f90 index f0b5c4e22..9d4a423d0 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -19,8 +19,9 @@ module phase implicit none private - type(Rotation), dimension(:,:,:), allocatable :: & - material_orientation0 !< initial orientation of each grain,IP,element + type(tRotationContainer), dimension(:), allocatable :: & + phase_orientation0, & + phase_orientation type(rotation), dimension(:,:,:), allocatable :: & crystallite_orientation !< current orientation @@ -346,7 +347,7 @@ contains subroutine phase_init integer :: & - ph + ph, ce, co, ma class (tNode), pointer :: & debug_constitutive, & materials, & @@ -367,6 +368,25 @@ subroutine phase_init materials => config_material%get('material') phases => config_material%get('phase') + + allocate(phase_orientation0(phases%length)) + do ph = 1,phases%length + allocate(phase_orientation0(ph)%data(count(material_phaseID==ph))) + enddo + + do ce = 1, size(material_phaseID,2) + ma = discretization_materialAt((ce-1)/discretization_nIPs+1) + do co = 1,homogenization_Nconstituents(material_homogenizationID(ce)) + ph = material_phaseID(co,ce) + phase_orientation0(ph)%data(material_phaseEntry(co,ce)) = material_orientation0_2(ma)%data(co) + enddo + enddo + + allocate(phase_orientation(phases%length)) + do ph = 1,phases%length + phase_orientation(ph)%data = phase_orientation0(ph)%data + enddo + call mechanical_init(materials,phases) call damage_init call thermal_init(phases) @@ -602,7 +622,7 @@ function crystallite_push33ToRef(co,ce, tensor33) ph = material_phaseID(co,ce) en = material_phaseEntry(co,ce) - T = matmul(material_orientation0(co,ph,en)%asMatrix(),transpose(math_inv33(phase_F(co,ce)))) ! ToDo: initial orientation correct? + T = matmul(phase_orientation0(ph)%data(en)%asMatrix(),transpose(math_inv33(phase_F(co,ce)))) ! ToDo: initial orientation correct? crystallite_push33ToRef = matmul(transpose(T),matmul(tensor33,T)) diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index cbaf149ec..de23d39fa 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -68,7 +68,7 @@ submodule(phase) mechanical dS_dFe, & !< derivative of 2nd P-K stress with respect to elastic deformation gradient dS_dFi !< derivative of 2nd P-K stress with respect to intermediate deformation gradient end subroutine phase_hooke_SandItsTangents - + module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,ph,en) real(pReal), dimension(3,3), intent(out) :: & Li !< inleastic velocity gradient @@ -204,7 +204,6 @@ module subroutine mechanical_init(materials,phases) num_crystallite, & material, & constituents, & - constituent, & phase, & mech @@ -228,8 +227,6 @@ module subroutine mechanical_init(materials,phases) allocate(phase_mechanical_P(phases%length)) allocate(phase_mechanical_S0(phases%length)) - allocate(material_orientation0(homogenization_maxNconstituents,phases%length,maxVal(material_phaseEntry))) - do ph = 1, phases%length Nmembers = count(material_phaseID == ph) @@ -260,15 +257,11 @@ module subroutine mechanical_init(materials,phases) do el = 1, size(material_phaseMemberAt,3); do ip = 1, size(material_phaseMemberAt,2) do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) material => materials%get(discretization_materialAt(el)) - constituents => material%get('constituents') - constituent => constituents%get(co) ph = material_phaseID(co,(el-1)*discretization_nIPs + ip) en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip) - call material_orientation0(co,ph,en)%fromQuaternion(constituent%get_as1dFloat('O',requiredSize=4)) - - phase_mechanical_Fp0(ph)%data(1:3,1:3,en) = material_orientation0(co,ph,en)%asMatrix() ! Fp reflects initial orientation (see 10.1016/j.actamat.2006.01.005) + phase_mechanical_Fp0(ph)%data(1:3,1:3,en) = phase_orientation0(ph)%data(en)%asMatrix() ! Fp reflects initial orientation (see 10.1016/j.actamat.2006.01.005) phase_mechanical_Fp0(ph)%data(1:3,1:3,en) = phase_mechanical_Fp0(ph)%data(1:3,1:3,en) & / math_det33(phase_mechanical_Fp0(ph)%data(1:3,1:3,en))**(1.0_pReal/3.0_pReal) phase_mechanical_Fi0(ph)%data(1:3,1:3,en) = math_I3 diff --git a/src/phase_mechanical_plastic_nonlocal.f90 b/src/phase_mechanical_plastic_nonlocal.f90 index 8eb90f547..03376721f 100644 --- a/src/phase_mechanical_plastic_nonlocal.f90 +++ b/src/phase_mechanical_plastic_nonlocal.f90 @@ -1450,8 +1450,8 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e) elseif (prm%chi_GB >= 0.0_pReal) then !* GRAIN BOUNDARY ! !* fixed transmissivity for adjacent ips with different texture (only if explicitly given in material.config) - if (any(dNeq(material_orientation0(1,ph,en)%asQuaternion(), & - material_orientation0(1,neighbor_phase,neighbor_me)%asQuaternion())) .and. & + if (any(dNeq(phase_orientation0(ph)%data(en)%asQuaternion(), & + phase_orientation0(neighbor_phase)%data(neighbor_me)%asQuaternion())) .and. & (.not. phase_localPlasticity(neighbor_phase))) & forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = sqrt(prm%chi_GB) else