From 6701af6425148ea2c1c7d621a98abb51ca4b9adf Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Sun, 1 Mar 2020 23:00:06 +0100
Subject: [PATCH] better readable, tested, and following standard notation

---
 src/homogenization_mech_RGC.f90 |  4 ++--
 src/math.f90                    | 32 +++++++++++++++++++++-----------
 2 files changed, 23 insertions(+), 13 deletions(-)

diff --git a/src/homogenization_mech_RGC.f90 b/src/homogenization_mech_RGC.f90
index 1f1d9d6ec..2152211b7 100644
--- a/src/homogenization_mech_RGC.f90
+++ b/src/homogenization_mech_RGC.f90
@@ -669,7 +669,7 @@ module procedure mech_RGC_updateState
        nDef = 0.0_pReal
        do i = 1,3; do j = 1,3
          do k = 1,3; do l = 1,3
-           nDef(i,j) = nDef(i,j) - nVect(k)*gDef(i,l)*math_civita(j,k,l)                            ! compute the interface mismatch tensor from the jump of deformation gradient
+           nDef(i,j) = nDef(i,j) - nVect(k)*gDef(i,l)*math_LeviCivita(j,k,l)                        ! compute the interface mismatch tensor from the jump of deformation gradient
          enddo; enddo
          nDefNorm = nDefNorm + nDef(i,j)**2.0_pReal                                                 ! compute the norm of the mismatch tensor
        enddo; enddo
@@ -689,7 +689,7 @@ module procedure mech_RGC_updateState
          rPen(i,j,iGrain) = rPen(i,j,iGrain) + 0.5_pReal*(muGrain*bgGrain + muGNghb*bgGNghb)*prm%xiAlpha &
                                                 *surfCorr(abs(intFace(1)))/prm%dAlpha(abs(intFace(1))) &
                                                 *cosh(prm%ciAlpha*nDefNorm) &
-                                                *0.5_pReal*nVect(l)*nDef(i,k)/nDefNorm*math_civita(k,l,j) &
+                                                *0.5_pReal*nVect(l)*nDef(i,k)/nDefNorm*math_LeviCivita(k,l,j) &
                                                 *tanh(nDefNorm/xSmoo_RGC)
        enddo; enddo;enddo; enddo
      enddo interfaceLoop
diff --git a/src/math.f90 b/src/math.f90
index 80c3ed54e..499190fac 100644
--- a/src/math.f90
+++ b/src/math.f90
@@ -271,31 +271,30 @@ end function math_identity4th
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief permutation tensor e_ijk used for computing cross product of two tensors
+!> @brief permutation tensor e_ijk
 ! e_ijk =  1 if even permutation of ijk
 ! e_ijk = -1 if odd permutation of ijk
 ! e_ijk =  0 otherwise
 !--------------------------------------------------------------------------------------------------
-real(pReal) pure function math_civita(i,j,k)
+real(pReal) pure function math_LeviCivita(i,j,k)
 
   integer, intent(in) :: i,j,k
 
-  math_civita = 0.0_pReal
-  if (((i == 1).and.(j == 2).and.(k == 3)) .or. &
-      ((i == 2).and.(j == 3).and.(k == 1)) .or. &
-      ((i == 3).and.(j == 1).and.(k == 2))) math_civita = 1.0_pReal
-  if (((i == 1).and.(j == 3).and.(k == 2)) .or. &
-      ((i == 2).and.(j == 1).and.(k == 3)) .or. &
-      ((i == 3).and.(j == 2).and.(k == 1))) math_civita = -1.0_pReal
+  if     (all([i,j,k] == [1,2,3]) .or. all([i,j,k] == [2,3,1]) .or. all([i,j,k] == [3,1,2])) then
+    math_LeviCivita = +1.0_pReal
+  elseif (all([i,j,k] == [3,2,1]) .or. all([i,j,k] == [2,1,3]) .or. all([i,j,k] == [1,3,2])) then
+    math_LeviCivita = -1.0_pReal
+  else
+    math_LeviCivita =  0.0_pReal
+  endif
 
-end function math_civita
+end function math_LeviCivita
 
 
 !--------------------------------------------------------------------------------------------------
 !> @brief kronecker delta function d_ij
 ! d_ij = 1 if i = j
 ! d_ij = 0 otherwise
-! inspired by http://fortraninacworld.blogspot.de/2012/12/ternary-operator.html
 !--------------------------------------------------------------------------------------------------
 real(pReal) pure function math_delta(i,j)
 
@@ -1309,6 +1308,7 @@ subroutine unitTest
     sort_out_  = reshape([-1,-1,  +1,+5,  +5,+6,  +3,-2],[2,4])
 
   integer,     dimension(5) :: range_out_ = [1,2,3,4,5]
+  integer,     dimension(3) :: ijk
 
   real(pReal)                 :: det
   real(pReal), dimension(3)   :: v3_1,v3_2,v3_3,v3_4
@@ -1425,6 +1425,16 @@ subroutine unitTest
   if(math_binomial(49,6) /= 13983816) &
     call IO_error(0,ext_msg='math_binomial')
 
+  ijk = cshift([1,2,3],int(r*1.0e2_pReal))
+  if(dNeq(math_LeviCivita(ijk(1),ijk(2),ijk(3)),+1.0_pReal)) &
+    call IO_error(0,ext_msg='math_LeviCivita(even)')
+  ijk = cshift([3,2,1],int(r*2.0e2_pReal))
+  if(dNeq(math_LeviCivita(ijk(1),ijk(2),ijk(3)),-1.0_pReal)) &
+    call IO_error(0,ext_msg='math_LeviCivita(odd)')
+  ijk = cshift([2,2,1],int(r*2.0e2_pReal))
+  if(dNeq0(math_LeviCivita(ijk(1),ijk(2),ijk(3))))&
+    call IO_error(0,ext_msg='math_LeviCivita')
+
 end subroutine unitTest
 
 end module math