I don't like loops

use language features and helper functions for shorter code

src/grid/spectral_utilities.f90

@@ -688,9 +688,10 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C)
   type(rotation), intent(in)                     :: rot_BC                                          !< rotation of load frame
   logical,        intent(in), dimension(3,3)     :: mask_stress                                     !< mask of stress BC
 
-  integer :: j, k, m, n
+  integer :: i, j
   logical, dimension(9)   :: mask_stressVector
-  real(pReal), dimension(9,9) :: temp99_Real
+  logical, dimension(9,9) :: mask
+  real(pReal), dimension(9,9) :: temp99_real
   integer :: size_reduced = 0
   real(pReal),              dimension(:,:), allocatable ::  &
     s_reduced, &                                                                                    !< reduced compliance matrix (depending on number of stress BC)
@@ -702,10 +703,7 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C)
   mask_stressVector = reshape(transpose(mask_stress), [9])
   size_reduced = count(mask_stressVector)
   if(size_reduced > 0) then
-    allocate (c_reduced(size_reduced,size_reduced), source =0.0_pReal)
+    temp99_real = math_3333to99(rot_BC%rotate(C))
-    allocate (s_reduced(size_reduced,size_reduced), source =0.0_pReal)
-    allocate (sTimesC(size_reduced,size_reduced),   source =0.0_pReal)
-    temp99_Real = math_3333to99(rot_BC%rotTensor4(C))
 
     if(debugGeneral) then
       write(6,'(/,a)') ' ... updating masked compliance ............................................'
@@ -713,42 +711,21 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C)
                                                    transpose(temp99_Real)*1.0e-9_pReal
       flush(6)
     endif
-    k = 0                                                                                           ! calculate reduced stiffness
-    do n = 1,9
-      if(mask_stressVector(n)) then
-        k = k + 1
-        j = 0
-        do m = 1,9
-          if(mask_stressVector(m)) then
-            j = j + 1
-            c_reduced(k,j) = temp99_Real(n,m)
-    endif; enddo; endif; enddo
 
+    do i = 1,9; do j = 1,9
+      mask(i,j) = mask_stressVector(i) .and. mask_stressVector(j)
+    enddo; enddo
+    c_reduced = reshape(pack(temp99_Real,mask),[size_reduced,size_reduced])
+
+    allocate(s_reduced,mold = c_reduced)
     call math_invert(s_reduced, errmatinv, c_reduced)                                               ! invert reduced stiffness
     if (any(IEEE_is_NaN(s_reduced))) errmatinv = .true.
     if (errmatinv) call IO_error(error_ID=400,ext_msg='utilities_maskedCompliance')
-    temp99_Real = 0.0_pReal                                                                         ! fill up compliance with zeros
-     k = 0
-     do n = 1,9
-       if(mask_stressVector(n)) then
-         k = k + 1
-         j = 0
-         do m = 1,9
-           if(mask_stressVector(m)) then
-             j = j + 1
-             temp99_Real(n,m) = s_reduced(k,j)
-    endif; enddo; endif; enddo
 
 !--------------------------------------------------------------------------------------------------
 ! check if inversion was successful
     sTimesC = matmul(c_reduced,s_reduced)
-    do m=1, size_reduced
+    errmatinv = errmatinv .or. any(dNeq(sTimesC,math_identity2nd(size_reduced),1.0e-12_pReal))
-      do n=1, size_reduced
-        errmatinv = errmatinv &
-               .or. (m==n .and. abs(sTimesC(m,n)-1.0_pReal) > 1.0e-12_pReal) &                      ! diagonal elements of S*C should be 1
-               .or. (m/=n .and. abs(sTimesC(m,n))           > 1.0e-12_pReal)                        ! off-diagonal elements of S*C should be 0
-      enddo
-    enddo
     if (debugGeneral .or. errmatinv) then
       write(formatString, '(i2)') size_reduced
       formatString = '(/,a,/,'//trim(formatString)//'('//trim(formatString)//'(2x,es9.2,1x)/))'
@@ -757,15 +734,18 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C)
       write(6,trim(formatString),advance='no') ' S (load) ', transpose(s_reduced)
       if(errmatinv) call IO_error(error_ID=400,ext_msg='utilities_maskedCompliance')
     endif
+    temp99_real = reshape(unpack(reshape(s_reduced,[size_reduced**2]),reshape(mask,[81]),0.0_pReal),[9,9])
   else
     temp99_real = 0.0_pReal
   endif
+
+  utilities_maskedCompliance = math_99to3333(temp99_Real)
+
   if(debugGeneral) then
     write(6,'(/,a,/,9(9(2x,f10.5,1x)/),/)',advance='no') &
       ' Masked Compliance (load) * GPa =', transpose(temp99_Real)*1.0e9_pReal
     flush(6)
   endif
-  utilities_maskedCompliance = math_99to3333(temp99_Real)
 
 end function utilities_maskedCompliance
 
@@ -862,7 +842,7 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,&
   write(6,'(/,a,/,3(3(2x,f12.4,1x)/))',advance='no') ' Piola--Kirchhoff stress (lab) / MPa =',&
                                                       transpose(P_av)*1.e-6_pReal
   if(present(rotation_BC)) &
-    P_av = rotation_BC%rotTensor2(P_av)
+    P_av = rotation_BC%rotate(P_av)
   write(6,'(/,a,/,3(3(2x,f12.4,1x)/))',advance='no') ' Piola--Kirchhoff stress       / MPa =',&
                                                       transpose(P_av)*1.e-6_pReal
   flush(6)