From dfe6d0a1957f516fe2c1d8d38329224ede53dcd9 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Thu, 18 Nov 2021 21:06:38 +0100
Subject: [PATCH] more support for Voigt notation

---
 src/math.f90                     | 31 ++++++++++++++++++++++++++++++-
 src/phase_mechanical_elastic.f90 |  1 -
 2 files changed, 30 insertions(+), 2 deletions(-)

diff --git a/src/math.f90 b/src/math.f90
index 31bd3b952..2515c64d7 100644
--- a/src/math.f90
+++ b/src/math.f90
@@ -854,12 +854,13 @@ end function math_66toSym3333
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief convert 66 Voigt matrix into symmetric 3x3x3x3 matrix
+!> @brief Convert 6x6 Voigt matrix into symmetric 3x3x3x3 matrix.
 !--------------------------------------------------------------------------------------------------
 pure function math_Voigt66to3333(m66)
 
   real(pReal), dimension(3,3,3,3) :: math_Voigt66to3333
   real(pReal), dimension(6,6), intent(in) :: m66                                                    !< 6x6 matrix
+
   integer :: i,j
 
 
@@ -873,6 +874,31 @@ pure function math_Voigt66to3333(m66)
 end function math_Voigt66to3333
 
 
+!--------------------------------------------------------------------------------------------------
+!> @brief Convert symmetric 3x3x3x3 matrix into 6x6 Voigt matrix.
+!--------------------------------------------------------------------------------------------------
+pure function math_3333toVoigt66(m3333)
+
+  real(pReal), dimension(6,6)                 :: math_3333toVoigt66
+  real(pReal), dimension(3,3,3,3), intent(in) :: m3333                                              !< symmetric 3x3x3x3 matrix (no internal check)
+
+  integer :: i,j
+
+
+#ifndef __INTEL_COMPILER
+  do concurrent(i=1:6, j=1:6)
+    math_3333toVoigt66(i,j) = m3333(MAPVOIGT(1,i),MAPVOIGT(2,i),MAPVOIGT(1,j),MAPVOIGT(2,j))
+  end do
+#else
+  do i=1,6; do j=1,6
+    math_3333toVoigt66(i,j) = m3333(MAPVOIGT(1,i),MAPVOIGT(2,i),MAPVOIGT(1,j),MAPVOIGT(2,j))
+  end do; end do
+#endif
+
+end function math_3333toVoigt66
+
+
+
 !--------------------------------------------------------------------------------------------------
 !> @brief draw a random sample from Gauss variable
 !--------------------------------------------------------------------------------------------------
@@ -1261,6 +1287,9 @@ subroutine selfTest
   if(any(dNeq(math_sym3333to66(math_66toSym3333(t66)),t66,1.0e-15_pReal))) &
     error stop 'math_sym3333to66/math_66toSym3333'
 
+  if(any(dNeq(math_3333toVoigt66(math_Voigt66to3333(t66)),t66,1.0e-15_pReal))) &
+    error stop 'math_3333toVoigt66/math_Voigt66to3333'
+
   call random_number(v6)
   if(any(dNeq0(math_6toSym33(v6) - math_symmetric33(math_6toSym33(v6))))) &
     error stop 'math_symmetric33'
diff --git a/src/phase_mechanical_elastic.f90 b/src/phase_mechanical_elastic.f90
index 8ab6884f9..0f7394e86 100644
--- a/src/phase_mechanical_elastic.f90
+++ b/src/phase_mechanical_elastic.f90
@@ -167,7 +167,6 @@ end subroutine phase_hooke_SandItsTangents
 
 !--------------------------------------------------------------------------------------------------
 !> @brief returns the homogenized elasticity matrix
-!> ToDo: homogenizedC66 would be more consistent
 !--------------------------------------------------------------------------------------------------
 module function phase_homogenizedC66(ph,en) result(C)