Merge branch 'MoreImprovements' into PGI-support

CMakeLists.txt

@@ -117,9 +117,7 @@ elseif (DAMASK_SOLVER STREQUAL "fem" OR DAMASK_SOLVER STREQUAL "mesh")
 else ()
     message (FATAL_ERROR "Build target (DAMASK_SOLVER) is not defined")
 endif ()
-
+list(APPEND CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/cmake)
-# set linker commands (needs to be done after defining the project)
-set (CMAKE_LINKER           "${PETSC_LINKER}")
 
 if (CMAKE_BUILD_TYPE STREQUAL "")
   set (CMAKE_BUILD_TYPE "RELEASE")
@@ -168,9 +166,6 @@ add_definitions (-DDAMASKVERSION="${DAMASK_V}")
 # definition of other macros
 add_definitions (-DPETSc)
 
-set (DAMASK_INCLUDE_FLAGS "${DAMASK_INCLUDE_FLAGS} ${PETSC_INCLUDES}")
-list(APPEND CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/cmake)
-
 if (CMAKE_Fortran_COMPILER_ID STREQUAL "Intel")
   include(Compiler-Intel)
 elseif(CMAKE_Fortran_COMPILER_ID STREQUAL "GNU")
@@ -183,14 +178,14 @@ endif ()
 
 
 set (CMAKE_Fortran_FLAGS_${CMAKE_BUILD_TYPE} "${BUILDCMD_PRE} ${OPENMP_FLAGS} ${STANDARD_CHECK} ${OPTIMIZATION_FLAGS} ${COMPILE_FLAGS} ${PRECISION_FLAGS}")
-set (CMAKE_Fortran_LINK_EXECUTABLE           "${BUILDCMD_PRE} ${CMAKE_LINKER} ${OPENMP_FLAGS}   ${OPTIMIZATION_FLAGS} ${LINKER_FLAGS}")
+set (CMAKE_Fortran_LINK_EXECUTABLE           "${BUILDCMD_PRE} ${PETSC_LINKER} ${OPENMP_FLAGS}   ${OPTIMIZATION_FLAGS} ${LINKER_FLAGS}")
 
 if (CMAKE_BUILD_TYPE STREQUAL "DEBUG")
   set (CMAKE_Fortran_FLAGS_${CMAKE_BUILD_TYPE} "${CMAKE_Fortran_FLAGS_${CMAKE_BUILD_TYPE}} ${DEBUG_FLAGS}")
   set (CMAKE_Fortran_LINK_EXECUTABLE           "${CMAKE_Fortran_LINK_EXECUTABLE} ${DEBUG_FLAGS}")
 endif ()
 
-set (CMAKE_Fortran_FLAGS_${CMAKE_BUILD_TYPE}   "${CMAKE_Fortran_FLAGS_${CMAKE_BUILD_TYPE}} ${DAMASK_INCLUDE_FLAGS} ${BUILDCMD_POST}")
+set (CMAKE_Fortran_FLAGS_${CMAKE_BUILD_TYPE}   "${CMAKE_Fortran_FLAGS_${CMAKE_BUILD_TYPE}} ${PETSC_INCLUDES} ${BUILDCMD_POST}")
 set (CMAKE_Fortran_LINK_EXECUTABLE             "${CMAKE_Fortran_LINK_EXECUTABLE} <OBJECTS> -o <TARGET> <LINK_LIBRARIES> ${PETSC_EXTERNAL_LIB} ${BUILDCMD_POST}")
 
 message ("Fortran Compiler Flags:\n${CMAKE_Fortran_FLAGS_${CMAKE_BUILD_TYPE}}\n")

PRIVATE

@@ -1 +1 @@
-Subproject commit cda69f9a59fe64223439a2c725e1a78cf22b28aa
+Subproject commit 99b076706a186ec7deb051ae181c2d9697c799fc

VERSION

@@ -1 +1 @@
-v2.0.3-1294-g034367fa
+v2.0.3-1354-gef5a8a3a

installation/patch/README.md

@@ -9,14 +9,6 @@ cd DAMASK_ROOT
 patch -p1 < installation/patch/nameOfPatch
 ```
 
-## Available patches
-
-  * **disable_HDF5** disables all HDF5 output.
-    HDF5 output is an experimental feature. Also, some routines not present in HDF5 1.8.x are removed to allow compilation of DAMASK with HDF5 < 1.10.x
-
-  * **disable_old_output** disables all non-HDF5 output.
-    Saves some memory when using only HDF5 output
-
 ## Create patch
 commit your changes
 

processing/post/addCompatibilityMismatch.py

@@ -182,8 +182,8 @@ for name in filenames:
   nodes = damask.grid_filters.node_coord(size,F)
   
   if options.shape:
-    centres = damask.grid_filters.cell_coord(size,F)
+    centers = damask.grid_filters.cell_coord(size,F)
-    shapeMismatch = shapeMismatch(  size,table.get(options.defgrad).reshape(grid[2],grid[1],grid[0],3,3),nodes,centres)
+    shapeMismatch = shapeMismatch(  size,table.get(options.defgrad).reshape(grid[2],grid[1],grid[0],3,3),nodes,centers)
     table.add('shapeMismatch(({}))'.format(options.defgrad),
               shapeMismatch.reshape((-1,1)),
               scriptID+' '+' '.join(sys.argv[1:]))

processing/post/addStrainTensors.py

@@ -2,10 +2,9 @@
 
 import os
 import sys
+from io import StringIO
 from optparse import OptionParser
 
-import numpy as np
-
 import damask
 
 

processing/post/addTable.py

@@ -2,6 +2,7 @@
 
 import os
 import sys
+from io import StringIO
 from optparse import OptionParser
 
 import damask

processing/post/averageDown.py

@@ -2,6 +2,7 @@
 
 import os
 import sys
+from io import StringIO
 from optparse import OptionParser
 
 import numpy as np

processing/post/blowUp.py

@@ -2,6 +2,7 @@
 
 import os
 import sys
+from io import StringIO
 from optparse import OptionParser
 
 from scipy import ndimage

processing/post/growTable.py

@@ -2,10 +2,9 @@
 
 import os
 import sys
+from io import StringIO
 from optparse import OptionParser
 
-import numpy as np
-
 import damask
 
 scriptName = os.path.splitext(os.path.basename(__file__))[0]

python/damask/colormaps.py

@@ -328,9 +328,9 @@ class Color():
     Msh = np.zeros(3,'d')
     Msh[0] = np.sqrt(np.dot(self.color,self.color))
     if (Msh[0] > 0.001):
-      Msh[1] = np.acos(self.color[0]/Msh[0])
+      Msh[1] = np.arccos(self.color[0]/Msh[0])
       if (self.color[1] != 0.0):
-        Msh[2] = np.atan2(self.color[2],self.color[1])
+        Msh[2] = np.arctan2(self.color[2],self.color[1])
 
     converted = Color('MSH', Msh)
     self.model = converted.model

python/damask/dadf5.py

@@ -880,7 +880,7 @@ class DADF5():
       else:
         
         nodes = vtk.vtkPoints()
-        with h5py.File(self.fname) as f:
+        with h5py.File(self.fname,'r') as f:
           nodes.SetData(numpy_support.numpy_to_vtk(f['/geometry/x_n'][()],deep=True))
           
           vtk_geom = vtk.vtkUnstructuredGrid()

python/damask/orientation.py

@@ -170,9 +170,18 @@ class Rotation:
     ################################################################################################
     # convert to different orientation representations (numpy arrays)
 
-    def asQuaternion(self):
+    def asQuaternion(self,
-      """Unit quaternion: (q, p_1, p_2, p_3)."""
+                     quaternion = False):
-      return self.quaternion.asArray()
+      """
+      Unit quaternion [q, p_1, p_2, p_3] unless quaternion == True: damask.quaternion object.
+
+      Parameters
+      ----------
+      quaternion : bool, optional
+          return quaternion as DAMASK object.
+      
+      """
+      return self.quaternion if quaternion else self.quaternion.asArray()
       
     def asEulers(self,
                  degrees = False):
@@ -190,33 +199,36 @@ class Rotation:
       return eu
     
     def asAxisAngle(self,
-                    degrees = False):
+                    degrees = False,
+                    pair = False):
       """
-      Axis angle pair: ([n_1, n_2, n_3], ω).
+      Axis angle representation [n_1, n_2, n_3, ω] unless pair == True: ([n_1, n_2, n_3], ω).
 
       Parameters
       ----------
       degrees : bool, optional
           return rotation angle in degrees.
+      pair : bool, optional
+          return tuple of axis and angle.
       
       """
       ax = qu2ax(self.quaternion.asArray())
       if degrees: ax[3] = np.degrees(ax[3])
-      return ax
+      return (ax[:3],np.degrees(ax[3])) if pair else ax
       
     def asMatrix(self):
       """Rotation matrix."""
       return qu2om(self.quaternion.asArray())
 
     def asRodrigues(self,
-                    vector=False):
+                    vector = False):
       """
-      Rodrigues-Frank vector: ([n_1, n_2, n_3], tan(ω/2)).
+      Rodrigues-Frank vector representation [n_1, n_2, n_3, tan(ω/2)] unless vector == True: [n_1, n_2, n_3] * tan(ω/2).
 
       Parameters
       ----------
       vector : bool, optional
-          return as array of length 3, i.e. scale the unit vector giving the rotation axis.
+          return as actual Rodrigues--Frank vector, i.e. rotation axis scaled by tan(ω/2).
 
       """
       ro = qu2ro(self.quaternion.asArray())
@@ -252,8 +264,8 @@ class Rotation:
                        acceptHomomorph = False,
                        P = -1):
 
-      qu = quaternion if isinstance(quaternion, np.ndarray) and quaternion.dtype == np.dtype(float) \
+      qu =   quaternion if isinstance(quaternion, np.ndarray) and quaternion.dtype == np.dtype(float) \
-                      else np.array(quaternion,dtype=float)
+                        else np.array(quaternion,dtype=float)
       if P > 0: qu[1:4] *= -1                                                                       # convert from P=1 to P=-1
       if qu[0] < 0.0:
         if acceptHomomorph:
@@ -1193,9 +1205,9 @@ class Orientation:
     ref = orientations[0]
     for o in orientations:
       closest.append(o.equivalentOrientations(
-                    ref.disorientation(o,
+                     ref.disorientation(o,
-                                       SST = False,                                                 # select (o[ther]'s) sym orientation
+                                        SST = False,                                                # select (o[ther]'s) sym orientation
-                                       symmetries = True)[2]).rotation)                             # with lowest misorientation
+                                        symmetries = True)[2]).rotation)                            # with lowest misorientation
 
     return Orientation(Rotation.fromAverage(closest,weights),ref.lattice)
 

python/damask/table.py

@@ -22,7 +22,7 @@ class Table():
             Additional, human-readable information.
         
         """
-        self.comments = ['table.py v {}'.format(version)] if not comments else [c for c in comments]
+        self.comments = [] if comments is None else [c for c in comments]
         self.data = pd.DataFrame(data=data)
         self.shapes = shapes
         self.__label_condensed()
@@ -77,10 +77,9 @@ class Table():
         if keyword == 'header':
             header = int(header)
         else:
-            raise Exception
+            raise TypeError
 
-        comments = ['table.py:from_ASCII v {}'.format(version)]
+        comments = [f.readline()[:-1] for i in range(1,header)]
-        comments+= [f.readline()[:-1] for i in range(1,header)]
         labels   = f.readline().split()
        
         shapes = {}
@@ -138,9 +137,12 @@ class Table():
                 break
        
         data = np.loadtxt(content)
+        for c in range(data.shape[1]-10):
+            shapes['n/a_{}'.format(c+1)] = (1,)
 
         return Table(data,shapes,comments)
 
+
     @property
     def labels(self):
         return list(self.shapes.keys())
@@ -271,7 +273,9 @@ class Table():
 
     def append(self,other):
         """
-        Append other table vertically (similar to numpy.vstack). Requires matching shapes and order.
+        Append other table vertically (similar to numpy.vstack).
+
+        Requires matching labels/shapes and order.
 
         Parameters
         ----------
@@ -287,8 +291,9 @@ class Table():
 
     def join(self,other):
         """
-        Append other table horizontally (similar to numpy.hstack). Requires matching number of rows
+        Append other table horizontally (similar to numpy.hstack).
-        and no common lables
+
+        Requires matching number of rows and no common labels.
 
         Parameters
         ----------

python/tests/test_mechanics.py

@@ -113,7 +113,7 @@ class TestMechanics:
    
     def test_strain_tensor_rotation_equivalence(self):
          """Ensure that left and right strain differ only by a rotation."""
-         F = np.random.random((self.n,3,3))
+         F = np.broadcast_to(np.eye(3),[self.n,3,3]) + (np.random.random((self.n,3,3))*0.5 - 0.25)
          m = np.random.random()*5.0-2.5
          assert np.allclose(np.linalg.det(mechanics.strain_tensor(F,'U',m)),
                             np.linalg.det(mechanics.strain_tensor(F,'V',m)))

src/CPFEM.f90

@@ -377,7 +377,7 @@ subroutine CPFEM_results(inc,time)
   call constitutive_results
   call crystallite_results
   call homogenization_results
-  call results_removeLink('current') ! ToDo: put this into closeJobFile
+  call results_finalizeIncrement
   call results_closeJobFile
 
 end subroutine CPFEM_results

src/CPFEM2.f90

@@ -201,7 +201,7 @@ subroutine CPFEM_results(inc,time)
   call crystallite_results
   call homogenization_results
   call discretization_results
-  call results_removeLink('current') ! ToDo: put this into closeJobFile?
+  call results_finalizeIncrement
   call results_closeJobFile
 
 end subroutine CPFEM_results

src/DAMASK_interface.f90

@@ -269,10 +269,10 @@ subroutine DAMASK_interface_init
     write(6,'(a,a)')    ' Working dir argument:   ', trim(workingDirArg)
   write(6,'(a,a)')      ' Geometry argument:      ', trim(geometryArg)
   write(6,'(a,a)')      ' Load case argument:     ', trim(loadcaseArg)
-  write(6,'(a,a)')      ' Working directory:      ', trim(getCWD())
+  write(6,'(a,a)')      ' Working directory:      ', getCWD()
   write(6,'(a,a)')      ' Geometry file:          ', trim(geometryFile)
   write(6,'(a,a)')      ' Loadcase file:          ', trim(loadCaseFile)
-  write(6,'(a,a)')      ' Solver job name:        ', trim(getSolverJobName())
+  write(6,'(a,a)')      ' Solver job name:        ', getSolverJobName()
   if (interface_restartInc > 0) &
     write(6,'(a,i6.6)') ' Restart from increment: ', interface_restartInc
 
@@ -308,7 +308,7 @@ subroutine setWorkingDirectory(workingDirectoryArg)
   workingDirectory = trim(rectifyPath(workingDirectory))
   error = setCWD(trim(workingDirectory))
   if(error) then
-    write(6,'(/,a)') ' ERROR: Working directory "'//trim(workingDirectory)//'" does not exist'
+    write(6,'(/,a)') ' ERROR: Invalid Working directory: '//trim(workingDirectory)
     call quit(1)
   endif
 
@@ -318,8 +318,9 @@ end subroutine setWorkingDirectory
 !--------------------------------------------------------------------------------------------------
 !> @brief solver job name (no extension) as combination of geometry and load case name
 !--------------------------------------------------------------------------------------------------
-character(len=1024) function getSolverJobName()
+function getSolverJobName()
 
+  character(len=:), allocatable :: getSolverJobName
   integer :: posExt,posSep
 
   posExt = scan(geometryFile,'.',back=.true.)
@@ -330,7 +331,7 @@ character(len=1024) function getSolverJobName()
   posExt = scan(loadCaseFile,'.',back=.true.)
   posSep = scan(loadCaseFile,'/',back=.true.)
 
-  getSolverJobName = trim(getSolverJobName)//'_'//loadCaseFile(posSep+1:posExt-1)
+  getSolverJobName = getSolverJobName//'_'//loadCaseFile(posSep+1:posExt-1)
 
 end function getSolverJobName
 
@@ -338,15 +339,16 @@ end function getSolverJobName
 !--------------------------------------------------------------------------------------------------
 !> @brief basename of geometry file with extension from command line arguments
 !--------------------------------------------------------------------------------------------------
-character(len=1024) function getGeometryFile(geometryParameter)
+function getGeometryFile(geometryParameter)
 
-  character(len=1024), intent(in) :: geometryParameter
+  character(len=:), allocatable :: getGeometryFile
-  logical                         :: file_exists
+  character(len=*),  intent(in) :: geometryParameter
-  external                        :: quit
+  logical                       :: file_exists
+  external                      :: quit
 
   getGeometryFile = trim(geometryParameter)
-  if (scan(getGeometryFile,'/') /= 1) getGeometryFile = trim(getCWD())//'/'//trim(getGeometryFile)
+  if (scan(getGeometryFile,'/') /= 1) getGeometryFile = getCWD()//'/'//trim(getGeometryFile)
-  getGeometryFile = makeRelativePath(trim(getCWD()), getGeometryFile)
+  getGeometryFile = makeRelativePath(getCWD(), getGeometryFile)
 
   inquire(file=trim(getGeometryFile), exist=file_exists)
   if (.not. file_exists) then
@@ -360,15 +362,16 @@ end function getGeometryFile
 !--------------------------------------------------------------------------------------------------
 !> @brief relative path of loadcase from command line arguments
 !--------------------------------------------------------------------------------------------------
-character(len=1024) function getLoadCaseFile(loadCaseParameter)
+function getLoadCaseFile(loadCaseParameter)
 
-  character(len=1024), intent(in) :: loadCaseParameter
+  character(len=:), allocatable :: getLoadCaseFile
-  logical                         :: file_exists
+  character(len=*),  intent(in) :: loadCaseParameter
-  external                        :: quit
+  logical                       :: file_exists
+  external                      :: quit
 
   getLoadCaseFile = trim(loadCaseParameter)
-  if (scan(getLoadCaseFile,'/') /= 1) getLoadCaseFile = trim(getCWD())//'/'//trim(getLoadCaseFile)
+  if (scan(getLoadCaseFile,'/') /= 1) getLoadCaseFile = getCWD()//'/'//trim(getLoadCaseFile)
-  getLoadCaseFile = makeRelativePath(trim(getCWD()), getLoadCaseFile)
+  getLoadCaseFile = makeRelativePath(getCWD(), getLoadCaseFile)
 
   inquire(file=trim(getLoadCaseFile), exist=file_exists)
   if (.not. file_exists) then

src/FEsolving.f90

@@ -46,10 +46,10 @@ subroutine FE_init
     call IO_open_inputFile(FILEUNIT)
     rewind(FILEUNIT)
     do
-      read (FILEUNIT,'(a256)',END=100) line
+      read (FILEUNIT,'(A)',END=100) line
       chunkPos = IO_stringPos(line)
       if(IO_lc(IO_stringValue(line,chunkPos,1)) == 'solver') then
-        read (FILEUNIT,'(a256)',END=100) line                                                       ! next line
+        read (FILEUNIT,'(A)',END=100) line                                                          ! next line
         chunkPos = IO_stringPos(line)
         symmetricSolver = (IO_intValue(line,chunkPos,2) /= 1)
       endif

src/IO.f90

@@ -23,7 +23,7 @@ module IO
   public :: &
     IO_init, &
     IO_read_ASCII, &
-    IO_open_file, &
+    IO_open_file, &                                                                                 ! deprecated, use IO_read_ASCII
     IO_open_jobFile_binary, &
     IO_isBlank, &
     IO_getTag, &
@@ -44,8 +44,7 @@ module IO
     IO_countDataLines
 #elif defined(Marc4DAMASK)
     IO_fixedNoEFloatValue, &
-    IO_fixedIntValue, &
+    IO_fixedIntValue
-    IO_countNumericalDataLines
 #endif
 #endif
 
@@ -245,7 +244,7 @@ subroutine IO_open_inputFile(fileUnit)
    
    
     do
-      read(unit2,'(A256)',END=220) line
+      read(unit2,'(A)',END=220) line
       chunkPos = IO_stringPos(line)
    
       if (IO_lc(IO_StringValue(line,chunkPos,1))=='*include') then
@@ -384,7 +383,7 @@ function IO_stringValue(string,chunkPos,myChunk,silent)
   logical                                            :: warn
 
   if (present(silent)) then
-    warn = silent
+    warn = .not. silent
   else
     warn = .false.
   endif
@@ -414,11 +413,10 @@ real(pReal) function IO_floatValue (string,chunkPos,myChunk)
 
   valuePresent: if (myChunk > chunkPos(1) .or. myChunk < 1) then
     call IO_warning(201,el=myChunk,ext_msg=MYNAME//trim(string))
-  else  valuePresent
+  else valuePresent
-    IO_floatValue = &
+    IO_floatValue = verifyFloatValue(trim(adjustl(string(chunkPos(myChunk*2):chunkPos(myChunk*2+1)))),&
-                verifyFloatValue(trim(adjustl(string(chunkPos(myChunk*2):chunkPos(myChunk*2+1)))),&
+                                     VALIDCHARACTERS,MYNAME)
-                                        VALIDCHARACTERS,MYNAME)
+  endif valuePresent
-  endif  valuePresent
 
 end function IO_floatValue
 
@@ -466,12 +464,12 @@ real(pReal) function IO_fixedNoEFloatValue (string,ends,myChunk)
   pos_exp = scan(string(ends(myChunk)+1:ends(myChunk+1)),'+-',back=.true.)
   hasExponent: if (pos_exp > 1) then
     base  = verifyFloatValue(trim(adjustl(string(ends(myChunk)+1:ends(myChunk)+pos_exp-1))),&
-                                VALIDBASE,MYNAME//'(base): ')
+                             VALIDBASE,MYNAME//'(base): ')
     expon = verifyIntValue(trim(adjustl(string(ends(myChunk)+pos_exp:ends(myChunk+1)))),&
-                                VALIDEXP,MYNAME//'(exp): ')
+                           VALIDEXP,MYNAME//'(exp): ')
   else hasExponent
     base  = verifyFloatValue(trim(adjustl(string(ends(myChunk)+1:ends(myChunk+1)))),&
-                                VALIDBASE,MYNAME//'(base): ')
+                             VALIDBASE,MYNAME//'(base): ')
     expon = 0
   endif hasExponent
   IO_fixedNoEFloatValue = base*10.0_pReal**real(expon,pReal)
@@ -530,8 +528,8 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg)
   character(len=*), optional, intent(in) :: ext_msg
 
   external                               :: quit
-  character(len=1024)                    :: msg
+  character(len=pStringLen)              :: msg
-  character(len=1024)                    :: formatString
+  character(len=pStringLen)              :: formatString
 
   select case (error_ID)
 
@@ -550,14 +548,8 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg)
       msg = 'could not read file:'
     case (103)
       msg = 'could not assemble input files'
-    case (104)
-      msg = '{input} recursion limit reached'
-    case (105)
-      msg = 'unknown output:'
     case (106)
       msg = 'working directory does not exist:'
-    case (107)
-      msg = 'line length exceeds limit of 256'
 
 !--------------------------------------------------------------------------------------------------
 ! lattice error messages
@@ -777,8 +769,8 @@ subroutine IO_warning(warning_ID,el,ip,g,ext_msg)
   integer,          optional, intent(in) :: el,ip,g
   character(len=*), optional, intent(in) :: ext_msg
  
-  character(len=1024)                    :: msg
+  character(len=pStringLen)              :: msg
-  character(len=1024)                    :: formatString
+  character(len=pStringLen)              :: formatString
  
   select case (warning_ID)
     case (1)
@@ -924,38 +916,6 @@ end function IO_countDataLines
 #endif
 
 
-#ifdef Marc4DAMASK
-!--------------------------------------------------------------------------------------------------
-!> @brief count lines containig data up to next *keyword
-!--------------------------------------------------------------------------------------------------
-integer function IO_countNumericalDataLines(fileUnit)
-
-  integer, intent(in)                :: fileUnit                                                    !< file handle
- 
- 
-  integer, allocatable, dimension(:) :: chunkPos
-  character(len=pStringLen)          :: line, &
-                                        tmp
- 
-  IO_countNumericalDataLines = 0
-  line = ''
- 
-  do while (trim(line) /= IO_EOF)
-    line = IO_read(fileUnit)
-    chunkPos = IO_stringPos(line)
-    tmp = IO_lc(IO_stringValue(line,chunkPos,1))
-    if (verify(trim(tmp),'0123456789') == 0) then                                                   ! numerical values
-      IO_countNumericalDataLines = IO_countNumericalDataLines + 1
-    else
-      exit
-    endif
-  enddo
-  backspace(fileUnit)
-
-end function IO_countNumericalDataLines
-#endif
-
-
 !--------------------------------------------------------------------------------------------------
 !> @brief count items in consecutive lines depending on lines
 !> @details Marc:      ints concatenated by "c" as last char or range of values a "to" b
@@ -1040,7 +1000,7 @@ function IO_continuousIntValues(fileUnit,maxN,lookupName,lookupMap,lookupMaxN)
 
 #if defined(Marc4DAMASK)
  do
-   read(fileUnit,'(A256)',end=100) line
+   read(fileUnit,'(A)',end=100) line
    chunkPos = IO_stringPos(line)
    if (chunkPos(1) < 1) then                                                                        ! empty line
      exit
@@ -1081,14 +1041,14 @@ function IO_continuousIntValues(fileUnit,maxN,lookupName,lookupMap,lookupMaxN)
 !--------------------------------------------------------------------------------------------------
 ! check if the element values in the elset are auto generated
  backspace(fileUnit)
- read(fileUnit,'(A256)',end=100) line
+ read(fileUnit,'(A)',end=100) line
  chunkPos = IO_stringPos(line)
  do i = 1,chunkPos(1)
    if (IO_lc(IO_stringValue(line,chunkPos,i)) == 'generate') rangeGeneration = .true.
  enddo
 
  do l = 1,c
-   read(fileUnit,'(A256)',end=100) line
+   read(fileUnit,'(A)',end=100) line
    chunkPos = IO_stringPos(line)
    if (verify(IO_stringValue(line,chunkPos,1),'0123456789') > 0) then                               ! a non-int, i.e. set names follow on this line
      do i = 1,chunkPos(1)                                                                           ! loop over set names in line

src/element.f90

@@ -39,7 +39,7 @@ module element
   integer, parameter, private :: &
     NELEMTYPE = 13
  
-  integer, dimension(NelemType), parameter, private :: NNODE = &
+  integer, dimension(NELEMTYPE), parameter, private :: NNODE = &
     [ &
         3, & ! 2D 3node 1ip
         6, & ! 2D 6node 3ip
@@ -57,7 +57,7 @@ module element
        20  & ! 3D 20node 27ip
     ]                                                                                               !< number of nodes that constitute a specific type of element
  
-  integer, dimension(NelemType), parameter, public :: GEOMTYPE = &
+  integer, dimension(NELEMTYPE), parameter, public :: GEOMTYPE = &
     [ &
         1, & 
         2, & 
@@ -74,8 +74,7 @@ module element
        10  & 
     ]                                                                                               !< geometry type of particular element type
  
-  !integer, dimension(maxval(geomType)), parameter, private :: NCELLNODE = &                        ! Intel 16.0 complains
+  integer, dimension(maxval(GEOMTYPE)), parameter, private :: NCELLNODE = &
-  integer, dimension(10), parameter, private :: NCELLNODE = &
     [ &
         3, &
         7, &
@@ -89,8 +88,7 @@ module element
        64  &
     ]                                                                                               !< number of cell nodes in a specific geometry type
  
-  !integer, dimension(maxval(geomType)), parameter, private :: NIP = &                              ! Intel 16.0 complains
+  integer, dimension(maxval(GEOMTYPE)), parameter, private :: NIP = &
-  integer, dimension(10), parameter, private :: NIP = &
     [ &
         1, & 
         3, & 
@@ -104,8 +102,7 @@ module element
        27  & 
     ]                                                                                               !< number of IPs in a specific geometry type
  
-  !integer, dimension(maxval(geomType)), parameter, private  :: CELLTYPE = &                        ! Intel 16.0 complains
+  integer, dimension(maxval(GEOMTYPE)), parameter, private  :: CELLTYPE = &
-  integer, dimension(10), parameter, private  :: CELLTYPE = &
     [ &
        1, & ! 2D 3node
        2, & ! 2D 4node
@@ -119,8 +116,7 @@ module element
        4  & ! 3D 8node
     ]                                                                                               !< cell type that is used by each geometry type
 
-  !integer, dimension(maxval(cellType)), parameter, private :: nIPNeighbor = &                      ! Intel 16.0 complains
+  integer, dimension(maxval(CELLTYPE)), parameter, private :: NIPNEIGHBOR = &
-  integer, dimension(4), parameter, private :: NIPNEIGHBOR = &
     [ &
        3, & ! 2D 3node
        4, & ! 2D 4node
@@ -128,8 +124,7 @@ module element
        6  & ! 3D 8node
     ]                                                                                               !< number of ip neighbors / cell faces in a specific cell type 
 
-  !integer, dimension(maxval(cellType)), parameter, private :: NCELLNODESPERCELLFACE = &            ! Intel 16.0 complains
+  integer, dimension(maxval(CELLTYPE)), parameter, private :: NCELLNODEPERCELLFACE = &
-  integer, dimension(4), parameter, private :: NCELLNODEPERCELLFACE = &
     [ &
        2, & ! 2D 3node
        2, & ! 2D 4node
@@ -137,8 +132,7 @@ module element
        4  & ! 3D 8node
     ]                                                                                               !< number of cell nodes in a specific cell type 
 
-  !integer, dimension(maxval(CELLTYPE)), parameter, private  :: NCELLNODEPERCELL = &                ! Intel 16.0 complains
+  integer, dimension(maxval(CELLTYPE)), parameter, private  :: NCELLNODEPERCELL = &
-  integer, dimension(4), parameter, private :: NCELLNODEPERCELL = &
     [ &
        3, & ! 2D 3node
        4, & ! 2D 4node

src/grid/grid_damage_spectral.f90

@@ -54,7 +54,7 @@ contains
 !--------------------------------------------------------------------------------------------------
 subroutine grid_damage_spectral_init
     
-  PetscInt, dimension(worldsize) :: localK  
+  PetscInt, dimension(0:worldsize-1) :: localK  
   integer :: i, j, k, cell
   DM :: damage_grid
   Vec :: uBound, lBound
@@ -78,8 +78,8 @@ subroutine grid_damage_spectral_init
 ! initialize solver specific parts of PETSc
   call SNESCreate(PETSC_COMM_WORLD,damage_snes,ierr); CHKERRQ(ierr)
   call SNESSetOptionsPrefix(damage_snes,'damage_',ierr);CHKERRQ(ierr) 
-  localK              = 0
+  localK            = 0
-  localK(worldrank+1) = grid3
+  localK(worldrank) = grid3
   call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,PETSC_COMM_WORLD,ierr)
   call DMDACreate3D(PETSC_COMM_WORLD, &
          DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, &                                    ! cut off stencil at boundary

src/grid/grid_mech_FEM.f90

@@ -53,7 +53,7 @@ module grid_mech_FEM
     F_aim_lastInc  = math_I3, &                                                                     !< previous average deformation gradient
     P_av = 0.0_pReal                                                                                !< average 1st Piola--Kirchhoff stress
  
-  character(len=1024), private :: incInfo                                                           !< time and increment information
+  character(len=pStringLen), private :: incInfo                                                     !< time and increment information
  
   real(pReal), private, dimension(3,3,3,3) :: &
     C_volAvg = 0.0_pReal, &                                                                         !< current volume average stiffness 
@@ -81,7 +81,7 @@ contains
 subroutine grid_mech_FEM_init
     
   real(pReal) :: HGCoeff = 0.0e-2_pReal
-  PetscInt, dimension(worldsize) :: localK
+  PetscInt, dimension(0:worldsize-1) :: localK
   real(pReal), dimension(3,3) :: &
     temp33_Real = 0.0_pReal
   real(pReal), dimension(4,8) :: &
@@ -94,7 +94,6 @@ subroutine grid_mech_FEM_init
                        1.0_pReal,-1.0_pReal,-1.0_pReal,-1.0_pReal, &
                        1.0_pReal, 1.0_pReal, 1.0_pReal, 1.0_pReal], [4,8])
   PetscErrorCode :: ierr
-  integer        :: rank
   integer(HID_T) :: fileHandle, groupHandle
   character(len=pStringLen) :: fileName
   real(pReal), dimension(3,3,3,3) :: devNull
@@ -121,8 +120,8 @@ subroutine grid_mech_FEM_init
 ! initialize solver specific parts of PETSc
   call SNESCreate(PETSC_COMM_WORLD,mech_snes,ierr); CHKERRQ(ierr)
   call SNESSetOptionsPrefix(mech_snes,'mech_',ierr);CHKERRQ(ierr) 
-  localK              = 0
+  localK            = 0
-  localK(worldrank+1) = grid3
+  localK(worldrank) = grid3
   call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,PETSC_COMM_WORLD,ierr)
   call DMDACreate3d(PETSC_COMM_WORLD, &
          DM_BOUNDARY_PERIODIC, DM_BOUNDARY_PERIODIC, DM_BOUNDARY_PERIODIC, &

src/grid/grid_mech_spectral_basic.f90

@@ -55,7 +55,7 @@ module grid_mech_spectral_basic
     F_aim_lastInc = math_I3, &                                                                      !< previous average deformation gradient
     P_av = 0.0_pReal                                                                                !< average 1st Piola--Kirchhoff stress
 
-  character(len=1024), private :: incInfo                                                           !< time and increment information 
+  character(len=pStringLen), private :: incInfo                                                     !< time and increment information 
   real(pReal), private, dimension(3,3,3,3) :: &
     C_volAvg = 0.0_pReal, &                                                                         !< current volume average stiffness 
     C_volAvgLastInc = 0.0_pReal, &                                                                  !< previous volume average stiffness

src/grid/grid_mech_spectral_polarisation.f90

@@ -59,7 +59,7 @@ module grid_mech_spectral_polarisation
     F_av = 0.0_pReal, &                                                                             !< average incompatible def grad field
     P_av = 0.0_pReal                                                                                !< average 1st Piola--Kirchhoff stress
   
-  character(len=1024), private :: incInfo                                                           !< time and increment information
+  character(len=pStringLen), private :: incInfo                                                     !< time and increment information
   real(pReal), private, dimension(3,3,3,3) :: &
     C_volAvg = 0.0_pReal, &                                                                         !< current volume average stiffness 
     C_volAvgLastInc = 0.0_pReal, &                                                                  !< previous volume average stiffness
@@ -100,7 +100,7 @@ subroutine grid_mech_spectral_polarisation_init
     FandF_tau, &                                                                                    ! overall pointer to solution data
     F, &                                                                                            ! specific (sub)pointer
     F_tau                                                                                           ! specific (sub)pointer
-  PetscInt, dimension(worldsize) :: localK 
+  PetscInt, dimension(0:worldsize-1) :: localK 
   integer(HID_T) :: fileHandle, groupHandle
   integer        :: fileUnit
   character(len=pStringLen) :: fileName
@@ -130,8 +130,8 @@ subroutine grid_mech_spectral_polarisation_init
 ! initialize solver specific parts of PETSc
   call SNESCreate(PETSC_COMM_WORLD,snes,ierr); CHKERRQ(ierr)
   call SNESSetOptionsPrefix(snes,'mech_',ierr);CHKERRQ(ierr) 
-  localK              = 0
+  localK            = 0
-  localK(worldrank+1) = grid3
+  localK(worldrank) = grid3
   call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,PETSC_COMM_WORLD,ierr)
   call DMDACreate3d(PETSC_COMM_WORLD, &
          DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, &                                    ! cut off stencil at boundary

src/grid/grid_thermal_spectral.f90

@@ -55,7 +55,7 @@ contains
 !--------------------------------------------------------------------------------------------------
 subroutine grid_thermal_spectral_init
 
-  PetscInt, dimension(worldsize) :: localK  
+  PetscInt, dimension(0:worldsize-1) :: localK  
   integer :: i, j, k, cell
   DM :: thermal_grid
   PetscScalar,  dimension(:,:,:), pointer :: x_scal
@@ -77,8 +77,8 @@ subroutine grid_thermal_spectral_init
 ! initialize solver specific parts of PETSc
   call SNESCreate(PETSC_COMM_WORLD,thermal_snes,ierr); CHKERRQ(ierr)
   call SNESSetOptionsPrefix(thermal_snes,'thermal_',ierr);CHKERRQ(ierr) 
-  localK              = 0
+  localK            = 0
-  localK(worldrank+1) = grid3
+  localK(worldrank) = grid3
   call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,PETSC_COMM_WORLD,ierr)
   call DMDACreate3D(PETSC_COMM_WORLD, &
          DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, &                                    ! cut off stencil at boundary

src/grid/spectral_utilities.f90

@@ -700,7 +700,7 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C)
     c_reduced, &                                                                                    !< reduced stiffness (depending on number of stress BC)
     sTimesC                                                                                         !< temp variable to check inversion
   logical :: errmatinv
-  character(len=1024):: formatString
+  character(len=pStringLen):: formatString
  
   mask_stressVector = reshape(transpose(mask_stress), [9])
   size_reduced = count(mask_stressVector)

src/kinematics_cleavage_opening.f90

@@ -65,7 +65,7 @@ subroutine kinematics_cleavage_opening_init
 
   integer :: maxNinstance,p,instance
 
-  write(6,'(/,a)')   ' <<<+-  kinematics_'//KINEMATICS_cleavage_opening_LABEL//' init  -+>>>'
+  write(6,'(/,a)') ' <<<+-  kinematics_'//KINEMATICS_cleavage_opening_LABEL//' init  -+>>>'; flush(6)
 
   maxNinstance = count(phase_kinematics == KINEMATICS_cleavage_opening_ID)
   if (maxNinstance == 0) return

src/kinematics_slipplane_opening.f90

@@ -51,7 +51,7 @@ subroutine kinematics_slipplane_opening_init
 
   integer :: maxNinstance,p,instance
 
-  write(6,'(/,a)')   ' <<<+-  kinematics_'//KINEMATICS_slipplane_opening_LABEL//' init  -+>>>'
+  write(6,'(/,a)') ' <<<+-  kinematics_'//KINEMATICS_slipplane_opening_LABEL//' init  -+>>>'; flush(6)
 
   maxNinstance = count(phase_kinematics == KINEMATICS_slipplane_opening_ID)
   if (maxNinstance == 0) return
@@ -134,51 +134,35 @@ subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc,
   dLd_dTstar = 0.0_pReal
   do i = 1, prm%totalNslip
 
-    projection_d = math_outer(prm%slip_direction(1:3,i),prm%slip_normal(1:3,i))
+    projection_d = math_outer(prm%slip_direction(1:3,i), prm%slip_normal(1:3,i))
     projection_t = math_outer(prm%slip_transverse(1:3,i),prm%slip_normal(1:3,i))
-    projection_n = math_outer(prm%slip_normal(1:3,i),prm%slip_normal(1:3,i))
+    projection_n = math_outer(prm%slip_normal(1:3,i),    prm%slip_normal(1:3,i))
     
     traction_d    = math_mul33xx33(S,projection_d)
     traction_t    = math_mul33xx33(S,projection_t)
     traction_n    = math_mul33xx33(S,projection_n)
     
-    traction_crit = prm%critLoad(i)* damage(homog)%p(damageOffset)                                                        ! degrading critical load carrying capacity by damage 
+    traction_crit = prm%critLoad(i)* damage(homog)%p(damageOffset)                                  ! degrading critical load carrying capacity by damage 
 
-    udotd = sign(1.0_pReal,traction_d)* &
+    udotd = sign(1.0_pReal,traction_d)* prm%sdot0* (  abs(traction_d)/traction_crit &
-      prm%sdot0* &
+                                                    - abs(traction_d)/prm%critLoad(i))**prm%n
-      (abs(traction_d)/traction_crit - &
+    udott = sign(1.0_pReal,traction_t)* prm%sdot0* (  abs(traction_t)/traction_crit &
-       abs(traction_d)/prm%critLoad(i))**prm%n
+                                                    - abs(traction_t)/prm%critLoad(i))**prm%n
-    if (abs(udotd) > tol_math_check) then
+    udotn = prm%sdot0* (  max(0.0_pReal,traction_n)/traction_crit &
-      Ld = Ld + udotd*projection_d
+                        - max(0.0_pReal,traction_n)/prm%critLoad(i))**prm%n
-      dudotd_dt = udotd*prm%n/traction_d
-      forall (k=1:3,l=1:3,m=1:3,n=1:3) &
-        dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + &
-          dudotd_dt*projection_d(k,l)*projection_d(m,n)
-    endif                
 
-    udott = sign(1.0_pReal,traction_t)* &
+    dudotd_dt = udotd*prm%n/traction_d
-      prm%sdot0* &
+    dudott_dt = udott*prm%n/traction_t
-      (abs(traction_t)/traction_crit - &
+    dudotn_dt = udotn*prm%n/traction_n
-       abs(traction_t)/prm%critLoad(i))**prm%n
-    if (abs(udott) > tol_math_check) then
-      Ld = Ld + udott*projection_t
-      dudott_dt = udott*prm%n/traction_t
-      forall (k=1:3,l=1:3,m=1:3,n=1:3) &
-        dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + &
-          dudott_dt*projection_t(k,l)*projection_t(m,n)
-    endif
 
-    udotn = &
+    forall (k=1:3,l=1:3,m=1:3,n=1:3) &
-      prm%sdot0* &
+      dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + dudotd_dt*projection_d(k,l)*projection_d(m,n) &
-      (max(0.0_pReal,traction_n)/traction_crit - &
+                                                + dudott_dt*projection_t(k,l)*projection_t(m,n) &
-       max(0.0_pReal,traction_n)/prm%critLoad(i))**prm%n
+                                                + dudotn_dt*projection_n(k,l)*projection_n(m,n)
-    if (abs(udotn) > tol_math_check) then
+
-      Ld = Ld + udotn*projection_n
+    Ld = Ld + udotd*projection_d &
-      dudotn_dt = udotn*prm%n/traction_n
+            + udott*projection_t &
-      forall (k=1:3,l=1:3,m=1:3,n=1:3) &
+            + udotn*projection_n
-        dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + &
-          dudotn_dt*projection_n(k,l)*projection_n(m,n)
-    endif 
   enddo
   
   end associate

src/kinematics_thermal_expansion.f90

@@ -42,7 +42,7 @@ subroutine kinematics_thermal_expansion_init
   real(pReal),  dimension(:), allocatable :: &
     temp
   
-  write(6,'(/,a)')   ' <<<+-  kinematics_'//KINEMATICS_thermal_expansion_LABEL//' init  -+>>>'
+  write(6,'(/,a)') ' <<<+-  kinematics_'//KINEMATICS_thermal_expansion_LABEL//' init  -+>>>'; flush(6)
  
   Ninstance = count(phase_kinematics == KINEMATICS_thermal_expansion_ID)
   

src/math.f90

@@ -398,19 +398,19 @@ pure function math_exp33(A,n)
   real(pReal), dimension(3,3), intent(in) :: A
   real(pReal), dimension(3,3) :: B, math_exp33
   real(pReal) :: invFac
-  integer     :: order
+  integer     :: n_
-
-  B = math_I3                                                                                       ! init
-  invFac = 1.0_pReal                                                                                ! 0!
-  math_exp33 = B                                                                                    ! A^0 = eye2
 
   if (present(n)) then
-    order = n
+    n_ = n
   else
-    order = 5
+    n_ = 5
   endif
 
-  do i = 1, order
+  invFac     = 1.0_pReal                                                                            ! 0!
+  B          = math_I3
+  math_exp33 = math_I3                                                                              ! A^0 = I
+
+  do i = 1, n_
     invFac = invFac/real(i,pReal)                                                                   ! invfac = 1/(i!)
     B = matmul(B,A)
     math_exp33 = math_exp33 + invFac*B                                                              ! exp = SUM (A^i)/(i!)
@@ -882,16 +882,20 @@ real(pReal) function math_sampleGaussVar(meanvalue, stddev, width)
   real(pReal), intent(in), optional ::  width                                                       ! width of considered values as multiples of standard deviation
   real(pReal), dimension(2) ::          rnd                                                         ! random numbers
   real(pReal) ::                        scatter, &                                                  ! normalized scatter around meanvalue
-                                        myWidth
+                                        width_
 
   if (abs(stddev) < tol_math_check) then
     math_sampleGaussVar = meanvalue
   else
-    myWidth = merge(width,3.0_pReal,present(width))                                                 ! use +-3*sigma as default value for scatter if not given
+    if (present(width)) then
+      width_ = width
+    else
+      width_ = 3.0_pReal                                                                            ! use +-3*sigma as default scatter
+    endif
 
     do
       call random_number(rnd)
-      scatter = myWidth * (2.0_pReal * rnd(1) - 1.0_pReal)
+      scatter = width_ * (2.0_pReal * rnd(1) - 1.0_pReal)
       if (rnd(2) <= exp(-0.5_pReal * scatter ** 2.0_pReal)) exit                                    ! test if scattered value is drawn
     enddo
 

src/mesh_grid.f90

@@ -296,7 +296,7 @@ end subroutine readGeom
 
 
 !---------------------------------------------------------------------------------------------------
-!> @brief Calculate undeformed position of IPs/cell centres (pretend to be an element)
+!> @brief Calculate undeformed position of IPs/cell centers (pretend to be an element)
 !---------------------------------------------------------------------------------------------------
 function IPcoordinates0(grid,geomSize,grid3Offset)
 

src/mesh_marc.f90

@@ -488,8 +488,7 @@ subroutine inputRead_mapNodes(fileContent)
     chunkPos = IO_stringPos(fileContent(l))
     if( IO_lc(IO_stringValue(fileContent(l),chunkPos,1)) == 'coordinates' ) then
       do i = 1,size(mesh_mapFEtoCPnode,2)
-        mesh_mapFEtoCPnode(1,i) = IO_fixedIntValue (fileContent(l+1+i),[0,10],1)
+        mesh_mapFEtoCPnode(1:2,i) = [IO_fixedIntValue (fileContent(l+1+i),[0,10],1),i]              ! ToDo: use IO_intValue
-        mesh_mapFEtoCPnode(2,i) = i
       enddo
       exit
     endif
@@ -520,9 +519,9 @@ subroutine inputRead_elemNodes(nodes, &
     chunkPos = IO_stringPos(fileContent(l))
     if( IO_lc(IO_stringValue(fileContent(l),chunkPos,1)) == 'coordinates' ) then
       do i=1,nNode
-        m = mesh_FEasCP('node',IO_fixedIntValue(fileContent(l+1+i),node_ends,1))
+        m = mesh_FEasCP('node',IO_fixedIntValue(fileContent(l+1+i),node_ends,1))                    !ToDo: use IO_intValue
         do j = 1,3
-          nodes(j,m) = mesh_unitlength * IO_fixedNoEFloatValue(fileContent(l+1+i),node_ends,j+1)
+          nodes(j,m) = mesh_unitlength * IO_fixedNoEFloatValue(fileContent(l+1+i),node_ends,j+1)    !ToDo: use IO_floatValue
         enddo
       enddo
       exit

src/quaternions.f90

@@ -1,37 +1,9 @@
-! ###################################################################
-! Copyright (c) 2013-2015, Marc De Graef/Carnegie Mellon University
-! Modified      2017-2019, Martin Diehl/Max-Planck-Institut für Eisenforschung GmbH
-! All rights reserved.
-!
-! Redistribution and use in source and binary forms, with or without modification, are
-! permitted provided that the following conditions are met:
-!
-!     - Redistributions of source code must retain the above copyright notice, this list
-!        of conditions and the following disclaimer.
-!     - Redistributions in binary form must reproduce the above copyright notice, this
-!        list of conditions and the following disclaimer in the documentation and/or
-!        other materials provided with the distribution.
-!     - Neither the names of Marc De Graef, Carnegie Mellon University nor the names
-!        of its contributors may be used to endorse or promote products derived from
-!        this software without specific prior written permission.
-!
-! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
-! AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
-! IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
-! ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
-! LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
-! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
-! SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
-! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
-! OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE
-! USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
-! ###################################################################
-
 !---------------------------------------------------------------------------------------------------
-!> @author Marc De Graef, Carnegie Mellon University
 !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
+!> @author Philip Eisenlohr, Michigan State University
 !> @brief general quaternion math, not limited to unit quaternions
-!> @details  w is the real part, (x, y, z) are the imaginary parts.
+!> @details w is the real part, (x, y, z) are the imaginary parts.
+!> @details https://en.wikipedia.org/wiki/Quaternion
 !---------------------------------------------------------------------------------------------------
 module quaternions
   use prec
@@ -78,15 +50,16 @@ module quaternions
 
     procedure, public  :: abs__
     procedure, public  :: dot_product__
-    procedure, public  :: conjg__
     procedure, public  :: exp__
     procedure, public  :: log__
-
+    procedure, public  :: conjg => conjg__
-    procedure, public  :: homomorphed => quat_homomorphed
-    procedure, public  :: asArray
     procedure, public  :: real  => real__
     procedure, public  :: aimag => aimag__
 
+    procedure, public  :: homomorphed
+    procedure, public  :: asArray
+    procedure, public  :: inverse
+
   end type
 
   interface assignment (=)
@@ -118,6 +91,14 @@ module quaternions
     module procedure log__
   end interface log
 
+  interface real
+    module procedure real__
+  end interface real
+
+  interface aimag
+    module procedure aimag__
+  end interface aimag
+  
   private :: &
     unitTest
 
@@ -125,49 +106,46 @@ contains
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief doing self test
+!> @brief do self test
 !--------------------------------------------------------------------------------------------------
 subroutine quaternions_init
 
-  write(6,'(/,a)')   ' <<<+-  quaternions init  -+>>>'
+  write(6,'(/,a)') ' <<<+-  quaternions init  -+>>>'; flush(6)
   call unitTest
 
 end subroutine quaternions_init
 
 
 !---------------------------------------------------------------------------------------------------
-!> constructor for a quaternion from a 4-vector
+!> construct a quaternion from a 4-vector
 !---------------------------------------------------------------------------------------------------
 type(quaternion) pure function init__(array)
 
   real(pReal), intent(in), dimension(4) :: array
 
-  init__%w=array(1)
+  init__%w = array(1)
-  init__%x=array(2)
+  init__%x = array(2)
-  init__%y=array(3)
+  init__%y = array(3)
-  init__%z=array(4)
+  init__%z = array(4)
 
 end function init__
 
 
 !---------------------------------------------------------------------------------------------------
-!> assing a quaternion
+!> assign a quaternion
 !---------------------------------------------------------------------------------------------------
 elemental pure subroutine assign_quat__(self,other)
 
   type(quaternion), intent(out) :: self
   type(quaternion), intent(in)  :: other
 
-  self%w = other%w
+  self = [other%w,other%x,other%y,other%z]
-  self%x = other%x
-  self%y = other%y
-  self%z = other%z
   
 end subroutine assign_quat__
 
 
 !---------------------------------------------------------------------------------------------------
-!> assing a 4-vector
+!> assign a 4-vector
 !---------------------------------------------------------------------------------------------------
 pure subroutine assign_vec__(self,other)
 
@@ -183,67 +161,57 @@ end subroutine assign_vec__
 
 
 !---------------------------------------------------------------------------------------------------
-!> addition of two quaternions
+!> add a quaternion
 !---------------------------------------------------------------------------------------------------
 type(quaternion) elemental pure function add__(self,other)
 
   class(quaternion), intent(in) :: self,other
 
-  add__%w = self%w + other%w
+  add__ = [ self%w,  self%x,  self%y ,self%z] &
-  add__%x = self%x + other%x
+        + [other%w, other%x, other%y,other%z]
-  add__%y = self%y + other%y
-  add__%z = self%z + other%z
   
 end function add__
 
 
 !---------------------------------------------------------------------------------------------------
-!> unary positive operator
+!> return (unary positive operator)
 !---------------------------------------------------------------------------------------------------
 type(quaternion) elemental pure function pos__(self)
 
   class(quaternion), intent(in) :: self
 
-  pos__%w = self%w
+  pos__ = self * (+1.0_pReal)
-  pos__%x = self%x
-  pos__%y = self%y
-  pos__%z = self%z
   
 end function pos__
 
 
 !---------------------------------------------------------------------------------------------------
-!> subtraction of two quaternions
+!> subtract a quaternion
 !---------------------------------------------------------------------------------------------------
 type(quaternion) elemental pure function sub__(self,other)
 
   class(quaternion), intent(in) :: self,other
 
-  sub__%w = self%w - other%w
+  sub__ = [ self%w,  self%x,  self%y ,self%z] &
-  sub__%x = self%x - other%x
+        - [other%w, other%x, other%y,other%z]
-  sub__%y = self%y - other%y
-  sub__%z = self%z - other%z
   
 end function sub__
 
 
 !---------------------------------------------------------------------------------------------------
-!> unary positive operator
+!> negate (unary negative operator)
 !---------------------------------------------------------------------------------------------------
 type(quaternion) elemental pure function neg__(self)
 
   class(quaternion), intent(in) :: self
 
-  neg__%w = -self%w
+  neg__ = self * (-1.0_pReal)
-  neg__%x = -self%x
-  neg__%y = -self%y
-  neg__%z = -self%z
   
 end function neg__
 
 
 !---------------------------------------------------------------------------------------------------
-!> multiplication of two quaternions
+!> multiply with a quaternion
 !---------------------------------------------------------------------------------------------------
 type(quaternion) elemental pure function mul_quat__(self,other)
 
@@ -258,23 +226,20 @@ end function mul_quat__
 
 
 !---------------------------------------------------------------------------------------------------
-!> multiplication of quaternions with scalar
+!> multiply with a scalar
 !---------------------------------------------------------------------------------------------------
 type(quaternion) elemental pure function mul_scal__(self,scal)
 
   class(quaternion), intent(in) :: self
-  real(pReal), intent(in) :: scal
+  real(pReal),       intent(in) :: scal
 
-  mul_scal__%w = self%w*scal
+  mul_scal__ = [self%w,self%x,self%y,self%z]*scal
-  mul_scal__%x = self%x*scal
-  mul_scal__%y = self%y*scal
-  mul_scal__%z = self%z*scal
   
 end function mul_scal__
 
 
 !---------------------------------------------------------------------------------------------------
-!> division of two quaternions
+!> divide by a quaternion
 !---------------------------------------------------------------------------------------------------
 type(quaternion) elemental pure function div_quat__(self,other)
 
@@ -286,12 +251,12 @@ end function div_quat__
 
 
 !---------------------------------------------------------------------------------------------------
-!> divisiont of quaternions by scalar
+!> divide by a scalar
 !---------------------------------------------------------------------------------------------------
 type(quaternion) elemental pure function div_scal__(self,scal)
 
   class(quaternion), intent(in) :: self
-  real(pReal), intent(in) :: scal
+  real(pReal),       intent(in) :: scal
 
   div_scal__ = [self%w,self%x,self%y,self%z]/scal
 
@@ -299,7 +264,7 @@ end function div_scal__
 
 
 !---------------------------------------------------------------------------------------------------
-!> equality of two quaternions
+!> test equality
 !---------------------------------------------------------------------------------------------------
 logical elemental pure function eq__(self,other)
 
@@ -312,7 +277,7 @@ end function eq__
 
 
 !---------------------------------------------------------------------------------------------------
-!> inequality of two quaternions
+!> test inequality
 !---------------------------------------------------------------------------------------------------
 logical elemental pure function neq__(self,other)
 
@@ -324,20 +289,7 @@ end function neq__
 
 
 !---------------------------------------------------------------------------------------------------
-!> quaternion to the power of a scalar
+!> raise to the power of a quaternion
-!---------------------------------------------------------------------------------------------------
-type(quaternion) elemental pure function pow_scal__(self,expon)
-
-  class(quaternion), intent(in) :: self
-  real(pReal), intent(in) :: expon
-
-  pow_scal__ = exp(log(self)*expon)
-
-end function pow_scal__
-
-
-!---------------------------------------------------------------------------------------------------
-!> quaternion to the power of a quaternion
 !---------------------------------------------------------------------------------------------------
 type(quaternion) elemental pure function pow_quat__(self,expon)
 
@@ -350,45 +302,60 @@ end function pow_quat__
 
 
 !---------------------------------------------------------------------------------------------------
-!> exponential of a quaternion
+!> raise to the power of a scalar
-!> ToDo: Lacks any check for invalid operations
+!---------------------------------------------------------------------------------------------------
+type(quaternion) elemental pure function pow_scal__(self,expon)
+
+  class(quaternion), intent(in) :: self
+  real(pReal),       intent(in) :: expon
+
+  pow_scal__ = exp(log(self)*expon)
+
+end function pow_scal__
+
+
+!---------------------------------------------------------------------------------------------------
+!> take exponential
 !---------------------------------------------------------------------------------------------------
 type(quaternion) elemental pure function exp__(self)
 
   class(quaternion), intent(in) :: self
   real(pReal)                   :: absImag
 
-  absImag = norm2([self%x, self%y, self%z])
+  absImag = norm2(aimag(self))
 
-  exp__ = exp(self%w) * [                 cos(absImag), &
+  exp__ = merge(exp(self%w) * [                 cos(absImag), &
-                         self%x/absImag * sin(absImag), &
+                               self%x/absImag * sin(absImag), &
-                         self%y/absImag * sin(absImag), &
+                               self%y/absImag * sin(absImag), &
-                         self%z/absImag * sin(absImag)]
+                               self%z/absImag * sin(absImag)], &
+                IEEE_value(1.0_pReal,IEEE_SIGNALING_NAN), &
+                dNeq0(absImag))
 
 end function exp__
 
 
 !---------------------------------------------------------------------------------------------------
-!> logarithm of a quaternion
+!> take logarithm
-!> ToDo: Lacks any check for invalid operations
 !---------------------------------------------------------------------------------------------------
 type(quaternion) elemental pure function log__(self)
 
   class(quaternion), intent(in) :: self
   real(pReal)                   :: absImag
 
-  absImag = norm2([self%x, self%y, self%z])
+  absImag = norm2(aimag(self))
 
-  log__ = [log(abs(self)), &
+  log__ = merge([log(abs(self)), &
-           self%x/absImag * acos(self%w/abs(self)), &
+                 self%x/absImag * acos(self%w/abs(self)), &
-           self%y/absImag * acos(self%w/abs(self)), &
+                 self%y/absImag * acos(self%w/abs(self)), &
-           self%z/absImag * acos(self%w/abs(self))]
+                 self%z/absImag * acos(self%w/abs(self))], &
+                IEEE_value(1.0_pReal,IEEE_SIGNALING_NAN), &
+                dNeq0(absImag))
 
 end function log__
 
 
 !---------------------------------------------------------------------------------------------------
-!> norm of a quaternion
+!> return norm
 !---------------------------------------------------------------------------------------------------
 real(pReal) elemental pure function abs__(a)
 
@@ -400,7 +367,7 @@ end function abs__
 
 
 !---------------------------------------------------------------------------------------------------
-!> dot product of two quaternions
+!> calculate dot product
 !---------------------------------------------------------------------------------------------------
 real(pReal) elemental pure function dot_product__(a,b)
 
@@ -412,31 +379,31 @@ end function dot_product__
 
 
 !---------------------------------------------------------------------------------------------------
-!> conjugate complex of a quaternion
+!> take conjugate complex
 !---------------------------------------------------------------------------------------------------
 type(quaternion) elemental pure function conjg__(a)
 
   class(quaternion), intent(in) :: a
 
-  conjg__ = quaternion([a%w, -a%x, -a%y, -a%z])
+  conjg__ = [a%w, -a%x, -a%y, -a%z]
 
 end function conjg__
 
 
 !---------------------------------------------------------------------------------------------------
-!> homomorphed quaternion of a quaternion
+!> homomorph
 !---------------------------------------------------------------------------------------------------
-type(quaternion) elemental pure function quat_homomorphed(self)
+type(quaternion) elemental pure function homomorphed(self)
 
   class(quaternion), intent(in) :: self
 
-  quat_homomorphed = quaternion(-[self%w,self%x,self%y,self%z])
+  homomorphed = - self
 
-end function quat_homomorphed
+end function homomorphed
 
 
 !---------------------------------------------------------------------------------------------------
-!> quaternion as plain array
+!> return as plain array
 !---------------------------------------------------------------------------------------------------
 pure function asArray(self)
 
@@ -449,7 +416,7 @@ end function asArray
 
 
 !---------------------------------------------------------------------------------------------------
-!> real part of a quaternion
+!> real part (scalar)
 !---------------------------------------------------------------------------------------------------
 pure function real__(self)
 
@@ -462,7 +429,7 @@ end function real__
 
 
 !---------------------------------------------------------------------------------------------------
-!> imaginary part of a quaternion
+!> imaginary part (3-vector)
 !---------------------------------------------------------------------------------------------------
 pure function aimag__(self)
 
@@ -474,45 +441,86 @@ pure function aimag__(self)
 end function aimag__
 
 
+!---------------------------------------------------------------------------------------------------
+!> inverse
+!---------------------------------------------------------------------------------------------------
+type(quaternion) elemental pure function inverse(self)
+
+  class(quaternion), intent(in) :: self
+
+  inverse = conjg(self)/abs(self)**2.0_pReal
+
+end function inverse
+
+
 !--------------------------------------------------------------------------------------------------
 !> @brief check correctness of (some) quaternions functions
 !--------------------------------------------------------------------------------------------------
 subroutine unitTest
 
   real(pReal), dimension(4) :: qu
-
   type(quaternion)          :: q, q_2
 
   call random_number(qu)
-  q = qu
+  qu = (qu-0.5_pReal) * 2.0_pReal
+  q  = quaternion(qu)
+  
+  q_2= qu
+  if(any(dNeq(q%asArray(),q_2%asArray())))             call IO_error(401,ext_msg='assign_vec__')
 
   q_2 = q + q
-  if(any(dNeq(q_2%asArray(),2.0_pReal*qu)))         call IO_error(401,ext_msg='add__')
+  if(any(dNeq(q_2%asArray(),2.0_pReal*qu)))            call IO_error(401,ext_msg='add__')
 
   q_2 = q - q
-  if(any(dNeq0(q_2%asArray())))                     call IO_error(401,ext_msg='sub__')
+  if(any(dNeq0(q_2%asArray())))                        call IO_error(401,ext_msg='sub__')
 
-  q_2 = q * 5.0_preal
+  q_2 = q * 5.0_pReal
-  if(any(dNeq(q_2%asArray(),5.0_pReal*qu)))         call IO_error(401,ext_msg='mul__')
+  if(any(dNeq(q_2%asArray(),5.0_pReal*qu)))            call IO_error(401,ext_msg='mul__')
 
-  q_2 = q / 0.5_preal
+  q_2 = q / 0.5_pReal
-  if(any(dNeq(q_2%asArray(),2.0_pReal*qu)))         call IO_error(401,ext_msg='div__')
+  if(any(dNeq(q_2%asArray(),2.0_pReal*qu)))            call IO_error(401,ext_msg='div__')
+  
+  q_2 = q * 0.3_pReal
+  if(dNeq0(abs(q)) .and. q_2 == q)                     call IO_error(401,ext_msg='eq__')
   
   q_2 = q
-  if(q_2 /= q)                                      call IO_error(401,ext_msg='eq__')
+  if(q_2 /= q)                                         call IO_error(401,ext_msg='neq__')
 
-  if(any(dNeq(q%asArray(),qu)))                     call IO_error(401,ext_msg='eq__')
+  if(dNeq(abs(q),norm2(qu)))                           call IO_error(401,ext_msg='abs__')
-  if(dNeq(q%real(),       qu(1)))                   call IO_error(401,ext_msg='real()')
+  if(dNeq(abs(q)**2.0_pReal, real(q*q%conjg()),1.0e-14_pReal)) &
-  if(any(dNeq(q%aimag(),  qu(2:4))))                call IO_error(401,ext_msg='aimag()')
+                                                       call IO_error(401,ext_msg='abs__/*conjg')
+
+  if(any(dNeq(q%asArray(),qu)))                        call IO_error(401,ext_msg='eq__')
+  if(dNeq(q%real(),       qu(1)))                      call IO_error(401,ext_msg='real()')
+  if(any(dNeq(q%aimag(),  qu(2:4))))                   call IO_error(401,ext_msg='aimag()')
 
   q_2 = q%homomorphed()
-  if(q                 /= q_2*    (-1.0_pReal))     call IO_error(401,ext_msg='homomorphed')
+  if(q                 /= q_2*    (-1.0_pReal))        call IO_error(401,ext_msg='homomorphed')
-  if(dNeq(q_2%real(),     qu(1)*  (-1.0_pReal)))    call IO_error(401,ext_msg='homomorphed/real')
+  if(dNeq(q_2%real(),     qu(1)*  (-1.0_pReal)))       call IO_error(401,ext_msg='homomorphed/real')
-  if(any(dNeq(q_2%aimag(),qu(2:4)*(-1.0_pReal))))   call IO_error(401,ext_msg='homomorphed/aimag')
+  if(any(dNeq(q_2%aimag(),qu(2:4)*(-1.0_pReal))))      call IO_error(401,ext_msg='homomorphed/aimag')
 
   q_2 = conjg(q)
-  if(dNeq(q_2%real(),     q%real()))                call IO_error(401,ext_msg='conjg/real')
+  if(dNeq(abs(q),abs(q_2)))                            call IO_error(401,ext_msg='conjg/abs')
-  if(any(dNeq(q_2%aimag(),q%aimag()*(-1.0_pReal)))) call IO_error(401,ext_msg='conjg/aimag')
+  if(q /= conjg(q_2))                                  call IO_error(401,ext_msg='conjg/involution')
+  if(dNeq(q_2%real(),     q%real()))                   call IO_error(401,ext_msg='conjg/real')
+  if(any(dNeq(q_2%aimag(),q%aimag()*(-1.0_pReal))))    call IO_error(401,ext_msg='conjg/aimag')
+  
+  if(abs(q) > 0.0_pReal) then
+    q_2 = q * q%inverse()
+    if(     dNeq(real(q_2), 1.0_pReal,1.0e-15_pReal))  call IO_error(401,ext_msg='inverse/real')
+    if(any(dNeq0(aimag(q_2),          1.0e-15_pReal))) call IO_error(401,ext_msg='inverse/aimag')
+    
+    q_2 = q/abs(q)
+    q_2 = conjg(q_2) - inverse(q_2)
+    if(any(dNeq0(q_2%asArray(),1.0e-15_pReal)))        call IO_error(401,ext_msg='inverse/conjg')
+  endif
+
+#if !(defined(__GFORTRAN__) &&  __GNUC__ < 9)
+  if (norm2(aimag(q)) > 0.0_pReal) then
+    if (dNeq0(abs(q-exp(log(q))),1.0e-13_pReal))       call IO_error(401,ext_msg='exp/log')
+    if (dNeq0(abs(q-log(exp(q))),1.0e-13_pReal))       call IO_error(401,ext_msg='log/exp')
+  endif
+#endif
 
 end subroutine unitTest
 

src/results.f90

@@ -47,6 +47,7 @@ module results
     results_openJobFile, &
     results_closeJobFile, &
     results_addIncrement, &
+    results_finalizeIncrement, &
     results_addGroup, &
     results_openGroup, &
     results_closeGroup, &
@@ -119,6 +120,17 @@ subroutine results_addIncrement(inc,time)
 end subroutine results_addIncrement
 
 
+!--------------------------------------------------------------------------------------------------
+!> @brief finalize increment
+!> @details remove soft link
+!--------------------------------------------------------------------------------------------------
+subroutine results_finalizeIncrement
+
+  call results_removeLink('current')
+
+end subroutine results_finalizeIncrement
+
+
 !--------------------------------------------------------------------------------------------------
 !> @brief open a group from the results file
 !--------------------------------------------------------------------------------------------------

src/system_routines.f90

@@ -93,21 +93,24 @@ end function isDirectory
 !--------------------------------------------------------------------------------------------------
 !> @brief gets the current working directory
 !--------------------------------------------------------------------------------------------------
-character(len=1024) function getCWD()
+function getCWD()
 
   character(kind=C_CHAR), dimension(1024) :: charArray                                              ! C string is an array
+  character(len=:),       allocatable     :: getCWD
   integer(C_INT) :: stat
   integer :: i
  
   call getCurrentWorkDir_C(charArray,stat)
+  
   if (stat /= 0_C_INT) then
     getCWD = 'Error occured when getting currend working directory'
   else
-    getCWD = repeat('',len(getCWD))
+    allocate(character(len=1024)::getCWD)
     arrayToString: do i=1,len(getCWD)
       if (charArray(i) /= C_NULL_CHAR) then
         getCWD(i:i)=charArray(i)
       else
+        getCWD = getCWD(:i-1)
         exit
       endif
     enddo arrayToString
@@ -119,21 +122,24 @@ end function getCWD
 !--------------------------------------------------------------------------------------------------
 !> @brief gets the current host name
 !--------------------------------------------------------------------------------------------------
-character(len=1024) function getHostName()
+function getHostName()
 
   character(kind=C_CHAR), dimension(1024) :: charArray                                              ! C string is an array
+  character(len=:),       allocatable     :: getHostName
   integer(C_INT) :: stat
   integer :: i
  
   call getHostName_C(charArray,stat)
+
   if (stat /= 0_C_INT) then
     getHostName = 'Error occured when getting host name'
   else
-    getHostName = repeat('',len(getHostName))
+    allocate(character(len=1024)::getHostName)
     arrayToString: do i=1,len(getHostName)
       if (charArray(i) /= C_NULL_CHAR) then
         getHostName(i:i)=charArray(i)
       else
+        getHostName = getHostName(:i-1)
         exit
       endif
     enddo arrayToString