diff --git a/src/material.f90 b/src/material.f90 index c62442fcd..75ddceba5 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -20,7 +20,7 @@ module material type(Rotation), dimension(:), allocatable :: data end type - type(tRotationContainer), dimension(:), allocatable :: material_orientation0_2 + type(tRotationContainer), dimension(:), allocatable :: material_orientation0 integer, dimension(:), allocatable, public, protected :: & homogenization_Nconstituents !< number of grains in each homogenization @@ -46,7 +46,7 @@ module material public :: & tRotationContainer, & - material_orientation0_2, & + material_orientation0, & material_init contains @@ -161,14 +161,15 @@ subroutine parse() enddo - allocate(material_orientation0_2(materials%length)) + allocate(material_orientation0(materials%length)) do ma = 1, materials%length material => materials%get(ma) constituents => material%get('constituents') - allocate(material_orientation0_2(ma)%data(constituents%length)) + allocate(material_orientation0(ma)%data(constituents%length)) do co = 1, constituents%length - call material_orientation0_2(ma)%data(co)%fromQuaternion(constituent%get_as1dFloat('O',requiredSize=4)) + constituent => constituents%get(co) + call material_orientation0(ma)%data(co)%fromQuaternion(constituent%get_as1dFloat('O',requiredSize=4)) enddo enddo diff --git a/src/phase.f90 b/src/phase.f90 index 9d4a423d0..9dd8df094 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -23,9 +23,6 @@ module phase phase_orientation0, & phase_orientation - type(rotation), dimension(:,:,:), allocatable :: & - crystallite_orientation !< current orientation - type :: tTensorContainer real(pReal), dimension(:,:,:), allocatable :: data end type @@ -259,8 +256,7 @@ module phase ph, & i, & e - type(rotation), dimension(1,discretization_nIPs,discretization_Nelems), intent(in) :: & - orientation !< crystal orientation + type(tRotationContainer), dimension(:), intent(in) :: orientation end subroutine plastic_nonlocal_updateCompatibility module subroutine plastic_dependentState(co,ip,el) @@ -378,7 +374,7 @@ subroutine phase_init 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) + phase_orientation0(ph)%data(material_phaseEntry(co,ce)) = material_orientation0(ma)%data(co) enddo enddo @@ -523,8 +519,6 @@ subroutine crystallite_init() iMax = discretization_nIPs eMax = discretization_Nelems - allocate(crystallite_orientation(cMax,iMax,eMax)) - num_crystallite => config_numerics%get('crystallite',defaultVal=emptyDict) num%subStepMinCryst = num_crystallite%get_asFloat ('subStepMin', defaultVal=1.0e-3_pReal) @@ -594,13 +588,16 @@ subroutine crystallite_orientations(co,ip,el) ip, & !< counter in integration point loop el !< counter in element loop + integer :: ph, en - call crystallite_orientation(co,ip,el)%fromMatrix(transpose(math_rotationalPart(& - mechanical_F_e(material_phaseAt(co,el),material_phaseMemberAt(co,ip,el))))) + + ph = material_phaseID(co,(el-1)*discretization_nIPs + ip) + en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip) + + call phase_orientation(ph)%data(en)%fromMatrix(transpose(math_rotationalPart(mechanical_F_e(ph,en)))) if (plasticState(material_phaseAt(1,el))%nonlocal) & - call plastic_nonlocal_updateCompatibility(crystallite_orientation, & - material_phaseAt(1,el),ip,el) + call plastic_nonlocal_updateCompatibility(phase_orientation,material_phaseAt(1,el),ip,el) end subroutine crystallite_orientations diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index de23d39fa..ed3a737b1 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -248,18 +248,14 @@ module subroutine mechanical_init(materials,phases) phase => phases%get(ph) mech => phase%get('mechanical') #if defined(__GFORTRAN__) - output_constituent(ph)%label = output_as1dString(mech) + output_constituent(ph)%label = output_as1dString(mech) #else - output_constituent(ph)%label = mech%get_as1dString('output',defaultVal=emptyStringArray) + output_constituent(ph)%label = mech%get_as1dString('output',defaultVal=emptyStringArray) #endif enddo - 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)) - - ph = material_phaseID(co,(el-1)*discretization_nIPs + ip) - en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip) + do ph = 1, phases%length + do en = 1, count(material_phaseID == ph) 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) & @@ -268,13 +264,12 @@ module subroutine mechanical_init(materials,phases) phase_mechanical_F0(ph)%data(1:3,1:3,en) = math_I3 phase_mechanical_Fe(ph)%data(1:3,1:3,en) = math_inv33(matmul(phase_mechanical_Fi0(ph)%data(1:3,1:3,en), & - phase_mechanical_Fp0(ph)%data(1:3,1:3,en))) ! assuming that euler angles are given in internal strain free configuration - phase_mechanical_Fp(ph)%data(1:3,1:3,en) = phase_mechanical_Fp0(ph)%data(1:3,1:3,en) - phase_mechanical_Fi(ph)%data(1:3,1:3,en) = phase_mechanical_Fi0(ph)%data(1:3,1:3,en) - phase_mechanical_F(ph)%data(1:3,1:3,en) = phase_mechanical_F0(ph)%data(1:3,1:3,en) - + phase_mechanical_Fp0(ph)%data(1:3,1:3,en))) ! assuming that euler angles are given in internal strain free configuration enddo - enddo; enddo + phase_mechanical_Fp(ph)%data = phase_mechanical_Fp0(ph)%data + phase_mechanical_Fi(ph)%data = phase_mechanical_Fi0(ph)%data + phase_mechanical_F(ph)%data = phase_mechanical_F0(ph)%data + enddo ! initialize elasticity @@ -903,80 +898,71 @@ subroutine crystallite_results(group,ph) character(len=:), allocatable :: structureLabel - call results_closeGroup(results_addGroup(group//'/mechanical')) + call results_closeGroup(results_addGroup(group//'/mechanical')) - do ou = 1, size(output_constituent(ph)%label) + do ou = 1, size(output_constituent(ph)%label) - select case (output_constituent(ph)%label(ou)) - case('F') - call results_writeDataset(group//'/mechanical/',phase_mechanical_F(ph)%data,'F',& - 'deformation gradient','1') - case('F_e') - call results_writeDataset(group//'/mechanical/',phase_mechanical_Fe(ph)%data,'F_e',& - 'elastic deformation gradient','1') - case('F_p') - call results_writeDataset(group//'/mechanical/',phase_mechanical_Fp(ph)%data,'F_p', & - 'plastic deformation gradient','1') - case('F_i') - call results_writeDataset(group//'/mechanical/',phase_mechanical_Fi(ph)%data,'F_i', & - 'inelastic deformation gradient','1') - case('L_p') - call results_writeDataset(group//'/mechanical/',phase_mechanical_Lp(ph)%data,'L_p', & - 'plastic velocity gradient','1/s') - case('L_i') - call results_writeDataset(group//'/mechanical/',phase_mechanical_Li(ph)%data,'L_i', & - 'inelastic velocity gradient','1/s') - case('P') - call results_writeDataset(group//'/mechanical/',phase_mechanical_P(ph)%data,'P', & - 'first Piola-Kirchhoff stress','Pa') - case('S') - call results_writeDataset(group//'/mechanical/',phase_mechanical_S(ph)%data,'S', & - '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//'/mechanical',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//'/mechanical/'//output_constituent(ph)%label(ou)) - end select - enddo + select case (output_constituent(ph)%label(ou)) + case('F') + call results_writeDataset(group//'/mechanical/',phase_mechanical_F(ph)%data,'F',& + 'deformation gradient','1') + case('F_e') + call results_writeDataset(group//'/mechanical/',phase_mechanical_Fe(ph)%data,'F_e',& + 'elastic deformation gradient','1') + case('F_p') + call results_writeDataset(group//'/mechanical/',phase_mechanical_Fp(ph)%data,'F_p', & + 'plastic deformation gradient','1') + case('F_i') + call results_writeDataset(group//'/mechanical/',phase_mechanical_Fi(ph)%data,'F_i', & + 'inelastic deformation gradient','1') + case('L_p') + call results_writeDataset(group//'/mechanical/',phase_mechanical_Lp(ph)%data,'L_p', & + 'plastic velocity gradient','1/s') + case('L_i') + call results_writeDataset(group//'/mechanical/',phase_mechanical_Li(ph)%data,'L_i', & + 'inelastic velocity gradient','1/s') + case('P') + call results_writeDataset(group//'/mechanical/',phase_mechanical_P(ph)%data,'P', & + 'first Piola-Kirchhoff stress','Pa') + case('S') + call results_writeDataset(group//'/mechanical/',phase_mechanical_S(ph)%data,'S', & + '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(phase_orientation(ph)%data) + call results_writeDataset(group//'/mechanical',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//'/mechanical/'//output_constituent(ph)%label(ou)) + end select + enddo contains !-------------------------------------------------------------------------------------------------- -!> @brief select rotations for output +!> @brief Convert orientation for output: ToDo: implement in HDF5/results !-------------------------------------------------------------------------------------------------- - function select_rotations(dataset,ph) + function select_rotations(dataset) - integer, intent(in) :: ph - type(rotation), dimension(:,:,:), intent(in) :: dataset - real(pReal), dimension(4,count(material_phaseID==ph)) :: select_rotations - integer :: el,ip,co,j + type(rotation), dimension(:), intent(in) :: dataset + real(pReal), dimension(4,size(dataset,1)) :: select_rotations + integer :: en - 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 + do en = 1, size(dataset,1) + select_rotations(:,en) = dataset(en)%asQuaternion() enddo end function select_rotations diff --git a/src/phase_mechanical_plastic_nonlocal.f90 b/src/phase_mechanical_plastic_nonlocal.f90 index 03376721f..0eccf767e 100644 --- a/src/phase_mechanical_plastic_nonlocal.f90 +++ b/src/phase_mechanical_plastic_nonlocal.f90 @@ -1394,7 +1394,7 @@ end function rhoDotFlux !-------------------------------------------------------------------------------------------------- module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e) - type(rotation), dimension(1,discretization_nIPs,discretization_Nelems), intent(in) :: & + type(tRotationContainer), dimension(:), intent(in) :: & orientation ! crystal orientation integer, intent(in) :: & ph, & @@ -1464,7 +1464,7 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e) !* the number of compatible slip systems is minimized with the sum of the original compatibility values exceeding one. !* Finally the smallest compatibility value is decreased until the sum is exactly equal to one. !* All values below the threshold are set to zero. - mis = orientation(1,i,e)%misorientation(orientation(1,neighbor_i,neighbor_e)) + mis = orientation(ph)%data(en)%misorientation(orientation(neighbor_phase)%data(neighbor_me)) mySlipSystems: do s1 = 1,ns neighborSlipSystems: do s2 = 1,ns my_compatibility(1,s2,s1,n) = math_inner(prm%slip_normal(1:3,s1), &