From 0396332b7b69c7b25ec76923715392d23246d602 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 14 Sep 2023 22:43:50 +0200 Subject: [PATCH] needed for fixed non-Schmid behavior --- src/crystal.f90 | 18 ++++++------- src/grid/spectral_utilities.f90 | 2 ++ src/math.f90 | 48 +++++++++++++++++++++++++++------ 3 files changed, 51 insertions(+), 17 deletions(-) diff --git a/src/crystal.f90 b/src/crystal.f90 index 7bf9dc2f7..437e45227 100644 --- a/src/crystal.f90 +++ b/src/crystal.f90 @@ -1445,10 +1445,10 @@ end function crystal_SchmidMatrix_slip !-------------------------------------------------------------------------------------------------- function crystal_SchmidMatrix_twin(Ntwin,lattice,cOverA) result(SchmidMatrix) - integer, dimension(:), intent(in) :: Ntwin !< number of active twin systems per family - character(len=*), intent(in) :: lattice !< Bravais lattice (Pearson symbol) - real(pREAL), intent(in) :: cOverA !< c/a ratio - real(pREAL), dimension(3,3,sum(Ntwin)) :: SchmidMatrix + integer, dimension(:), intent(in) :: Ntwin !< number of active twin systems per family + character(len=*), intent(in) :: lattice !< Bravais lattice (Pearson symbol) + real(pREAL), intent(in) :: cOverA !< c/a ratio + real(pREAL), dimension(3,3,sum(Ntwin)) :: SchmidMatrix real(pREAL), dimension(3,3,sum(Ntwin)) :: coordinateSystem real(pREAL), dimension(:,:), allocatable :: twinSystems @@ -1521,10 +1521,10 @@ end function crystal_SchmidMatrix_trans !-------------------------------------------------------------------------------------------------- function crystal_SchmidMatrix_cleavage(Ncleavage,lattice,cOverA) result(SchmidMatrix) - integer, dimension(:), intent(in) :: Ncleavage !< number of active cleavage systems per family - character(len=*), intent(in) :: lattice !< Bravais lattice (Pearson symbol) - real(pREAL), intent(in) :: cOverA !< c/a ratio - real(pREAL), dimension(3,3,3,sum(Ncleavage)) :: SchmidMatrix + integer, dimension(:), intent(in) :: Ncleavage !< number of active cleavage systems per family + character(len=*), intent(in) :: lattice !< Bravais lattice (Pearson symbol) + real(pREAL), intent(in) :: cOverA !< c/a ratio + real(pREAL), dimension(3,3,3,sum(Ncleavage)) :: SchmidMatrix real(pREAL), dimension(3,3,sum(Ncleavage)) :: coordinateSystem real(pREAL), dimension(:,:), allocatable :: cleavageSystems @@ -1904,7 +1904,7 @@ end function buildInteraction !-------------------------------------------------------------------------------------------------- function buildCoordinateSystem(active,potential,system,lattice,cOverA) result(coordinateSystem) - integer, dimension(:), intent(in) :: & + integer, dimension(:), intent(in) :: & active, & !< # of active systems per family potential !< # of potential systems per family real(pREAL), dimension(:,:), intent(in) :: & diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index 03d945bee..de9518b8c 100644 --- a/src/grid/spectral_utilities.f90 +++ b/src/grid/spectral_utilities.f90 @@ -674,6 +674,7 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C) logical :: errmatinv character(len=pSTRLEN):: formatString + mask_stressVector = .not. reshape(transpose(mask_stress), [9]) size_reduced = count(mask_stressVector) if (size_reduced > 0) then @@ -696,6 +697,7 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C) write(formatString, '(i2)') size_reduced formatString = '(/,1x,a,/,'//trim(formatString)//'('//trim(formatString)//'(2x,es9.2,1x)/))' print trim(formatString), 'C * S (load) ', transpose(matmul(c_reduced,s_reduced)) + print trim(formatString), 'C (load) ', transpose(c_reduced) print trim(formatString), 'S (load) ', transpose(s_reduced) if (errmatinv) error stop 'matrix inversion error' end if diff --git a/src/math.f90 b/src/math.f90 index 4122f5170..24141f4e9 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -24,6 +24,12 @@ module math implicit none(type,external) public + + interface math_expand + module procedure math_expand_int + module procedure math_expand_real + end interface math_expand + #if __INTEL_COMPILER >= 1900 ! do not make use of associated entities available to other modules private :: & @@ -136,7 +142,7 @@ end subroutine math_init pure recursive subroutine math_sort(a, istart, iend, sortDim) integer, dimension(:,:), intent(inout) :: a - integer, intent(in),optional :: istart,iend, sortDim + integer, optional, intent(in) :: istart,iend, sortDim integer :: ipivot,s,e,d @@ -199,11 +205,11 @@ end subroutine math_sort !> @details takes a set of numbers (a,b,c,...) and corresponding multiples (x,y,z,...) !> to return a vector of x times a, y times b, z times c, ... !-------------------------------------------------------------------------------------------------- -pure function math_expand(what,how) +pure function math_expand_int(what,how) - real(pREAL), dimension(:), intent(in) :: what - integer, dimension(:), intent(in) :: how - real(pREAL), dimension(sum(how)) :: math_expand + integer, dimension(:), intent(in) :: what + integer, dimension(:), intent(in) :: how + integer, dimension(sum(how)) :: math_expand_int integer :: i @@ -211,10 +217,33 @@ pure function math_expand(what,how) if (sum(how) == 0) return do i = 1, size(how) - math_expand(sum(how(1:i-1))+1:sum(how(1:i))) = what(mod(i-1,size(what))+1) + math_expand_int(sum(how(1:i-1))+1:sum(how(1:i))) = what(mod(i-1,size(what))+1) end do -end function math_expand +end function math_expand_int + + +!-------------------------------------------------------------------------------------------------- +!> @brief vector expansion +!> @details takes a set of numbers (a,b,c,...) and corresponding multiples (x,y,z,...) +!> to return a vector of x times a, y times b, z times c, ... +!-------------------------------------------------------------------------------------------------- +pure function math_expand_real(what,how) + + real(pREAL), dimension(:), intent(in) :: what + integer, dimension(:), intent(in) :: how + real(pREAL), dimension(sum(how)) :: math_expand_real + + integer :: i + + + if (sum(how) == 0) return + + do i = 1, size(how) + math_expand_real(sum(how(1:i-1))+1:sum(how(1:i))) = what(mod(i-1,size(what))+1) + end do + +end function math_expand_real !-------------------------------------------------------------------------------------------------- @@ -1309,7 +1338,10 @@ subroutine math_selfTest() if (any(abs([1.0_pREAL,2.0_pREAL,2.0_pREAL,1.0_pREAL,1.0_pREAL,1.0_pREAL] - & math_expand([1.0_pREAL,2.0_pREAL],[1,2,3])) > tol_math_check)) & - error stop 'math_expand [1,2] by [1,2,3] => [1,2,2,1,1,1]' + error stop 'math_expand_real [1,2] by [1,2,3] => [1,2,2,1,1,1]' + + if (any(abs([1,2,2,1,1,1] - math_expand([1,2],[1,2,3])) /= 0)) & + error stop 'math_expand_int [1,2] by [1,2,3] => [1,2,2,1,1,1]' call math_sort(sort_in_,1,3,2) if (any(sort_in_ /= sort_out_)) &