From 372536d57ed0f571ba334f6c75795eb7c2f7e06a Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Sun, 22 Sep 2019 06:53:03 -0700
Subject: [PATCH] unit test for rotation class

---
 src/CPFEM.f90          |  2 ++
 src/CPFEM2.f90         |  2 ++
 src/homogenization.f90 |  2 +-
 src/rotations.f90      | 60 +++++++++++++++++++++++++++++++++++++++---
 4 files changed, 62 insertions(+), 4 deletions(-)

diff --git a/src/CPFEM.f90 b/src/CPFEM.f90
index 42ebf4459..6649de542 100644
--- a/src/CPFEM.f90
+++ b/src/CPFEM.f90
@@ -9,6 +9,7 @@ module CPFEM
   use debug
   use FEsolving
   use math
+  use rotations
   use mesh
   use material
   use config
@@ -84,6 +85,7 @@ subroutine CPFEM_initAll(el,ip)
      call debug_init
      call config_init
      call math_init
+     call rotations_init
      call FE_init
 #ifdef DAMASK_HDF5
      call HDF5_utilities_init
diff --git a/src/CPFEM2.f90 b/src/CPFEM2.f90
index dce2bb747..7e7949ec8 100644
--- a/src/CPFEM2.f90
+++ b/src/CPFEM2.f90
@@ -10,6 +10,7 @@ module CPFEM2
   use config
   use FEsolving
   use math
+  use rotations
   use mesh
   use material
   use lattice
@@ -52,6 +53,7 @@ subroutine CPFEM_initAll
   call debug_init
   call config_init
   call math_init
+  call rotations_init
   call mesh_init
   call lattice_init
   call HDF5_utilities_init
diff --git a/src/homogenization.f90 b/src/homogenization.f90
index 64019b880..c7675ccba 100644
--- a/src/homogenization.f90
+++ b/src/homogenization.f90
@@ -801,7 +801,7 @@ subroutine homogenization_results
   integer :: p
   character(len=256) :: group
   
-  real(pReal), dimension(:,:,:), allocatable :: temp
+  !real(pReal), dimension(:,:,:), allocatable :: temp
                                              
   do p=1,size(config_name_homogenization)
     group = trim('current/materialpoint')//'/'//trim(config_name_homogenization(p))
diff --git a/src/rotations.f90 b/src/rotations.f90
index 4624f965c..1a0a8bca2 100644
--- a/src/rotations.f90
+++ b/src/rotations.f90
@@ -78,10 +78,24 @@ module rotations
       procedure, public  :: misorientation
   end type rotation
   
+  public :: &
+    rotations_init
+
 
 contains
 
 
+!--------------------------------------------------------------------------------------------------
+!> @brief doing self test
+!--------------------------------------------------------------------------------------------------
+subroutine rotations_init()
+
+  write(6,'(/,a)')   ' <<<+-  rotations init  -+>>>'
+  call unitTest()
+
+end subroutine rotations_init
+
+
 !---------------------------------------------------------------------------------------------------
 ! Return rotation in different represenations
 !---------------------------------------------------------------------------------------------------
@@ -553,7 +567,7 @@ function om2ax(om) result(ax)
   real(pReal), dimension(3)        :: Wr, Wi
   real(pReal), dimension((64+2)*3) :: work
   real(pReal), dimension(3,3)      :: VR, devNull, om_
-    integer                        :: ierr, i
+  integer                          :: ierr, i
   
   external :: dgeev
   
@@ -567,13 +581,13 @@ function om2ax(om) result(ax)
     ax(1:3) = [ 0.0, 0.0, 1.0 ]
   else
     call dgeev('N','V',3,om_,3,Wr,Wi,devNull,3,VR,3,work,size(work,1),ierr)
-    if (ierr /= 0) call IO_error(0,ext_msg='Error in om2ax DGEEV return not zero')
+    if (ierr /= 0) call IO_error(0,ext_msg='Error in om2ax: DGEEV return not zero')
 #if defined(__GFORTRAN__) &&  __GNUC__ < 9 || __INTEL_COMPILER < 1800
     i = maxloc(merge(1,0,cEq(cmplx(Wr,Wi,pReal),cmplx(1.0_pReal,0.0_pReal,pReal),tol=1.0e-14_pReal)),dim=1)
 #else
     i = findloc(cEq(cmplx(Wr,Wi,pReal),cmplx(1.0_pReal,0.0_pReal,pReal),tol=1.0e-14_pReal),.true.,dim=1) !find eigenvalue (1,0)
 #endif
-    if (i == 0) call IO_error(0,ext_msg='Error in om2ax Real eigenvalue not found')
+    if (i == 0) call IO_error(0,ext_msg='Error in om2ax Real: eigenvalue not found')
     ax(1:3) = VR(1:3,i)
     where (                dNeq0([om(2,3)-om(3,2), om(3,1)-om(1,3), om(1,2)-om(2,1)])) &
       ax(1:3) = sign(ax(1:3),-P *[om(2,3)-om(3,2), om(3,1)-om(1,3), om(1,2)-om(2,1)])
@@ -1180,4 +1194,44 @@ pure function cu2ho(cu) result(ho)
 end function cu2ho
 
 
+!--------------------------------------------------------------------------------------------------
+!> @brief check correctness of (some) rotations functions
+!--------------------------------------------------------------------------------------------------
+subroutine unitTest
+  
+  real(pReal), dimension(4)   :: qu,ro,ax
+  real(pReal), dimension(3)   :: eu,r
+  real(pReal), dimension(3,3) :: om
+  character(len=pStringLen)   :: msg
+  real    :: A,B
+  integer :: i
+
+  do i=1,5
+  
+    msg = ''
+
+    call random_number(r)
+    A = sqrt(r(3))
+    B = sqrt(1-0_pReal -r(3))
+    qu = [cos(2.0_pReal*PI*r(1))*A,&
+          sin(2.0_pReal*PI*r(2))*B,&
+          cos(2.0_pReal*PI*r(2))*B,&
+          sin(2.0_pReal*PI*r(1))*A]
+    if(qu(1)<0.0_pReal) qu = qu * (-1.0_pReal)
+    
+    if(any(dNeq(om2qu(qu2om(qu)),qu,1.0e-10_pReal))) msg = trim(msg)//'om2qu/qu2om,'
+    if(any(dNeq(eu2qu(qu2eu(qu)),qu,1.0e-10_pReal))) msg = trim(msg)//'eu2qu/qu2eu,'
+    if(any(dNeq(ax2qu(qu2ax(qu)),qu,1.0e-10_pReal))) msg = trim(msg)//'ax2qu/qu2ax,'
+    if(any(dNeq(ro2qu(qu2ro(qu)),qu,1.0e-10_pReal))) msg = trim(msg)//'ro2qu/qu2ro,'
+    if(any(dNeq(ho2qu(qu2ho(qu)),qu,1.0e-7_pReal)))  msg = trim(msg)//'ho2qu/qu2ho,'
+    if(any(dNeq(cu2qu(qu2cu(qu)),qu,1.0e-7_pReal)))  msg = trim(msg)//'cu2qu/qu2cu,'
+    
+    if(len_trim(msg) /= 0) &
+      call IO_error(401,ext_msg=msg)
+
+  enddo
+      
+end subroutine unitTest
+
+
 end module rotations