Merge branch 'misc-improvements' into 'development'

fixes for a few small glitches

See merge request damask/DAMASK!503

.gitlab-ci.yml

@@ -46,7 +46,7 @@ variables:
   # ++++++++++++ PETSc ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   PETSC_GNU:       "Libraries/PETSc/3.16.1/GNU-10-OpenMPI-4.1.1"
   PETSC_INTELLLVM: "Libraries/PETSc/3.16.3/oneAPI-2022.0.1-IntelMPI-2021.5.0"
-  PETSC_INTEL:    "Libraries/PETSc/3.16.2/Intel-2022.0.1-IntelMPI-2021.5.0"
+  PETSC_INTEL:    "Libraries/PETSc/3.16.3/Intel-2022.0.1-IntelMPI-2021.5.0"
   # ++++++++++++ MSC Marc +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   MSC:            "FEM/MSC/2021.3.1"
   IntelMarc:      "Compiler/Intel/19.1.2 Libraries/IMKL/2020"
@@ -79,63 +79,63 @@ mypy:
 
 
 ###################################################################################################
-test_grid_GNU:
+grid_GNU:
   stage: compile
   script:
     - module load ${COMPILER_GNU} ${MPI_GNU} ${PETSC_GNU}
     - cd PRIVATE/testing/pytest
     - pytest -k 'compile and grid' --basetemp ${TESTROOT}/compile_grid_GNU
 
-test_mesh_GNU:
+mesh_GNU:
   stage: compile
   script:
     - module load ${COMPILER_GNU} ${MPI_GNU} ${PETSC_GNU}
     - cd PRIVATE/testing/pytest
     - pytest -k 'compile and mesh' --basetemp ${TESTROOT}/compile_mesh_GNU
 
-test_grid_GNU-64bit:
+grid_GNU-64bit:
   stage: compile
   script:
     - module load Compiler/GNU/10 Libraries/PETSc/3.16.2/64bit
     - cd PRIVATE/testing/pytest
     - pytest -k 'compile and grid' --basetemp ${TESTROOT}/compile_grid_GNU-64bit
 
-test_mesh_GNU-64bit:
+mesh_GNU-64bit:
   stage: compile
   script:
     - module load Compiler/GNU/10 Libraries/PETSc/3.16.2/64bit
     - cd PRIVATE/testing/pytest
     - pytest -k 'compile and mesh' --basetemp ${TESTROOT}/compile_mesh_GNU-64bit
 
-test_grid_IntelLLVM:
+grid_IntelLLVM:
   stage: compile
   script:
     - module load ${COMPILER_INTELLLVM} ${MPI_INTELLLVM} ${PETSC_INTELLLVM}
     - cd PRIVATE/testing/pytest
     - pytest -k 'compile and grid' --basetemp ${TESTROOT}/compile_grid_IntelLLVM
 
-test_mesh_IntelLLVM:
+mesh_IntelLLVM:
   stage: compile
   script:
     - module load ${COMPILER_INTELLLVM} ${MPI_INTELLLVM} ${PETSC_INTELLLVM}
     - cd PRIVATE/testing/pytest
     - pytest -k 'compile and mesh' --basetemp ${TESTROOT}/compile_mesh_IntelLLVM
 
-test_grid_Intel:
+grid_Intel:
   stage: compile
   script:
     - module load ${COMPILER_INTEL} ${MPI_INTEL} ${PETSC_INTEL}
     - cd PRIVATE/testing/pytest
     - pytest -k 'compile and grid' --basetemp ${TESTROOT}/compile_grid_Intel
 
-test_mesh_Intel:
+mesh_Intel:
   stage: compile
   script:
     - module load ${COMPILER_INTEL} ${MPI_INTEL} ${PETSC_INTEL}
     - cd PRIVATE/testing/pytest
     - pytest -k 'compile and mesh' --basetemp ${TESTROOT}/compile_mesh_Intel
 
-test_Marc:
+Marc_Intel:
   stage: compile
   script:
     - module load $IntelMarc $HDF5Marc $MSC
@@ -158,7 +158,7 @@ setup_mesh:
     - cmake -DDAMASK_SOLVER=MESH -DCMAKE_INSTALL_PREFIX=${TESTROOT} ${CI_PROJECT_DIR}
     - make -j2 all install
 
-compile_Marc:
+setup_Marc:
   stage: compile
   script:
     - module load $IntelMarc $HDF5Marc $MSC

PRIVATE

@@ -1 +1 @@
-Subproject commit e2cf0c5513b575f2775d2c1b9bfba8c22cd40715
+Subproject commit d570b4a1fac5b9b99d026d3ff9bf593615e22ce5

python/damask/__init__.py

@@ -21,8 +21,8 @@ from ._rotation        import Rotation         # noqa
 from ._crystal         import Crystal          # noqa
 from ._orientation     import Orientation      # noqa
 from ._table           import Table            # noqa
-from ._vtk             import VTK              # noqa
 from ._colormap        import Colormap         # noqa
+from ._vtk             import VTK              # noqa
 from ._config          import Config           # noqa
 from ._configmaterial  import ConfigMaterial   # noqa
 from ._grid            import Grid             # noqa

python/damask/_grid.py

@@ -674,7 +674,7 @@ class Grid:
 
     def show(self) -> None:
         """Show on screen."""
-        VTK.from_rectilinear_grid(self.cells,self.size,self.origin).show()
+        VTK.from_image_data(self.cells,self.size,self.origin).show()
 
 
     def add_primitive(self,

python/damask/_result.py

@@ -981,7 +981,7 @@ class Result:
             t = 'tensor'
             if o is None: o = 'fro'
         else:
-            raise ValueError(f'invalid norm order {ord}')
+            raise ValueError(f'invalid shape of {x["label"]}')
 
         return {
                 'data':  np.linalg.norm(x['data'],ord=o,axis=axis,keepdims=True),

python/damask/_rotation.py

@@ -812,8 +812,7 @@ class Rotation:
         return Rotation.from_basis(R)
 
     @staticmethod
-    def from_parallel(a,b,
+    def from_parallel(a,b):
-                      **kwargs):
         """
         Initialize from pairs of two orthogonal lattice basis vectors.
 
@@ -963,8 +962,7 @@ class Rotation:
                  N = 500,
                  degrees = True,
                  fractions = True,
-                 rng_seed = None,
+                 rng_seed = None):
-                 **kwargs):
         """
         Sample discrete values from a binned orientation distribution function (ODF).
 

python/damask/_vtk.py

@@ -442,22 +442,24 @@ class VTK:
 
         mapper = vtk.vtkDataSetMapper()
         mapper.SetInputData(self.vtk_data)
+
         actor = vtk.vtkActor()
         actor.SetMapper(mapper)
+        actor.GetProperty().SetColor(230/255,150/255,68/255)
 
         ren = vtk.vtkRenderer()
+        ren.AddActor(actor)
+        ren.SetBackground(67/255,128/255,208/255)
 
         window = vtk.vtkRenderWindow()
         window.AddRenderer(ren)
-
-        ren.AddActor(actor)
-        ren.SetBackground(0.2,0.2,0.2)
-
         window.SetSize(width,height)
+        window.SetWindowName(util.execution_stamp('VTK','show'))
 
         iren = vtk.vtkRenderWindowInteractor()
         iren.SetRenderWindow(window)
-
+        if os.name == 'posix' and 'DISPLAY' not in os.environ:
-        iren.Initialize()
+            print('Found no rendering device')
-        window.Render()
+        else:
-        iren.Start()
+            window.Render()
+            iren.Start()

python/tests/test_Colormap.py

@@ -81,6 +81,7 @@ class TestColormap:
         assert Colormap.from_predefined('strain') == Colormap.from_predefined('strain')
         assert Colormap.from_predefined('strain') != Colormap.from_predefined('stress')
         assert Colormap.from_predefined('strain',N=128) != Colormap.from_predefined('strain',N=64)
+        assert not Colormap.from_predefined('strain',N=128) == 1
 
     @pytest.mark.parametrize('low,high',[((0,0,0),(1,1,1)),
                                          ([0,0,0],[1,1,1])])

python/tests/test_Crystal.py

@@ -40,6 +40,9 @@ class TestCrystal:
                     alpha=alpha,beta=beta,gamma=gamma)
         assert np.allclose(np.eye(3),np.einsum('ik,jk',c.basis_real,c.basis_reciprocal))
 
+    def test_basis_invalid(self):
+        with pytest.raises(KeyError):
+            Crystal(family='cubic').basis_real
 
     @pytest.mark.parametrize('keyFrame,keyLattice',[('uvw','direction'),('hkl','plane'),])
     @pytest.mark.parametrize('vector',np.array([

python/tests/test_Grid.py

@@ -44,6 +44,7 @@ class TestGrid:
 
     def test_equal(self,default):
         assert default == default
+        assert not default == 42
 
     def test_repr(self,default):
         print(default)

python/tests/test_Orientation.py

@@ -364,6 +364,11 @@ class TestOrientation:
                 table.save(reference)
             assert np.allclose(P,Table.load(reference).get('Schmid'))
 
+    def test_Schmid_invalid(self):
+        with pytest.raises(KeyError):
+            Orientation(lattice='fcc').Schmid()
+
+
 ### vectorization tests ###
 
     @pytest.mark.parametrize('lattice',['hP','cI','cF']) # tI not included yet
@@ -505,3 +510,7 @@ class TestOrientation:
         for loc in np.random.randint(0,blend,(10,len(blend))):
             assert np.allclose(o[tuple(loc[:len(o.shape)])].to_pole(uvw=v[tuple(loc[-len(v.shape[:-1]):])]),
                                o.to_pole(uvw=v)[tuple(loc)])
+
+    def test_mul_invalid(self):
+        with pytest.raises(TypeError):
+            Orientation.from_random(lattice='cF')*np.ones(3)

python/tests/test_Result.py

@@ -102,6 +102,9 @@ class TestResult:
         with pytest.raises(AttributeError):
             default.view('invalid',True)
 
+    def test_add_invalid(self,default):
+        default.add_absolute('xxxx')
+
     def test_add_absolute(self,default):
         default.add_absolute('F_e')
         in_memory = np.abs(default.place('F_e'))

python/tests/test_Rotation.py

@@ -792,6 +792,11 @@ class TestRotation:
         R = Rotation.from_random(shape,rng_seed=1)
         assert R == R if shape is None else (R == R).all()
 
+    @pytest.mark.parametrize('shape',[None,5,(4,6)])
+    def test_allclose(self,shape):
+        R = Rotation.from_random(shape,rng_seed=1)
+        assert R.allclose(R)
+
     @pytest.mark.parametrize('shape',[None,5,(4,6)])
     def test_unequal(self,shape):
         R = Rotation.from_random(shape,rng_seed=1)
@@ -1124,3 +1129,7 @@ class TestRotation:
         weights_r = np.histogramdd(Eulers_r,steps,rng)[0].flatten(order='F')/N * np.sum(weights)
 
         assert np.sqrt(((weights_r - weights) ** 2).mean()) < 5
+
+    def test_mul_invalid(self):
+        with pytest.raises(TypeError):
+            Rotation.from_random()*np.ones(3)

python/tests/test_VTK.py

@@ -28,6 +28,10 @@ class TestVTK:
     def _patch_execution_stamp(self, patch_execution_stamp):
         print('patched damask.util.execution_stamp')
 
+    def test_show(sef,default,monkeypatch):
+        monkeypatch.delenv('DISPLAY',raising=False)
+        default.show()
+
     def test_rectilinearGrid(self,tmp_path):
         cells  = np.random.randint(5,10,3)*2
         size   = np.random.random(3) + 1.0

src/grid/spectral_utilities.f90

@@ -156,7 +156,7 @@ subroutine spectral_utilities_init
   integer(C_INTPTR_T) :: alloc_local, local_K, local_K_offset
   integer(C_INTPTR_T), parameter :: &
     scalarSize = 1_C_INTPTR_T, &
-    vecSize    = 3_C_INTPTR_T, &
+    vectorSize = 3_C_INTPTR_T, &
     tensorSize = 9_C_INTPTR_T
   character(len=*), parameter :: &
     PETSCDEBUG = ' -snes_view -snes_monitor '
@@ -274,7 +274,7 @@ subroutine spectral_utilities_init
   call c_f_pointer(tensorField, tensorField_fourier, [3_C_INTPTR_T,3_C_INTPTR_T, &
                    gridFFTW(1)/2_C_INTPTR_T + 1_C_INTPTR_T ,              gridFFTW(2),local_K])     ! place a pointer for a fourier tensor representation
 
-  vectorField = fftw_alloc_complex(vecSize*alloc_local)
+  vectorField = fftw_alloc_complex(vectorSize*alloc_local)
   call c_f_pointer(vectorField, vectorField_real,   [3_C_INTPTR_T,&
                    2_C_INTPTR_T*(gridFFTW(1)/2_C_INTPTR_T + 1_C_INTPTR_T),gridFFTW(2),local_K])     ! place a pointer for a real vector representation
   call c_f_pointer(vectorField, vectorField_fourier,[3_C_INTPTR_T,&
@@ -288,42 +288,42 @@ subroutine spectral_utilities_init
 
 !--------------------------------------------------------------------------------------------------
 ! tensor MPI fftw plans
-  planTensorForth = fftw_mpi_plan_many_dft_r2c(3, [gridFFTW(3),gridFFTW(2),gridFFTW(1)], &          ! dimension, logical length in each dimension in reversed order
+  planTensorForth = fftw_mpi_plan_many_dft_r2c(3,gridFFTW(3:1:-1),tensorSize, &
-                                        tensorSize, FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, &! no. of transforms, default iblock and oblock
+                                               FFTW_MPI_DEFAULT_BLOCK,FFTW_MPI_DEFAULT_BLOCK, &
-                                                       tensorField_real, tensorField_fourier, &     ! input data, output data
+                                               tensorField_real,tensorField_fourier, &
-                                                               PETSC_COMM_WORLD, FFTW_planner_flag) ! use all processors, planer precision
+                                               PETSC_COMM_WORLD,FFTW_planner_flag)
   if (.not. c_associated(planTensorForth)) error stop 'FFTW error'
-  planTensorBack  = fftw_mpi_plan_many_dft_c2r(3, [gridFFTW(3),gridFFTW(2),gridFFTW(1)], &          ! dimension, logical length in each dimension in reversed order
+  planTensorBack  = fftw_mpi_plan_many_dft_c2r(3,gridFFTW(3:1:-1),tensorSize, &
-                                       tensorSize,  FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, &! no. of transforms, default iblock and oblock
+                                               FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, &
-                                                        tensorField_fourier,tensorField_real, &     ! input data, output data
+                                               tensorField_fourier,tensorField_real, &
-                                                               PETSC_COMM_WORLD, FFTW_planner_flag) ! all processors, planer precision
+                                               PETSC_COMM_WORLD, FFTW_planner_flag)
   if (.not. c_associated(planTensorBack)) error stop 'FFTW error'
 
 !--------------------------------------------------------------------------------------------------
 ! vector MPI fftw plans
-  planVectorForth = fftw_mpi_plan_many_dft_r2c(3, [gridFFTW(3),gridFFTW(2),gridFFTW(1)], &          ! dimension, logical length in each dimension in reversed order
+  planVectorForth = fftw_mpi_plan_many_dft_r2c(3,gridFFTW(3:1:-1),vectorSize, &
-                                           vecSize, FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK,&! no. of transforms, default iblock and oblock
+                                               FFTW_MPI_DEFAULT_BLOCK,FFTW_MPI_DEFAULT_BLOCK, &
-                                                       vectorField_real, vectorField_fourier, &     ! input data, output data
+                                               vectorField_real,vectorField_fourier, &
-                                                               PETSC_COMM_WORLD, FFTW_planner_flag) ! use all processors, planer precision
+                                               PETSC_COMM_WORLD,FFTW_planner_flag)
-  if (.not. C_ASSOCIATED(planVectorForth)) error stop 'FFTW error'
+  if (.not. c_associated(planVectorForth)) error stop 'FFTW error'
-  planVectorBack  = fftw_mpi_plan_many_dft_c2r(3, [gridFFTW(3),gridFFTW(2),gridFFTW(1)],  &         ! dimension, logical length in each dimension in reversed order
+  planVectorBack  = fftw_mpi_plan_many_dft_c2r(3,gridFFTW(3:1:-1),vectorSize, &
-                                         vecSize, FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, & ! no. of transforms, default iblock and oblock
+                                               FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, &
-                                                      vectorField_fourier,vectorField_real, &       ! input data, output data
+                                               vectorField_fourier,vectorField_real, &
-                                                               PETSC_COMM_WORLD, FFTW_planner_flag) ! all processors, planer precision
+                                               PETSC_COMM_WORLD, FFTW_planner_flag)
-  if (.not. C_ASSOCIATED(planVectorBack)) error stop 'FFTW error'
+  if (.not. c_associated(planVectorBack)) error stop 'FFTW error'
 
 !--------------------------------------------------------------------------------------------------
 ! scalar MPI fftw plans
-  planScalarForth = fftw_mpi_plan_many_dft_r2c(3, [gridFFTW(3),gridFFTW(2),gridFFTW(1)], &          ! dimension, logical length in each dimension in reversed order
+  planScalarForth = fftw_mpi_plan_many_dft_r2c(3,gridFFTW(3:1:-1),scalarSize, &
-                                       scalarSize, FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, &! no. of transforms, default iblock and oblock
+                                               FFTW_MPI_DEFAULT_BLOCK,FFTW_MPI_DEFAULT_BLOCK, &
-                                                       scalarField_real, scalarField_fourier, &     ! input data, output data
+                                               scalarField_real,scalarField_fourier, &
-                                                               PETSC_COMM_WORLD, FFTW_planner_flag) ! use all processors, planer precision
+                                               PETSC_COMM_WORLD,FFTW_planner_flag)
-  if (.not. C_ASSOCIATED(planScalarForth)) error stop 'FFTW error'
+  if (.not. c_associated(planScalarForth)) error stop 'FFTW error'
-  planScalarBack  = fftw_mpi_plan_many_dft_c2r(3, [gridFFTW(3),gridFFTW(2),gridFFTW(1)], &          ! dimension, logical length in each dimension in reversed order, no. of transforms
+  planScalarBack  = fftw_mpi_plan_many_dft_c2r(3,gridFFTW(3:1:-1),scalarSize, &
-                                       scalarSize, FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, &! no. of transforms, default iblock and oblock
+                                               FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, &
-                                                       scalarField_fourier,scalarField_real, &      ! input data, output data
+                                               scalarField_fourier,scalarField_real, &
-                                                               PETSC_COMM_WORLD, FFTW_planner_flag) ! use all processors, planer precision
+                                               PETSC_COMM_WORLD, FFTW_planner_flag)
-  if (.not. C_ASSOCIATED(planScalarBack)) error stop 'FFTW error'
+  if (.not. c_associated(planScalarBack)) error stop 'FFTW error'
 
 !--------------------------------------------------------------------------------------------------
 ! calculation of discrete angular frequencies, ordered as in FFTW (wrap around)

src/math.f90

@@ -84,38 +84,34 @@ contains
 subroutine math_init
 
   real(pReal), dimension(4) :: randTest
-  integer :: &
+  integer :: randSize
-    randSize, &
+  integer, dimension(:), allocatable :: seed
-    randomSeed                                                                                      !< fixed seeding for pseudo-random number generator, Default 0: use random seed
-  integer, dimension(:), allocatable :: randInit
   class(tNode), pointer :: &
     num_generic
 
+
   print'(/,1x,a)', '<<<+-  math init  -+>>>'; flush(IO_STDOUT)
 
   num_generic => config_numerics%get('generic',defaultVal=emptyDict)
-  randomSeed  = num_generic%get_asInt('random_seed', defaultVal = 0)
 
   call random_seed(size=randSize)
-  allocate(randInit(randSize))
+  allocate(seed(randSize))
-  if (randomSeed > 0) then
+
-    randInit = randomSeed
+  if (num_generic%contains('random_seed')) then
+    seed = num_generic%get_as1dInt('random_seed',requiredSize=randSize)
   else
     call random_seed()
-    call random_seed(get = randInit)
+    call random_seed(get = seed)
-    randInit(2:randSize) = randInit(1)
+  end if
-  endif
 
-  call random_seed(put = randInit)
+  call random_seed(put = seed)
   call random_number(randTest)
 
-  print'(/,a,i2)',                ' size  of random seed:     ', randSize
+  print'(/,a,i2)',              ' size  of random seed:     ', randSize
-  print'(  a,i0)',                ' value of random seed:     ', randInit(1)
+  print*,                       'value of random seed:     ', seed
-  print'(  a,4(/,26x,f17.14),/)', ' start of random sequence: ', randTest
+  print'(  a,4(/,26x,f17.14))', ' start of random sequence: ', randTest
 
-  call random_seed(put = randInit)
+  call selfTest()
-
-  call selfTest
 
 end subroutine math_init
 

src/rotations.f90

@@ -270,7 +270,7 @@ pure elemental subroutine standardize(self)
 
   class(rotation), intent(inout) :: self
 
-  if (self%q(1) < 0.0_pReal) self%q = - self%q
+  if (sign(1.0_pReal,self%q(1)) < 0.0_pReal) self%q = - self%q
 
 end subroutine standardize
 
@@ -450,7 +450,7 @@ pure function qu2om(qu) result(om)
   om(3,2) = 2.0_pReal*(qu(4)*qu(3)+qu(1)*qu(2))
   om(1,3) = 2.0_pReal*(qu(2)*qu(4)+qu(1)*qu(3))
 
-  if (P < 0.0_pReal) om = transpose(om)
+  if (sign(1.0_pReal,P) < 0.0_pReal) om = transpose(om)
 
 end function qu2om
 
@@ -480,7 +480,7 @@ pure function qu2eu(qu) result(eu)
           atan2( 2.0_pReal*chi, q03-q12 ), &
           atan2(( P*qu(1)*qu(3)+qu(2)*qu(4))*chi, (-P*qu(1)*qu(2)+qu(3)*qu(4))*chi )]
   endif degenerated
-  where(eu<0.0_pReal) eu = mod(eu+2.0_pReal*PI,[2.0_pReal*PI,PI,2.0_pReal*PI])
+  where(sign(1.0_pReal,eu)<0.0_pReal) eu = mod(eu+2.0_pReal*PI,[2.0_pReal*PI,PI,2.0_pReal*PI])
 
 end function qu2eu
 
@@ -602,7 +602,7 @@ pure function om2qu(om) result(qu)
           qu = [ (om(2,1) - om(1,2)) /s,(om(1,3) + om(3,1)) / s,(om(2,3) + om(3,2)) / s,0.25_pReal * s]
       endif
   endif
-  if(qu(1)<0._pReal) qu =-1.0_pReal * qu
+  if(sign(1.0_pReal,qu(1))<0.0_pReal) qu =-1.0_pReal * qu
   qu = qu*[1.0_pReal,P,P,P]
 
 end function om2qu
@@ -628,7 +628,7 @@ pure function om2eu(om) result(eu)
     eu = [atan2(om(1,2),om(1,1)), 0.5_pReal*PI*(1.0_pReal-om(3,3)),0.0_pReal ]
   end if
   where(abs(eu) < 1.e-8_pReal) eu = 0.0_pReal
-  where(eu<0.0_pReal)          eu = mod(eu+2.0_pReal*PI,[2.0_pReal*PI,PI,2.0_pReal*PI])
+  where(sign(1.0_pReal,eu)<0.0_pReal) eu = mod(eu+2.0_pReal*PI,[2.0_pReal*PI,PI,2.0_pReal*PI])
 
 end function om2eu
 
@@ -735,7 +735,7 @@ pure function eu2qu(eu) result(qu)
         -P*sPhi*cos(ee(1)-ee(3)), &
         -P*sPhi*sin(ee(1)-ee(3)), &
         -P*cPhi*sin(ee(1)+ee(3))]
-  if(qu(1) < 0.0_pReal) qu = qu * (-1.0_pReal)
+  if(sign(1.0_pReal,qu(1)) < 0.0_pReal) qu = qu * (-1.0_pReal)
 
 end function eu2qu
 
@@ -792,7 +792,7 @@ pure function eu2ax(eu) result(ax)
   else
     ax(1:3) = -P/tau * [ t*cos(delta), t*sin(delta), sin(sigma) ]                                   ! passive axis-angle pair so a minus sign in front
     ax(4) = alpha
-    if (alpha < 0.0_pReal) ax = -ax                                                                 ! ensure alpha is positive
+    if (sign(1.0_pReal,alpha) < 0.0_pReal) ax = -ax                                                 ! ensure alpha is positive
   end if
 
 end function eu2ax