needed for fixed non-Schmid behavior

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
+  integer,     dimension(:),              intent(in) :: Ntwin                                       !< number of active twin systems per family
-  character(len=*),                         intent(in) :: lattice                                   !< Bravais lattice (Pearson symbol)
+  character(len=*),                       intent(in) :: lattice                                     !< Bravais lattice (Pearson symbol)
-  real(pREAL),                              intent(in) :: cOverA                                    !< c/a ratio
+  real(pREAL),                            intent(in) :: cOverA                                      !< c/a ratio
-  real(pREAL),     dimension(3,3,sum(Ntwin))           :: SchmidMatrix
+  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
+  integer,     dimension(:),                  intent(in) :: Ncleavage                               !< number of active cleavage systems per family
-  character(len=*),                               intent(in) :: lattice                             !< Bravais lattice (Pearson symbol)
+  character(len=*),                           intent(in) :: lattice                                 !< Bravais lattice (Pearson symbol)
-  real(pREAL),                                    intent(in) :: cOverA                              !< c/a ratio
+  real(pREAL),                                intent(in) :: cOverA                                  !< c/a ratio
-  real(pREAL),     dimension(3,3,3,sum(Ncleavage))           :: SchmidMatrix
+  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) :: &

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

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) :: what
-  integer,       dimension(:), intent(in) :: how
+  integer, dimension(:), intent(in) :: how
-  real(pREAL), dimension(sum(how)) ::  math_expand
+  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_)) &