new version of AL spectral solver (seems to work, but still experimental)

new concept of DAMASK for Python: Now using also IO.f90, debug.f90 etc. instead of mimicing their necessary function. This needs DAMASK_python.f90 and DAMASK_python_interface.f90 (At the moment more or less copies of respective spectral files).
polished and renamed the scripts for converging ang files
renamed voronoi_randomSeeding.py to spectral_randomSeeding.py

code/DAMASK_python.f90

@@ -0,0 +1,40 @@
+! Copyright 2012 Max-Planck-Institut für Eisenforschung GmbH
+!
+! This file is part of DAMASK,
+! the Düsseldorf Advanced Material Simulation Kit.
+!
+! DAMASK is free software: you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation, either version 3 of the License, or
+! (at your option) any later version.
+!
+! DAMASK is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License
+! along with DAMASK. If not, see <http://www.gnu.org/licenses/>.
+!
+!##################################################################################################
+!* $Id$
+!##################################################################################################
+
+!********************************************************************
+! quit subroutine to satisfy IO_error
+!
+!********************************************************************
+subroutine quit(stop_id)
+ use prec, only: &
+   pInt
+   
+ implicit none
+ integer(pInt), intent(in) :: stop_id
+ 
+ if (stop_id ==    0_pInt) stop 0                                                                   ! normal termination
+ if (stop_id <= 9000_pInt) then                                                                     ! trigger regridding
+   write(6,'(i4)') stop_id
+   stop 1
+ endif
+ stop 'abnormal termination of DAMASK_spectral'
+end subroutine

code/DAMASK_python_interface.f90

@@ -0,0 +1,392 @@
+! Copyright 2011 Max-Planck-Institut für Eisenforschung GmbH
+!
+! This file is part of DAMASK,
+! the Düsseldorf Advanced Material Simulation Kit.
+!
+! DAMASK is free software: you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation, either version 3 of the License, or
+! (at your option) any later version.
+!
+! DAMASK is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License
+! along with DAMASK. If not, see <http://www.gnu.org/licenses/>.
+!
+!--------------------------------------------------------------------------------------------------
+!* $Id$
+!--------------------------------------------------------------------------------------------------
+!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
+!> @brief Interfacing between the spectral solver and the material subroutines provided
+!! by DAMASK
+!--------------------------------------------------------------------------------------------------
+module DAMASK_interface
+
+ implicit none
+ private
+ character(len=64),   parameter, public  :: FEsolver = 'Spectral'                                   !< Keyword for spectral solver
+ character(len=5),    parameter, public  :: inputFileExtension = '.geom'                            !< File extension for geometry description
+ character(len=4),    parameter, public  :: logFileExtension = '.log'                               !< Dummy variable as the spectral solver has no log
+ character(len=1024),            private :: geometryParameter, &                                    !< Interpretated parameter given at command line
+                                            loadcaseParameter                                       !< Interpretated parameter given at command line
+
+ public  :: getSolverWorkingDirectoryName, & !< Interpretated parameter given at command line
+            getSolverJobName, &
+            getLoadCase, &
+            getLoadCaseName, &
+            getModelName, &
+            DAMASK_interface_init 
+ private :: rectifyPath, &
+            makeRelativePath, &
+            getPathSep
+
+contains
+
+!--------------------------------------------------------------------------------------------------
+!> @brief initializes the solver by interpreting the command line arguments. Also writes
+!! information on computation on screen
+!--------------------------------------------------------------------------------------------------
+subroutine DAMASK_interface_init
+ use, intrinsic :: iso_fortran_env                                                                  ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
+ use prec,   only: pInt
+
+ implicit none
+ character(len=1024)   :: commandLine, &                                                            !< command line call as string
+                          hostName, &                                                               !< name of computer
+                          userName                                                                  !< name of user calling the executable
+ integer               :: i, &
+                          start ,&
+                          length
+ integer, dimension(8) :: dateAndTime                                                               ! type default integer
+ 
+ call get_command(commandLine)
+ call date_and_time(values = dateAndTime)
+ do i = 1,len(commandLine)                                                                          ! remove capitals
+   if(64<iachar(commandLine(i:i)) .and. iachar(commandLine(i:i))<91) & 
+                   commandLine(i:i) = achar(iachar(commandLine(i:i))+32)
+ enddo
+
+ if(index(commandLine,' -h ',.true.) > 0 .or. index(commandLine,' --help ',.true.) > 0) then        ! search for ' -h ' or '--help'
+   write(6,'(a)') '$Id$'
+#include "compilation_info.f90"
+   write(6,'(a)')  '#############################################################'
+   write(6,'(a)')  'DAMASK spectral:'
+   write(6,'(a)')  'The spectral method boundary value problem solver for'
+   write(6,'(a)')  'the Duesseldorf Advanced Material Simulation Kit'
+   write(6,'(a)')  '#############################################################'
+   write(6,'(a)')  'Valid command line switches:'
+   write(6,'(a)')  '   --geom    (-g, --geometry)'
+   write(6,'(a)')  '   --load    (-l, --loadcase)'
+   write(6,'(a)')  '   --restart (-r)'
+   write(6,'(a)')  '   --help    (-h)'
+   write(6,'(a)')  ' '
+   write(6,'(a)')  'Mandatory Arguments:'
+   write(6,'(a)')  '  --load PathToLoadFile/NameOfLoadFile.load'
+   write(6,'(a)')  '       "PathToGeomFile" will be the working directory.'
+   write(6,'(a)')  '       Make sure the file "material.config" exists in the working'
+   write(6,'(a)')  '           directory'   
+   write(6,'(a)')  '       For further configuration place "numerics.config"'
+   write(6,'(a)')  '           and "numerics.config" in that directory.'
+   write(6,'(a)')  ' '
+   write(6,'(a)')  '  --geom PathToGeomFile/NameOfGeom.geom'
+   write(6,'(a)')  ' '
+   write(6,'(a)')  'Optional Argument:'
+   write(6,'(a)')  '  --restart XX'
+   write(6,'(a)')  '       Reads in total increment No. XX-1 and continous to'
+   write(6,'(a)')  '           calculate total increment No. XX.'
+   write(6,'(a)')  '       Attention: Overwrites existing results file '
+   write(6,'(a)')  '           "NameOfGeom_NameOfLoadFile_spectralOut".'
+   write(6,'(a)')  '       Works only if the restart information for total increment'
+   write(6,'(a)')  '            No. XX-1 is available in the working directory.'
+   write(6,'(a)')  'Help:'
+   write(6,'(a)')  '  --help'
+   write(6,'(a)')  '       Prints this message and exits'
+   write(6,'(a)')  ' '
+   call quit(0_pInt)
+ endif
+ if (.not.(command_argument_count()==4 .or. command_argument_count()==6)) &                         ! check for correct number of given arguments (no --help)
+   stop 'Wrong Nr. of Arguments. Run DAMASK_spectral.exe --help'                                    ! Could not find valid keyword (position 0 +3). Functions from IO.f90 are not available
+ start = index(commandLine,'-g',.true.) + 3                                                         ! search for '-g' and jump to first char of geometry
+ if (index(commandLine,'--geom',.true.)>0) then                                                     ! if '--geom' is found, use that (contains '-g')
+   start = index(commandLine,'--geom',.true.) + 7
+ endif               
+ if (index(commandLine,'--geometry',.true.)>0) then                                                 ! again, now searching for --geometry'
+   start = index(commandLine,'--geometry',.true.) + 11
+ endif
+ if(start==3_pInt) then                                                                             ! Could not find valid keyword (position 0 +3). Functions from IO.f90 are not available
+   write(6,'(a)') 'No Geometry specified'
+   call quit(9999)
+ endif
+ length = index(commandLine(start:len(commandLine)),' ',.false.)
+
+ call get_command(commandLine)                                                                      ! may contain capitals
+ geometryParameter = ''                                                                             ! should be empty
+ geometryParameter(1:length)=commandLine(start:start+length)
+
+ do i=1,len(commandLine)                                                                            ! remove capitals
+   if(64<iachar(commandLine(i:i)) .and. iachar(commandLine(i:i))<91) commandLine(i:i)&
+                                                               = achar(iachar(commandLine(i:i))+32)
+ enddo
+ 
+ start = index(commandLine,'-l',.true.) + 3                                                         ! search for '-l' and jump forward iby 3 to given name
+ if (index(commandLine,'--load',.true.)>0) then                                                     ! if '--load' is found, use that (contains '-l')
+   start = index(commandLine,'--load',.true.) + 7
+ endif               
+ if (index(commandLine,'--loadcase',.true.)>0) then                                                 ! again, now searching for --loadcase'
+   start = index(commandLine,'--loadcase',.true.) + 11
+ endif
+ if(start==3_pInt) then                                                                             ! Could not find valid keyword (position 0 +3). Functions from IO.f90 are not available
+   write(6,'(a)') 'No Loadcase specified'
+   call quit(9999)
+ endif
+ length = index(commandLine(start:len(commandLine)),' ',.false.)
+
+ call get_command(commandLine)                                                                      ! may contain capitals
+ loadcaseParameter = ''                                                                             ! should be empty
+ loadcaseParameter(1:length)=commandLine(start:start+length)
+ 
+ do i=1,len(commandLine)                                                                            ! remove capitals
+   if(64<iachar(commandLine(i:i)) .and. iachar(commandLine(i:i))<91) commandLine(i:i)&
+                                                               = achar(iachar(commandLine(i:i))+32)
+ enddo
+
+ start = index(commandLine,'-r',.true.) + 3                                                         ! search for '-r' and jump forward iby 3 to given name
+ if (index(commandLine,'--restart',.true.)>0) then                                                  ! if '--restart' is found, use that (contains '-l')
+   start = index(commandLine,'--restart',.true.) + 7
+ endif 
+ length = index(commandLine(start:len(commandLine)),' ',.false.)
+
+ call get_command(commandLine)                                                                      ! may contain capitals
+ call GET_ENVIRONMENT_VARIABLE('HOST',hostName)
+ call GET_ENVIRONMENT_VARIABLE('USER',userName)
+
+ write(6,*)
+ write(6,*) '<<<+-  DAMASK_spectral_interface init  -+>>>'
+ write(6,*) '$Id$'
+#include "compilation_info.f90"
+ write(6,'(a,2(i2.2,a),i4.4)') ' Date:               ',dateAndTime(3),'/',&
+                                                       dateAndTime(2),'/',&
+                                                       dateAndTime(1) 
+ write(6,'(a,2(i2.2,a),i2.2)') ' Time:               ',dateAndTime(5),':',&
+                                                       dateAndTime(6),':',&
+                                                       dateAndTime(7)  
+ write(6,*) 'Host Name:          ', trim(hostName)
+ write(6,*) 'User Name:          ', trim(userName)
+ write(6,*) 'Path Separator:     ', getPathSep()
+ write(6,*) 'Command line call:  ', trim(commandLine)
+ write(6,*) 'Geometry Parameter: ', trim(geometryParameter)
+ write(6,*) 'Loadcase Parameter: ', trim(loadcaseParameter)
+ if (start/=3_pInt) write(6,*) 'Restart Parameter:  ', trim(commandLine(start:start+length))
+
+end subroutine DAMASK_interface_init
+
+!--------------------------------------------------------------------------------------------------
+!> @brief  extract working directory from loadcase file possibly based on current working dir
+!--------------------------------------------------------------------------------------------------
+ character(len=1024) function getSolverWorkingDirectoryName()
+
+ implicit none
+ character(len=1024) :: cwd
+ character           :: pathSep
+
+ pathSep = getPathSep()
+
+ if (geometryParameter(1:1) == pathSep) then                                                        ! absolute path given as command line argument
+   getSolverWorkingDirectoryName = geometryParameter(1:scan(geometryParameter,pathSep,back=.true.))
+ else
+   call getcwd(cwd)
+   getSolverWorkingDirectoryName = trim(cwd)//pathSep//geometryParameter(1:scan(geometryParameter,pathSep,back=.true.))
+ endif
+
+ getSolverWorkingDirectoryName = rectifyPath(getSolverWorkingDirectoryName)
+ 
+end function getSolverWorkingDirectoryName
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief  basename of geometry file from command line arguments
+!--------------------------------------------------------------------------------------------------
+character(len=1024) function getSolverJobName()
+
+ implicit none
+ getSolverJobName = trim(getModelName())//'_'//trim(getLoadCase())
+
+end function getSolverJobName
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief  basename of geometry file from command line arguments
+!--------------------------------------------------------------------------------------------------
+character(len=1024) function getModelName()
+
+ use prec, only: pInt
+
+ implicit none
+ character(len=1024) :: cwd
+ integer :: posExt,posSep
+ character :: pathSep
+
+ pathSep = getPathSep()
+ posExt = scan(geometryParameter,'.',back=.true.)
+ posSep = scan(geometryParameter,pathSep,back=.true.)
+
+ if (posExt <= posSep) posExt = len_trim(geometryParameter)+1                                       ! no extension present
+ getModelName = geometryParameter(1:posExt-1_pInt)                                                  ! path to geometry file (excl. extension)
+
+ if (scan(getModelName,pathSep) /= 1) then                                                          ! relative path given as command line argument
+   call getcwd(cwd)
+   getModelName = rectifyPath(trim(cwd)//'/'//getModelName)
+ else
+   getModelName = rectifyPath(getModelName)
+ endif
+
+ getModelName = makeRelativePath(getSolverWorkingDirectoryName(),&
+                                 getModelName)
+                                 
+end function getModelName
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief  name of load case file exluding extension
+!--------------------------------------------------------------------------------------------------
+character(len=1024) function getLoadCase()
+
+ implicit none
+ integer :: posExt,posSep
+ character :: pathSep
+
+ pathSep = getPathSep()
+ posExt = scan(loadcaseParameter,'.',back=.true.)
+ posSep = scan(loadcaseParameter,pathSep,back=.true.)
+
+ if (posExt <= posSep) posExt = len_trim(loadcaseParameter)+1                                       ! no extension present
+ getLoadCase = loadcaseParameter(posSep+1:posExt-1)                                                 ! name of load case file exluding extension
+
+end function getLoadCase
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief relative path of loadcase from command line arguments
+!--------------------------------------------------------------------------------------------------
+character(len=1024) function getLoadcaseName()
+
+ implicit none
+ character(len=1024) :: cwd
+ integer :: posExt = 0, posSep
+ character :: pathSep
+
+ pathSep = getPathSep()
+ getLoadcaseName = loadcaseParameter
+ posExt = scan(getLoadcaseName,'.',back=.true.)
+ posSep = scan(getLoadcaseName,pathSep,back=.true.)
+
+ if (posExt <= posSep) getLoadcaseName = trim(getLoadcaseName)//('.load')                           ! no extension present
+ if (scan(getLoadcaseName,pathSep) /= 1) then                                                       ! relative path given as command line argument
+   call getcwd(cwd)
+   getLoadcaseName = rectifyPath(trim(cwd)//pathSep//getLoadcaseName)
+ else
+   getLoadcaseName = rectifyPath(getLoadcaseName)
+ endif
+
+ getLoadcaseName = makeRelativePath(getSolverWorkingDirectoryName(),&
+                                    getLoadcaseName)
+end function getLoadcaseName
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief remove ../ and ./ from path
+!--------------------------------------------------------------------------------------------------
+function rectifyPath(path)
+
+ implicit none
+ character(len=*) :: path
+ character(len=len_trim(path)) :: rectifyPath
+ character :: pathSep
+ integer :: i,j,k,l !no pInt
+
+ pathSep = getPathSep()
+
+ !remove ./ from path
+ l = len_trim(path)
+ rectifyPath = path
+ do i = l,3,-1
+    if ( rectifyPath(i-1:i) == '.'//pathSep .and. rectifyPath(i-2:i-2) /= '.' ) &
+      rectifyPath(i-1:l) = rectifyPath(i+1:l)//'  '
+ enddo
+
+ !remove ../ and corresponding directory from rectifyPath
+ l = len_trim(rectifyPath)
+ i = index(rectifyPath(i:l),'..'//pathSep)
+ j = 0
+ do while (i > j)
+    j = scan(rectifyPath(1:i-2),pathSep,back=.true.)
+    rectifyPath(j+1:l) = rectifyPath(i+3:l)//repeat(' ',2+i-j)
+    if (rectifyPath(j+1:j+1) == pathSep) then                                                       !search for '//' that appear in case of XXX/../../XXX
+      k = len_trim(rectifyPath)
+      rectifyPath(j+1:k-1) = rectifyPath(j+2:k)
+      rectifyPath(k:k) = ' '
+    endif
+    i = j+index(rectifyPath(j+1:l),'..'//pathSep)
+ enddo
+ if(len_trim(rectifyPath) == 0) rectifyPath = pathSep
+
+end function rectifyPath
+
+ 
+!--------------------------------------------------------------------------------------------------
+!> @brief relative path from absolute a to absolute b
+!--------------------------------------------------------------------------------------------------
+character(len=1024) function makeRelativePath(a,b)
+
+ implicit none
+ character (len=*) :: a,b
+ character :: pathSep
+ integer :: i,posLastCommonSlash,remainingSlashes !no pInt
+
+ pathSep = getPathSep()
+ posLastCommonSlash = 0
+ remainingSlashes = 0
+
+ do i = 1, min(1024,len_trim(a),len_trim(b))
+   if (a(i:i) /= b(i:i)) exit
+   if (a(i:i) == pathSep) posLastCommonSlash = i
+ enddo
+ do i = posLastCommonSlash+1,len_trim(a)
+   if (a(i:i) == pathSep) remainingSlashes = remainingSlashes + 1
+ enddo
+ makeRelativePath = repeat('..'//pathSep,remainingSlashes)//b(posLastCommonSlash+1:len_trim(b))
+
+end function makeRelativePath
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief counting / and \ in $PATH System variable the character occuring more often is assumed
+!! to be the path separator
+!--------------------------------------------------------------------------------------------------
+character function getPathSep()
+
+ use prec, only: pInt
+
+ implicit none
+ character(len=2048) path
+ integer(pInt) :: backslash = 0_pInt, slash = 0_pInt
+ integer :: i
+
+ call get_environment_variable('PATH',path)
+ do i=1, len(trim(path))
+   if (path(i:i)=='/') slash     =     slash + 1_pInt
+   if (path(i:i)=='\') backslash = backslash + 1_pInt
+ enddo
+
+ if (backslash>slash) then
+   getPathSep = '\'
+ else
+   getPathSep = '/'
+ endif
+
+end function
+
+end module

code/DAMASK_spectral.f90

@@ -358,7 +358,7 @@ program DAMASK_spectral
      end select
  enddo; enddo
 101 close(myUnit)
-if (sum(bc(1:N_Loadcases)%incs)>9000_pInt) stop !discuss with Philip, stop code trouble. suggesting warning
+if (sum(bc(1:N_Loadcases)%incs)>9000_pInt) stop 'to many incs' !discuss with Philip, stop code trouble. suggesting warning
 
 !-------------------------------------------------------------------------------------------------- ToDo: if temperature at CPFEM is treated properly, move this up immediately after interface init
 ! initialization of all related DAMASK modules (e.g. mesh.f90 reads in geometry)
@@ -961,7 +961,7 @@ C_ref = C * wgt
            write(6,'(a,es11.4)')        'error divergence Real max = ',err_real_div_max
          endif
          write(6,'(a,f6.2,a,es11.4,a)') 'error divergence = ', err_div/err_div_tol,&
-                                                           ' (',err_div_RMS,' N/m³)'
+                                                           ' (',err_div,' N/m³)'
 
 !--------------------------------------------------------------------------------------------------
 ! to the actual spectral method calculation (mechanical equilibrium)

code/DAMASK_spectral_AL.f90

@@ -17,7 +17,7 @@
 ! along with DAMASK. If not, see <http://www.gnu.org/licenses/>.
 !
 !##################################################################################################
-!* $Id: DAMASK_spectral.f90 1321 2012-02-15 18:58:38Z MPIE\u.diehl $
+!* $Id$
 !##################################################################################################
 ! Material subroutine for BVP solution using spectral method
 !
@@ -111,7 +111,8 @@ program DAMASK_spectral_AL
 
 !--------------------------------------------------------------------------------------------------
 ! stress, stiffness and compliance average etc.
- real(pReal), dimension(3,3) ::           P_av, P, F_aim = math_I3, F_aim_lastInc = math_I3,&
+ real(pReal), dimension(3,3) ::           P_av = 0.0_pReal, P_star_av = 0.0_pReal, P, &
+                                          F_aim = math_I3, F_aim_lastInc = math_I3, lambda_av, &
                                           mask_stress, mask_defgrad, deltaF, F_star_av, &
                                           F_aim_lab                                                 ! quantities rotated to other coordinate system
  real(pReal), dimension(3,3,3,3) ::       dPdF, C_inc0, C=0.0_pReal, S_lastInc, C_lastInc           ! stiffness and compliance
@@ -141,7 +142,7 @@ program DAMASK_spectral_AL
 !--------------------------------------------------------------------------------------------------
 ! loop variables, convergence etc.
  real(pReal) :: time = 0.0_pReal, time0 = 0.0_pReal, timeinc = 1.0_pReal, timeinc_old = 0.0_pReal   ! elapsed time, begin of interval, time interval 
- real(pReal) :: guessmode, err_stress, err_stress_tol, err_f, err_p, err_crit        
+ real(pReal) :: guessmode, err_stress, err_stress_tol, err_f, err_p, err_crit, err_f_point, pstress_av_L2, err_div_rms, err_div     
  real(pReal), dimension(3,3), parameter ::  ones = 1.0_pReal, zeroes = 0.0_pReal
  complex(pReal), dimension(3,3) ::          temp33_Complex
  real(pReal),    dimension(3,3) ::          temp33_Real
@@ -163,8 +164,8 @@ program DAMASK_spectral_AL
  call DAMASK_interface_init
 
  print '(a)', ''
- print '(a)', ' <<<+-  DAMASK_spectral init  -+>>>'
+ print '(a)', ' <<<+-  DAMASK_spectral_AL init  -+>>>'
- print '(a)', ' $Id: DAMASK_spectral.f90 1321 2012-02-15 18:58:38Z MPIE\u.diehl $'
+ print '(a)', ' $Id$'
 #include "compilation_info.f90"
  print '(a,a)', ' Working Directory:    ',trim(getSolverWorkingDirectoryName())
  print '(a,a)', ' Solver Job Name:      ',trim(getSolverJobName())
@@ -344,8 +345,8 @@ program DAMASK_spectral_AL
 ! output of geometry
  print '(a)',  ''
  print '(a)',   '#############################################################'
- print '(a)',   'DAMASK spectral:'
+ print '(a)',   'DAMASK spectral_AL:'
- print '(a)',   'The spectral method boundary value problem solver for'
+ print '(a)',   'The AL spectral method boundary value problem solver for'
  print '(a)',   'the Duesseldorf Advanced Material Simulation Kit'
  print '(a)',   '#############################################################'
  print '(a,a)', 'geometry file:        ',trim(path)//'.geom'
@@ -506,9 +507,8 @@ program DAMASK_spectral_AL
    C = C + dPdF
  enddo; enddo; enddo
  C_inc0 = C * wgt                                                                     ! linear reference material stiffness
- P_av = 0.0_pReal
 
- !--------------------------------------------------------------------------------------------------
+!--------------------------------------------------------------------------------------------------
 ! possible restore deformation gradient from saved state
  if (restartInc > 1_pInt) then                                                                      ! using old values from file                                                      
    if (debugRestart) print '(a,i6,a)' , 'Reading values of increment ',&
@@ -603,8 +603,7 @@ program DAMASK_spectral_AL
 
 !--------------------------------------------------------------------------------------------------
 ! coordinates at beginning of inc
-       call deformed_fft(res,geomdim,math_rotate_backward33(F_aim,bc(loadcase)%rotation),&     ! calculate current coordinates
+       !call deformed_fft(res,geomdim,1.0_pReal,F_real(1:res(1),1:res(2),1:res(3),1:3,1:3),coordinates)! calculate current coordinates
-                             1.0_pReal,F_real(1:res(1),1:res(2),1:res(3),1:3,1:3),coordinates)
 
 !--------------------------------------------------------------------------------------------------
 ! winding forward of deformation aim in loadcase system
@@ -614,6 +613,15 @@ program DAMASK_spectral_AL
                   + deltaF
        F_aim_lastInc = temp33_Real
        F_star_av  = F_aim
+
+!--------------------------------------------------------------------------------------------------
+! Initialize / Update lambda to useful value
+       temp33_real = math_mul3333xx33(C*wgt, F_aim-F_aim_lastInc)
+       P_av = P_av + temp33_real        
+       do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
+          lambda(i,j,k,1:3,1:3) = lambda(i,j,k,1:3,1:3) + temp33_real
+       enddo; enddo; enddo
+
 !--------------------------------------------------------------------------------------------------
 ! update local deformation gradient
        deltaF = math_rotate_backward33(deltaF,bc(loadcase)%rotation)
@@ -630,11 +638,9 @@ program DAMASK_spectral_AL
 !--------------------------------------------------------------------------------------------------
 !Initialize pointwise data for AL scheme: ToDo: good choice?
        F_star(1:res(1),1:res(2),1:res(3),1:3,1:3) = F_real(1:res(1),1:res(2),1:res(3),1:3,1:3)
-       lambda = 0.0_pReal
 
 !--------------------------------------------------------------------------------------------------
 ! calculate reduced compliance
-
        if(size_reduced > 0_pInt) then                                                               ! calculate compliance in case stress BC is applied
          C_lastInc = math_rotate_forward3333(C*wgt,bc(loadcase)%rotation)                      ! calculate stiffness from former inc
          c_prev99 = math_Plain3333to99(C_lastInc)
@@ -677,7 +683,7 @@ program DAMASK_spectral_AL
 !##################################################################################################
 ! convergence loop (looping over iterations)
 !##################################################################################################
-       do while((iter < itmax .and. (err_crit > err_div_tol .or. err_stress > err_stress_tol))&
+       do while((iter < itmax .and. (err_div > err_div_tol .or. err_stress > err_stress_tol .or. err_crit > 5.0e-4))&
                 .or. iter < itmin)
          iter = iter + 1_pInt
 
@@ -700,10 +706,8 @@ program DAMASK_spectral_AL
            err_stress_tol = + huge(1.0_pReal)
          endif
          F_aim_lab = math_rotate_backward33(F_aim,bc(loadcase)%rotation)
-         write (*,'(a,/,3(3(es14.7,1x)/))',advance='no') 'F aim  =',&
+         write (*,'(a,/,3(3(f12.7,1x)/))',advance='no') 'F aim  =',&
                                                 math_transpose33(F_aim)
-         write (*,'(a,/,3(3(es14.7,1x)/))',advance='no') 'F* =',&
-                                                math_transpose33(F_star_av)
 
 !--------------------------------------------------------------------------------------------------
 ! doing Fourier transform
@@ -717,7 +721,32 @@ program DAMASK_spectral_AL
          if(res(3)>1_pInt) &
           lambda_fourier(1:res1_red,1:res(2),                res(3)/2_pInt+1_pInt,1:3,1:3)&
                                                              = cmplx(0.0_pReal,0.0_pReal,pReal)
+!--------------------------------------------------------------------------------------------------
+! calculating RMS divergence criterion in Fourier space
+         pstress_av_L2 = sqrt(maxval(math_eigenvalues33(math_mul33x33(lambda_av,&                    ! L_2 norm of average stress (http://mathworld.wolfram.com/SpectralNorm.html)
+                                                     math_transpose33(lambda_av)))))
+         err_div_RMS = 0.0_pReal
+         do k = 1_pInt, res(3); do j = 1_pInt, res(2)
+           do i = 2_pInt, res1_red -1_pInt                                                          ! Has somewhere a conj. complex counterpart. Therefore count it twice.
+             err_div_RMS = err_div_RMS &
+                   + 2.0_pReal*(sum (real(math_mul33x3_complex(lambda_fourier(i,j,k,1:3,1:3),&           ! (sqrt(real(a)**2 + aimag(a)**2))**2 = real(a)**2 + aimag(a)**2. do not take square root and square again
+                                                   xi(1:3,i,j,k))*TWOPIIMG)**2.0_pReal)&            ! --> sum squared L_2 norm of vector 
+                               +sum(aimag(math_mul33x3_complex(lambda_fourier(i,j,k,1:3,1:3),& 
+                                                                  xi(1:3,i,j,k))*TWOPIIMG)**2.0_pReal))
+           enddo
+           err_div_RMS = err_div_RMS &                                                              ! Those two layers (DC and Nyquist) do not have a conjugate complex counterpart
+                         + sum( real(math_mul33x3_complex(lambda_fourier(1       ,j,k,1:3,1:3),&
+                                                             xi(1:3,1       ,j,k))*TWOPIIMG)**2.0_pReal)&
+                         + sum(aimag(math_mul33x3_complex(lambda_fourier(1       ,j,k,1:3,1:3),&
+                                                             xi(1:3,1       ,j,k))*TWOPIIMG)**2.0_pReal)&
+                         + sum( real(math_mul33x3_complex(lambda_fourier(res1_red,j,k,1:3,1:3),&
+                                                             xi(1:3,res1_red,j,k))*TWOPIIMG)**2.0_pReal)&
+                         + sum(aimag(math_mul33x3_complex(lambda_fourier(res1_red,j,k,1:3,1:3),&
+                                                             xi(1:3,res1_red,j,k))*TWOPIIMG)**2.0_pReal)
+         enddo; enddo
 
+         err_div_RMS = sqrt(err_div_RMS)*wgt
+         err_div = err_div_RMS/pstress_av_L2
 !--------------------------------------------------------------------------------------------------
 ! using gamma operator to update F 
          if(memory_efficient) then                                                                  ! memory saving version, on-the-fly calculation of gamma_hat
@@ -745,36 +774,31 @@ program DAMASK_spectral_AL
            enddo; enddo; enddo
          endif      
          F_fourier(1,1,1,1:3,1:3) = cmplx((F_aim_lab - F_star_av)*real(Npoints,pReal),0.0_pReal,pReal)
+
 !--------------------------------------------------------------------------------------------------
 ! doing inverse Fourier transform
          call fftw_execute_dft_c2r(plan_correction,F_fourier,F_real)                                ! back transform of fluct deformation gradient
-         ! do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
-         ! write (*,'(a,/,3(3(es14.7,1x)/))',advance='no') 'delta F real =',&
-                                                ! math_transpose33(F_real(i,j,k,1:3,1:3)*wgt)
-         ! enddo; enddo; enddo
-
          F_real(1:res(1),1:res(2),1:res(3),1:3,1:3) = F_real(1:res(1),1:res(2),1:res(3),1:3,1:3) * wgt + &
                                                       F_star(1:res(1),1:res(2),1:res(3),1:3,1:3)
 
 !--------------------------------------------------------------------------------------------------
 !
-         print '(a)', '... update stress field P(F*) .....................................'
+         print '(a)', '... update stress field P(F*) and update F* and λ..........................'
          ielem = 0_pInt
-         temp33_Real = 0.0_pReal
          do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
            ielem = ielem + 1_pInt
            call CPFEM_general(3_pInt,&                                                              ! collect cycle
                               coordinates(i,j,k,1:3), F_lastInc(i,j,k,1:3,1:3),&
                               F_star(i,j,k,1:3,1:3),temperature(i,j,k),timeinc,ielem,1_pInt,&
                               sigma,dsde, P, dPdF)
-           temp33_Real = temp33_Real + F_real(i,j,k,1:3,1:3)
          enddo; enddo; enddo
 
-         write (*,'(a,/,3(3(es14.7,1x)/))',advance='no') 'F =',&
-                                                math_transpose33(temp33_Real*wgt)
          ielem = 0_pInt
          err_f = 0.0_pReal
+         err_f_point = 0.0_pReal
          F_star_av = 0.0_pReal
+         P_star_av = 0.0_pReal
+         lambda_av = 0.0_pReal
          do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
            ielem = ielem + 1_pInt
            call CPFEM_general(CPFEM_mode,&
@@ -782,57 +806,72 @@ program DAMASK_spectral_AL
                               F_star(i,j,k,1:3,1:3),temperature(i,j,k),timeinc,ielem,1_pInt,&       
                               sigma,dsde, P,dPdF)
            CPFEM_mode = 2_pInt                                                                  ! winding forward
-
-           if (iter == 1_pInt) lambda(i,j,k,1:3,1:3) = P
            temp33_Real = lambda(i,j,k,1:3,1:3) - P &
                          + math_mul3333xx33(C_inc0,F_real(i,j,k,1:3,1:3)- F_star(i,j,k,1:3,1:3))
-           ! write (*,'(a,/,3(3(es14.7,1x)/))',advance='no') 'F - F* =',&
+
-                                    ! math_transpose33(F_real(i,j,k,1:3,1:3)- F_star(i,j,k,1:3,1:3))
            F_star(i,j,k,1:3,1:3) =  F_star(i,j,k,1:3,1:3) +  &
                                     math_mul3333xx33(math_invSym3333(C_inc0 + dPdF), temp33_Real)
            lambda(i,j,k,1:3,1:3) = lambda(i,j,k,1:3,1:3) + math_mul3333xx33(C_inc0,F_real(i,j,k,1:3,1:3) &
                                                                                  - F_star(i,j,k,1:3,1:3))
            F_star_av = F_star_av + F_star(i,j,k,1:3,1:3)
+           lambda_av = lambda_av + lambda(i,j,k,1:3,1:3)
+           P_star_av = P_star_av + P
            temp33_real = F_star(i,j,k,1:3,1:3) - F_real(i,j,k,1:3,1:3)
+           err_f_point = max(err_f_point, maxval(temp33_real))
            err_f = max(err_f, sqrt(math_mul33xx33(temp33_real,temp33_real)))
-         enddo; enddo; enddo
            
-         F_star_av = F_star_av *wgt
-         write (*,'(a,/,3(3(es14.7,1x)/))',advance='no') 'F* =',&
-                                              math_transpose33(F_star_av)
-
-         print '(a)', '... update stress field P(F) .....................................'
-         ielem = 0_pInt
-         do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
-           ielem = ielem + 1_pInt
-           call CPFEM_general(3_pInt,&                                                              ! collect cycle
-                              coordinates(i,j,k,1:3), F_lastInc(i,j,k,1:3,1:3),&
-                              F_real(i,j,k,1:3,1:3),temperature(i,j,k),timeinc,ielem,1_pInt,&
-                              sigma,dsde,P,dPdF)
-         enddo; enddo; enddo
-         
-         ielem = 0_pInt
-         err_p = 0.0_pReal
-         P_av =0.0_pReal
-         do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
-           ielem = ielem + 1_pInt
-           call CPFEM_general(2_pInt,&
-                              coordinates(i,j,k,1:3),F_lastInc(i,j,k,1:3,1:3), &
-                              F_real(i,j,k,1:3,1:3),temperature(i,j,k),timeinc,ielem,1_pInt,&       
-                              sigma,dsde,P,dPdF)
-           P_av = P_av + P
            temp33_real = lambda(i,j,k,1:3,1:3) - P
            err_p = max(err_p, sqrt(math_mul33xx33(temp33_real,temp33_real)))
          enddo; enddo; enddo
 
-         P_av = math_rotate_forward33(P_av * wgt,bc(loadcase)%rotation)
+         F_star_av = F_star_av *wgt
-         write (*,'(a,/,3(3(es14.7,1x)/))',advance='no') 'P(F) =',&
+         write (*,'(a,/,3(3(f12.7,1x)/))',advance='no') 'F* =',&
-                                              math_transpose33(P_av)/1.e6_pReal
+                                              math_transpose33(F_star_av)
+         P_star_av = P_star_av *wgt
+         write (*,'(a,/,3(3(es14.7,1x)/))',advance='no') 'P(F*) / GPa =',&
+                                              math_transpose33(P_star_av) /1.e6_pReal
+         lambda_av = lambda_av *wgt
+         write (*,'(a,/,3(3(es14.7,1x)/))',advance='no') 'λ / GPa =',&
+                                              math_transpose33(lambda_av) /1.e6_pReal
+         ! print '(a)', '... update stress field P(F) .....................................'
+         ! ielem = 0_pInt
+         ! do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
+           ! ielem = ielem + 1_pInt
+           ! call CPFEM_general(3_pInt,&                                                              ! collect cycle
+                              ! coordinates(i,j,k,1:3), F_lastInc(i,j,k,1:3,1:3),&
+                              ! F_real(i,j,k,1:3,1:3),temperature(i,j,k),timeinc,ielem,1_pInt,&
+                              ! sigma,dsde,P,dPdF)
+         ! enddo; enddo; enddo
+         
+         ! ielem = 0_pInt
+         ! err_p = 0.0_pReal
+         ! P_av =0.0_pReal
+         ! do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
+           ! ielem = ielem + 1_pInt
+           ! call CPFEM_general(2_pInt,&
+                              ! coordinates(i,j,k,1:3),F_lastInc(i,j,k,1:3,1:3), &
+                              ! F_real(i,j,k,1:3,1:3),temperature(i,j,k),timeinc,ielem,1_pInt,&       
+                              ! sigma,dsde,P,dPdF)
+           ! P_av = P_av + P
+           ! temp33_real = lambda(i,j,k,1:3,1:3) - P
+           ! err_p = max(err_p, sqrt(math_mul33xx33(temp33_real,temp33_real)))
+         ! enddo; enddo; enddo
+         
+         ! P_av = math_rotate_forward33(P_av * wgt,bc(loadcase)%rotation)
+         ! write (*,'(a,/,3(3(es14.7,1x)/))',advance='no') 'P(F) / GPa =',&
+                                              ! math_transpose33(P_av)/1.e6_pReal
+         ! write (*,'(a,/,3(3(es14.7,1x)/))',advance='no') 'P(F*) - P(F) =',&
+                                              ! math_transpose33(P_star_av - P_av)
+         P_av = lambda_av
                                               
          err_f = err_f/sqrt(math_mul33xx33(F_star_av,F_star_av))
          err_p = err_p/sqrt(math_mul33xx33(P_av,P_av))
-         write(6,'(a,es14.7,1x)') 'error f', err_f
+         write(6,'(a,es14.7,1x)') 'error F', err_f
-         write(6,'(a,es14.7,1x)') 'error p', err_p
+         write(6,'(a,es14.7,1x)') 'max abs err F', err_f_point
+         write(6,'(a,es14.7,1x)') 'error P', err_p
+         write(6,'(a,es11.4)')        'error divergence  FT  RMS = ',err_div_RMS
+         write(6,'(a,es11.4)')        'error div = ',err_div
+         write(6,'(a,es11.4)')        'error stress = ',err_stress/err_stress_tol
          
        err_crit = max(err_p, err_f)
        
@@ -883,11 +922,16 @@ end program DAMASK_spectral_AL
 !
 !********************************************************************
 subroutine quit(stop_id)
- use prec
+ use prec, only: &
- implicit none
+   pInt
    
+ implicit none
  integer(pInt), intent(in) :: stop_id
  
- ! if (stop_id ==    0_pInt) stop 0_pInt                                                              ! normal termination
+ if (stop_id ==    0_pInt) stop 0                                                                   ! normal termination
+ if (stop_id <= 9000_pInt) then                                                                     ! trigger regridding
+   write(6,'(i4)') stop_id
+   stop 1
+ endif
  stop 'abnormal termination of DAMASK_spectral'
 end subroutine

code/math.f90

@@ -35,47 +35,52 @@ module math
  complex(pReal), parameter, public :: TWOPIIMG = (0.0_pReal,2.0_pReal)* pi
 
 ! *** 3x3 Identity ***
- real(pReal), dimension(3,3), parameter, public :: math_I3 = &
+ real(pReal), dimension(3,3), parameter, public :: &
- reshape( (/ &
+   math_I3 = reshape([&
      1.0_pReal,0.0_pReal,0.0_pReal, &
      0.0_pReal,1.0_pReal,0.0_pReal, &
- 0.0_pReal,0.0_pReal,1.0_pReal /),(/3,3/))
+     0.0_pReal,0.0_pReal,1.0_pReal  &
+     ],[3,3])
 
 ! *** Mandel notation ***
- integer(pInt), dimension (2,6), parameter :: mapMandel = &
+ integer(pInt), dimension (2,6), parameter, private :: &
- reshape((/&
+   mapMandel = reshape([&
      1_pInt,1_pInt, &
      2_pInt,2_pInt, &
      3_pInt,3_pInt, &
      1_pInt,2_pInt, &
      2_pInt,3_pInt, &
      1_pInt,3_pInt  &
- /),(/2,6/))
+     ],[2,6])
 
- real(pReal), dimension(6), parameter :: nrmMandel = &
+ real(pReal), dimension(6), parameter, private :: &
- (/1.0_pReal,1.0_pReal,1.0_pReal, 1.414213562373095_pReal, 1.414213562373095_pReal, 1.414213562373095_pReal/)
+   nrmMandel = [&
- real(pReal), dimension(6), parameter :: invnrmMandel = &
+     1.0_pReal,                1.0_pReal,                1.0_pReal,&
- (/1.0_pReal,1.0_pReal,1.0_pReal,0.7071067811865476_pReal,0.7071067811865476_pReal,0.7071067811865476_pReal/)
+     1.414213562373095_pReal,  1.414213562373095_pReal,  1.414213562373095_pReal]
+     
+ real(pReal), dimension(6), parameter , public :: &
+   invnrmMandel = [&
+     1.0_pReal,                1.0_pReal,                1.0_pReal,&
+     0.7071067811865476_pReal, 0.7071067811865476_pReal, 0.7071067811865476_pReal]
 
 ! *** Voigt notation ***
- integer(pInt), dimension (2,6), parameter :: mapVoigt = &
+ integer(pInt), dimension (2,6), parameter, private :: &
- reshape((/&
+   mapVoigt = reshape([&
      1_pInt,1_pInt, &
      2_pInt,2_pInt, &
      3_pInt,3_pInt, &
      2_pInt,3_pInt, &
      1_pInt,3_pInt, &
      1_pInt,2_pInt  &
- /),(/2,6/))
+     ],[2,6])
 
- real(pReal), dimension(6), parameter :: nrmVoigt = &
+ real(pReal), dimension(6), parameter, private :: & 
- (/1.0_pReal,1.0_pReal,1.0_pReal, 1.0_pReal, 1.0_pReal, 1.0_pReal/)
+   nrmVoigt    = 1.0_pReal, &
- real(pReal), dimension(6), parameter :: invnrmVoigt = &
+   invnrmVoigt = 1.0_pReal
- (/1.0_pReal,1.0_pReal,1.0_pReal, 1.0_pReal, 1.0_pReal, 1.0_pReal/)
 
 ! *** Plain notation ***
- integer(pInt), dimension (2,9), parameter :: mapPlain = &
+ integer(pInt), dimension (2,9), parameter, private :: &
- reshape((/&
+   mapPlain = reshape([&
      1_pInt,1_pInt, &
      1_pInt,2_pInt, &
      1_pInt,3_pInt, &
@@ -85,13 +90,15 @@ module math
      3_pInt,1_pInt, &
      3_pInt,2_pInt, &
      3_pInt,3_pInt  &
- /),(/2,9/))
+     ],[2,9])
 
 ! Symmetry operations as quaternions
 ! 24 for cubic, 12 for hexagonal = 36
-integer(pInt), dimension(2), parameter :: math_NsymOperations = (/24_pInt,12_pInt/)
+integer(pInt), dimension(2), parameter, private :: &
-real(pReal), dimension(4,36), parameter :: math_symOperations = &
+   math_NsymOperations = [24_pInt,12_pInt]
-  reshape((/&
+   
+real(pReal), dimension(4,36), parameter, private :: &
+  math_symOperations = reshape([&
      1.0_pReal,                 0.0_pReal,                 0.0_pReal,                 0.0_pReal, &                      ! cubic symmetry operations
      0.0_pReal,                 0.0_pReal,                 0.7071067811865476_pReal,  0.7071067811865476_pReal, &       !     2-fold symmetry
      0.0_pReal,                 0.7071067811865476_pReal,  0.0_pReal,                 0.7071067811865476_pReal, &
@@ -128,12 +135,20 @@ real(pReal), dimension(4,36), parameter :: math_symOperations = &
      0.5_pReal,                 0.0_pReal,                 0.0_pReal,                 0.866025403784439_pReal, &
     -0.5_pReal,                 0.0_pReal,                 0.0_pReal,                 0.866025403784439_pReal, &
      0.0_pReal,                 0.0_pReal,                 0.0_pReal,                 1.0_pReal &
-  /),(/4,36/))
+     ],[4,36])
   
  include 'fftw3.f03'
  
  public  :: math_init, &
-            math_range
+            qsort, &
+            math_range, &
+            math_identity2nd, &
+            math_civita
+            
+ private :: math_partition, &
+            math_delta, &
+            Gauss
+ 
 contains
  
 !**************************************************************************
@@ -308,7 +323,6 @@ end function math_partition
 pure function math_range(N)  
 
  implicit none
-
  integer(pInt), intent(in) :: N
  integer(pInt) :: i
  integer(pInt), dimension(N) :: math_range
@@ -324,7 +338,6 @@ end function math_range
 pure function math_identity2nd(dimen)  
 
  implicit none
-
  integer(pInt), intent(in) :: dimen
  integer(pInt) :: i
  real(pReal), dimension(dimen,dimen) :: math_identity2nd
@@ -344,7 +357,6 @@ end function math_identity2nd
 pure function math_civita(i,j,k)
 
  implicit none
-
  integer(pInt), intent(in) :: i,j,k
  real(pReal) math_civita
 
@@ -367,7 +379,6 @@ end function math_civita
 pure function math_delta(i,j)
 
  implicit none
-
  integer(pInt), intent (in) :: i,j
  real(pReal) :: math_delta
 
@@ -383,7 +394,6 @@ end function math_delta
 pure function math_identity4th(dimen)  
 
  implicit none
-
  integer(pInt), intent(in) :: dimen
  integer(pInt) :: i,j,k,l
  real(pReal), dimension(dimen,dimen,dimen,dimen) ::  math_identity4th
@@ -400,7 +410,6 @@ end function math_identity4th
 pure function math_vectorproduct(A,B)  
 
  implicit none
-
  real(pReal), dimension(3), intent(in) ::  A,B
  real(pReal), dimension(3) ::  math_vectorproduct
 
@@ -505,7 +514,6 @@ end function math_mul3333xx33
 pure function math_mul3333xx3333(A,B)  
 
  implicit none
-
  integer(pInt) :: i,j,k,l
  real(pReal), dimension(3,3,3,3), intent(in) ::  A
  real(pReal), dimension(3,3,3,3), intent(in) ::  B
@@ -527,7 +535,6 @@ end function math_mul3333xx3333
 pure function math_mul33x33(A,B)  
 
  implicit none
-
  integer(pInt) :: i,j
  real(pReal), dimension(3,3), intent(in) ::  A,B
  real(pReal), dimension(3,3) ::  math_mul33x33
@@ -544,7 +551,6 @@ end function math_mul33x33
 pure function math_mul66x66(A,B)  
 
  implicit none
-
  integer(pInt) :: i,j
  real(pReal), dimension(6,6), intent(in) ::  A,B
  real(pReal), dimension(6,6) ::  math_mul66x66
@@ -562,8 +568,8 @@ end function math_mul66x66
 pure function math_mul99x99(A,B)  
 
  use prec, only: pReal, pInt
- implicit none
 
+ implicit none
  integer(pInt)  i,j
  real(pReal), dimension(9,9), intent(in) ::  A,B
 
@@ -584,7 +590,6 @@ end function math_mul99x99
 pure function math_mul33x3(A,B)  
 
  implicit none
-
  integer(pInt) :: i
  real(pReal), dimension(3,3), intent(in) ::  A
  real(pReal), dimension(3),   intent(in) ::  B
@@ -600,7 +605,6 @@ end function math_mul33x3
 pure function math_mul33x3_complex(A,B)  
 
  implicit none
-
  integer(pInt) :: i
  complex(pReal), dimension(3,3), intent(in) ::  A
  real(pReal),    dimension(3),   intent(in) ::  B
@@ -636,7 +640,6 @@ end function math_mul66x6
 function math_qRnd()  
 
  implicit none
-
  real(pReal), dimension(4) :: math_qRnd
  real(pReal), dimension(3) :: rnd
  
@@ -655,7 +658,6 @@ end function math_qRnd
 pure function math_qMul(A,B)  
 
  implicit none
-
  real(pReal), dimension(4), intent(in) ::  A, B
  real(pReal), dimension(4) ::  math_qMul
 
@@ -673,7 +675,6 @@ end function math_qMul
 pure function math_qDot(A,B)  
 
  implicit none
-
  real(pReal), dimension(4), intent(in) :: A, B
  real(pReal) :: math_qDot
 
@@ -688,7 +689,6 @@ end function math_qDot
 pure function math_qConj(Q)  
 
  implicit none
-
  real(pReal), dimension(4), intent(in) ::  Q
  real(pReal), dimension(4) ::  math_qConj
 
@@ -704,7 +704,6 @@ end function math_qConj
 pure function math_qNorm(Q)  
 
  implicit none
-
  real(pReal), dimension(4), intent(in) ::  Q
  real(pReal) :: math_qNorm
  
@@ -719,7 +718,6 @@ end function math_qNorm
 pure function math_qInv(Q)  
 
  implicit none
-
  real(pReal), dimension(4), intent(in) ::  Q
  real(pReal), dimension(4) ::  math_qInv
  real(pReal) :: squareNorm
@@ -739,7 +737,6 @@ end function math_qInv
 pure function math_qRot(Q,v)  
 
  implicit none
-
  real(pReal), dimension(4), intent(in) :: Q
  real(pReal), dimension(3), intent(in) :: v
  real(pReal), dimension(3) :: math_qRot
@@ -767,7 +764,6 @@ end function math_qRot
 pure function math_transpose33(A)
 
  implicit none
-
  real(pReal),dimension(3,3),intent(in)  :: A
  real(pReal),dimension(3,3) :: math_transpose33
  integer(pInt) :: i,j
@@ -3098,7 +3094,7 @@ end subroutine shape_compare
 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 
 subroutine mesh_regular_grid(res,geomdim,defgrad_av,centroids,nodes)
 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 
-! Routine to build mesh of (distoreted) cubes for given coordinates (= center of the cubes)
+! Routine to build mesh of (distorted) cubes for given coordinates (= center of the cubes)
 !
  use debug, only: debug_math, &
                   debug_what, &
@@ -3336,7 +3332,6 @@ subroutine deformed_fft(res,geomdim,defgrad_av,scaling,defgrad,coords)
  coords_fftw =          fftw_alloc_complex(int(res1_red     *res(2)*res(3)*3_pInt,C_SIZE_T))        !C_SIZE_T is of type integer(8)
  call c_f_pointer(coords_fftw, coords_real,     [res(1)+2_pInt,res(2),res(3),3_pInt])
  call c_f_pointer(coords_fftw, coords_fourier,  [res1_red     ,res(2),res(3),3_pInt])
-
  fftw_forth = fftw_plan_many_dft_r2c(3_pInt,(/res(3),res(2) ,res(1)/),9_pInt,&                      ! dimensions , length in each dimension in reversed order
                           defgrad_real,(/res(3),res(2) ,res(1)+2_pInt/),&                               ! input data , physical length in each dimension in reversed order
                                      1_pInt,  res(3)*res(2)*(res(1)+2_pInt),&                                ! striding   , product of physical lenght in the 3 dimensions
@@ -3368,7 +3363,6 @@ subroutine deformed_fft(res,geomdim,defgrad_av,scaling,defgrad,coords)
                                            1:3,1:3) = cmplx(0.0_pReal,0.0_pReal,pReal)
                                                
  coords_fourier = cmplx(0.0_pReal,0.0_pReal,pReal)
- 
  do k = 1_pInt, res(3)
    k_s(3) = k-1_pInt
    if(k > res(3)/2_pInt+1_pInt) k_s(3) = k_s(3)-res(3)

processing/pre/OIMang_hex2cub.py

@@ -0,0 +1,108 @@
+#!/usr/bin/env python
+# -*- coding: UTF-8 no BOM -*-
+
+import string,os,sys
+from optparse import OptionParser, Option
+
+# -----------------------------
+class extendableOption(Option):
+# -----------------------------
+# used for definition of new option parser action 'extend', which enables to take multiple option arguments
+# taken from online tutorial http://docs.python.org/library/optparse.html
+  
+  ACTIONS = Option.ACTIONS + ("extend",)
+  STORE_ACTIONS = Option.STORE_ACTIONS + ("extend",)
+  TYPED_ACTIONS = Option.TYPED_ACTIONS + ("extend",)
+  ALWAYS_TYPED_ACTIONS = Option.ALWAYS_TYPED_ACTIONS + ("extend",)
+
+  def take_action(self, action, dest, opt, value, values, parser):
+    if action == "extend":
+      lvalue = value.split(",")
+      values.ensure_value(dest, []).extend(lvalue)
+    else:
+      Option.take_action(self, action, dest, opt, value, values, parser)
+
+ 
+# --------------------------------------------------------------------
+#                                MAIN
+# --------------------------------------------------------------------
+
+parser = OptionParser(option_class=extendableOption, usage='%prog options [file[s]]', description = """
+Converts ang files (EBSD Data) from hexagonal grid to a pixel grid
+
+""" + string.replace('$Id$','\n','\\n')
+)
+
+parser.add_option('-x', dest='columnX', action='store', type='int', \
+                        help='column containing x coordinates [%default]')
+
+parser.set_defaults(columnX = 3)
+
+(options,filenames) = parser.parse_args()
+
+counterX  = 0
+counterY  = 0
+addPoints = -1      # No of doubled points (must be the same for each odd/even line, initializing with -1 make countin easy!)
+
+# ------------------------------------------ setup file handles ---------------------------------------  
+
+files = []
+if filenames == []:
+  files.append({'name':'STDIN', 'input':sys.stdin, 'output':sys.stdout})
+else:
+  for name in filenames:
+    if os.path.exists(name):
+      files.append(   {'name':name, 'input':open(name),'output':open(os.path.splitext(name)[0]\
+                                                             +'_cub'+os.path.splitext(name)[1], 'w')})
+
+# ------------------------------------------ loop over input files ---------------------------------------  
+
+for file in files:
+  print file['name']
+  x = 0
+  for line in file['input']:
+    lineSplit=line.split()
+    
+    if lineSplit[0]=='#':
+      if len(lineSplit)>2:                                                                          # possibly interesting information
+        if line.split()[2]=='SqrGrid': 
+          print 'The file is already a square grid file.'
+          sys.exit()
+        if lineSplit[1]=='XSTEP:':      stepSizeX = float(lineSplit[2]) 
+        if lineSplit[1]=='YSTEP:':      stepSizeY = float(lineSplit[2])
+        
+        if lineSplit[2]=='HexGrid':     line='# GRID: SqrGrid\n'                                    # comments are not read by OIM, but still better to be correct
+        if lineSplit[1]=='NCOLS_ODD:':
+          NCols = int(int(lineSplit[2])*stepSizeX/stepSizeY)
+          line='# NCOLS_ODD: %d\n'% NCols
+        if lineSplit[1]=='NCOLS_EVEN:':
+          line='# NCOLS_EVEN: %d\n'% NCols
+          
+      file['output'].write(line)
+    else:                                                                                           # finished reading of header      
+      xOld = x
+      x = float(lineSplit[options.columnX])                                                         # current (original) x positions
+      
+      if x >  xOld:                                                                                 # same line, increase X
+        counterX+=1
+      else:                                                                                         # new line, increase in Y, reset X
+        counterY+=1
+        addPoints = -1                                                                              # to start at zero
+        counterX=0
+      
+      lineFirstPart =''                                                                             # split line around x and y coordinate
+      for i in xrange(options.columnX):
+        lineFirstPart =lineFirstPart+' '+lineSplit[i]
+      lineLastPart =''
+      for i in xrange(len(lineSplit)- (options.columnX+2)):
+        lineLastPart =lineLastPart+' '+lineSplit[i+options.columnX+2]
+        
+      if counterX+addPoints < NCols:
+        file['output'].write('%s %.6f %.6f  %s\n' %(lineFirstPart,(counterX+addPoints)*stepSizeY,   # write with new x and y position 
+                                                                  counterY*stepSizeY,lineLastPart))                     
+
+      if x - (counterX+addPoints)*stepSizeY > 0.5*stepSizeY and counterX+addPoints+1 < NCols:       # double point (interpolation error)
+                   
+        addPoints+=1
+        file['output'].write('%s %.6f %.6f  %s\n' %(lineFirstPart,(counterX+addPoints)*stepSizeY,\
+                                                                   counterY*stepSizeY,lineLastPart))

processing/pre/ang_hex2cub.py

@@ -1,84 +0,0 @@
-#!/usr/bin/env python
-# -*- coding: UTF-8 no BOM -*-
-
-import string,os,sys
-from optparse import OptionParser, OptionGroup, Option, SUPPRESS_HELP
-
-# -----------------------------
-class extendedOption(Option):
-# -----------------------------
-# used for definition of new option parser action 'extend', which enables to take multiple option arguments
-# taken from online tutorial http://docs.python.org/library/optparse.html
-    
-    ACTIONS = Option.ACTIONS + ("extend",)
-    STORE_ACTIONS = Option.STORE_ACTIONS + ("extend",)
-    TYPED_ACTIONS = Option.TYPED_ACTIONS + ("extend",)
-    ALWAYS_TYPED_ACTIONS = Option.ALWAYS_TYPED_ACTIONS + ("extend",)
-
-    def take_action(self, action, dest, opt, value, values, parser):
-        if action == "extend":
-            lvalue = value.split(",")
-            values.ensure_value(dest, []).extend(lvalue)
-        else:
-            Option.take_action(self, action, dest, opt, value, values, parser)
-
-parser = OptionParser(option_class=extendedOption, usage='%prog [geomfile[s]]', description = """
-Converts EBSD Data stored in *.ang files from hex to cub
-
-""" + string.replace('$Id: spectral_geomCheck.py 1084 2011-11-09 15:37:45Z MPIE\c.zambaldi $','\n','\\n')
-)
-
-(options, filenames) = parser.parse_args()
-
-# ------------------------------------------ setup file handles ---------------------------------------  
-
-columnX   = 3 # 0,1,2,3 (python notation!)
-counterX  = 0
-counterY  = 0
-addPoints = 0
-
-files = []
-for name in filenames:
-  if os.path.exists(name):
-    files.append(   {'name':name, 'input':open(name),'output':open('cub_'+name, 'w')})
-
-# ------------------------------------------ loop over input files ---------------------------------------  
-
-for file in files:
-  print file['name']
-  for line in file['input']:
-    if line.split()[0]=='#':
-      if len(line.split())>2:                                                                               # possibly interesting information
-        if line.split()[2]=='SqrGrid': 
-          print 'The file is already a square grid file.'
-          sys.exit()
-        if line.split()[2]=='HexGrid': line='# GRID: SqrGrid\n'
-        if line.split()[1]=='XSTEP:':      stepSizeX = float(line.split()[2]) 
-        if line.split()[1]=='YSTEP:':      stepSizeY = float(line.split()[2]) 
-        if line.split()[1]=='NCOLS_EVEN:': NColsEven = int(line.split()[2]) 
-      file['output'].write(line)
-    else:                                                                                                   # finished reading of header
-      lineSplit=line.split()
-      x = float(lineSplit[columnX])
-      y = float(lineSplit[columnX+1])
-      lineFirstPart =''
-      lineLastPart =''
-      for i in xrange(columnX):
-        lineFirstPart =lineFirstPart+' '+lineSplit[i]
-      for i in xrange(len(lineSplit)- (columnX+2)):
-        lineLastPart =lineLastPart+' '+lineSplit[i+columnX+2]
-
-      file['output'].write(lineFirstPart+' '+\
-                   str((counterX+addPoints)*stepSizeY)+' '+str(y)+' '+\
-                   lineLastPart+'\n')
-      if x + stepSizeX - (counterX+addPoints+1)*stepSizeY > 0.5*stepSizeY:                                 # double point (interpolation error)
-        addPoints+=1
-        file['output'].write(lineFirstPart+' '+\
-                   str((counterX+addPoints)*stepSizeY)+' '+str(y)+' '+\
-                   lineLastPart+'\n')
-
-      if(counterX == NColsEven + counterY%2):                                                              # new row (odd and even differ by 1)                   
-        counterY+=1
-        counterX=0
-        addPoints=0
-      counterX+=1

processing/pre/spectral_ang2geom.py

@@ -0,0 +1,115 @@
+#!/usr/bin/env python
+
+import os,sys,math,string,numpy
+from optparse import OptionParser, Option
+
+# -----------------------------
+class extendableOption(Option):
+# -----------------------------
+# used for definition of new option parser action 'extend', which enables to take multiple option arguments
+# taken from online tutorial http://docs.python.org/library/optparse.html
+  
+  ACTIONS = Option.ACTIONS + ("extend",)
+  STORE_ACTIONS = Option.STORE_ACTIONS + ("extend",)
+  TYPED_ACTIONS = Option.TYPED_ACTIONS + ("extend",)
+  ALWAYS_TYPED_ACTIONS = Option.ALWAYS_TYPED_ACTIONS + ("extend",)
+
+  def take_action(self, action, dest, opt, value, values, parser):
+    if action == "extend":
+      lvalue = value.split(",")
+      values.ensure_value(dest, []).extend(lvalue)
+    else:
+      Option.take_action(self, action, dest, opt, value, values, parser)
+
+
+# --------------------------------------------------------------------
+#                                MAIN
+# --------------------------------------------------------------------
+
+parser = OptionParser(option_class=extendableOption, usage='%prog options [file[s]]', description = """
+Converts EBSD data from cubic ang files into description for spectral solver (*.geom + material.config)
+Can discriminate two phases depending on threshold value 
+
+""" + string.replace('$Id$','\n','\\n')
+)
+
+
+parser.add_option('-c','--column',     dest='column', type='int', \
+                                       help='column containgin value to disciminate phase 1 and 2 [%default]')
+parser.add_option('-t','--threshold',  dest='threshold', type='float', \
+                                       help='threshold to disciminate phases. if (value < treshold) phase = 1')
+parser.add_option('-l','--long',       dest='useNoRange', action='store_true',\
+                                       help='write number for each point instead of "1 to N" in geom file [%default]')
+
+parser.set_defaults(column = 1)
+parser.set_defaults(threshold = sys.maxint)
+parser.set_defaults(useRange = False)
+
+(options,filenames) = parser.parse_args()
+
+# ------------------------------------------ setup file handles ---------------------------------------  
+eulers=numpy.array([0.0,0.0,0.0],'f')
+geomdim=numpy.array([0.0,0.0,0.0],'f')
+res=numpy.array([0,0,1],'i')      
+
+files = []
+if filenames == []:
+  files.append({'name':'STDIN', 'input':sys.stdin, 'material':sys.stdout, 'geom':sys.stdout})
+else:
+  for name in filenames:
+    if os.path.exists(name):
+      files.append(   {'name':name, 'input':open(name),\
+                                    'material':open(os.path.splitext(name)[0]+'_material.config', 'w'),\
+                                    'geom':open(os.path.splitext(name)[0]+'.geom', 'w')})
+
+# ------------------------------------------ loop over input files ---------------------------------------  
+for file in files:
+  point=0
+  if file['name'] != 'STDIN': print file['name']
+  file['material'].write('#---\n<homogenization>\n#---\n'+
+                         '\n[SX]\ntype isostrain\nNgrains 1\n\n'+
+                         '#---\n<microstructure>\n#---\n\n')
+  tempPart2=             '#---\n<texture>\n---\n\n'
+  for line in file['input']:
+    lineSplit=line.split()
+
+    if line[0]=='#':
+      if len(lineSplit)>2:
+        if line.split()[2]=='HexGrid': 
+          print 'The file is a hex grid file. Convert it first to sqr grid'
+          sys.exit()
+        if lineSplit[1]=='XSTEP:':      stepSizeX = float(lineSplit[2]) 
+        if lineSplit[1]=='YSTEP:':      stepSizeY = float(lineSplit[2]) 
+        if lineSplit[1]=='NCOLS_ODD:':  res[0]    = int(lineSplit[2]) 
+        if lineSplit[1]=='NROWS:':      res[1]    = int(lineSplit[2]) 
+    else:
+      point+=1
+      eulers = (float(lineSplit[0])/2.0/math.pi*360.0, \
+                   float(lineSplit[1])/2.0/math.pi*360.0, \
+                   float(lineSplit[2])/2.0/math.pi*360.0)
+               
+      if float(lineSplit[options.column-1])<options.threshold:
+        phase=1
+      else:
+        phase=2
+        
+      file['material'].write(\
+        '[Grain%08d]\ncrystallite 1\n(constituent) phase %1d texture %08d fraction 1.0\n' \
+                  %(point,phase,point))
+      tempPart2+=\
+        '[Grain%08d]\n(gauss) phi1 %4.2f Phi %4.2f phi2 %4.2f scatter 0.0 fraction 1.0\n'\
+                  %(point,eulers[0],eulers[1],eulers[2])
+  geomdim[0] = stepSizeX*res[0]
+  geomdim[1] = stepSizeY*res[1]
+  geomdim[2] = min(stepSizeX,stepSizeY)
+    
+  file['material'].write(tempPart2)
+  file['geom'].write(
+    'resolution a %4d b %4d c %1d \ndimension %5.3f y %5.3f z %5.3f\nhomogenization 1\n'\
+               %(res[0],res[1],res[2],geomdim[0],geomdim[1],geomdim[2]))
+  if options.useNoRange:
+    for x in xrange(res[0]*res[1]):
+      file['geom'].write('%08d' %x)
+  else:
+    file['geom'].write('1 to %d\n'%(res[0]*res[1]))
+

processing/pre/spectral_randomSeeding.py

@@ -39,7 +39,7 @@ mappings = {
 parser = OptionParser(option_class=extendedOption, usage='%prog options [file[s]]', description = """
 Offset microstructure index for points which see a microstructure different from themselves within a given (cubic) vicinity,
 i.e. within the region close to a grain/phase boundary.
-""" + string.replace('$Id: spectral_geomCheck 994 2011-09-05 13:38:10Z MPIE\p.eisenlohr $','\n','\\n')
+""" + string.replace('$Id$','\n','\\n')
 )
 
 parser.add_option('-f', '--file', dest='filename', type="string", \

processing/setup/setup_processing.py

@@ -50,9 +50,9 @@ if options.compiler not in compilers:
 f2py_compiler = {
                   'gfortran': 'gnu95   --f90flags="-fno-range-check -xf95-cpp-input -std=f2008"',
                   'gnu95':    'gnu95   --f90flags="-fno-range-check -xf95-cpp-input -std=f2008"',
-                  'intel32':  'intel   --f90flags="-fpp -stand f03"',
+                  'intel32':  'intel   --f90flags="-fpp -stand f03 -diag-disable 5268"',
-                  'intel':    'intelem --f90flags="-fpp -stand f03"',
+                  'intel':    'intelem --f90flags="-fpp -stand f03 -diag-disable 5268"',
-                  'ifort':    'intelem --f90flags="-fpp -stand f03"',
+                  'ifort':    'intelem --f90flags="-fpp -stand f03 -diag-disable 5268"',
                 }[options.compiler]
 compiler = {
                   'gfortran': 'gfortran',
@@ -85,8 +85,10 @@ bin_link = { \
                 'spectral_geomCrop.py',
                 'spectral_minimalSurface.py',
                 'spectral_vicinityOffset.py',
-                'voronoi_randomSeeding.py',
+                'spectral_ang2geom.py',
+                'spectral_randomSeeding.py',
                 'voronoi_tessellation.exe',
+                'OIMang_hex2cub',
                 ],
         'post' : [
                 '3Dvisualize.py',
@@ -128,11 +130,16 @@ execute = { \
                         # It uses the fortran wrapper f2py that is included in the numpy package to construct the
                         # module postprocessingMath.so out of the fortran code postprocessingMath.f90
                         # for the generation of the pyf file:
-                        #f2py -m DAMASK -h DAMASK.pyf --overwrite-signature ../../code/math.f90 \
+                        # f2py -m DAMASK -h DAMASK.pyf --overwrite-signature ../../code/math.f90 \
                         'f2py %s'%(os.path.join(codeDir,'damask.core.pyf')) +\
                         ' -c --fcompiler=%s'%(f2py_compiler) +\
-                        ' %s'%(os.path.join(codeDir,'core_modules.f90'))+\
+                        ' %s'%(os.path.join(codeDir,'prec.f90'))+\
+                        ' %s'%(os.path.join(codeDir,'DAMASK_python_interface.f90'))+\
+                        ' %s'%(os.path.join(codeDir,'IO.f90'))+\
+                        ' %s'%(os.path.join(codeDir,'numerics.f90'))+\
+                        ' %s'%(os.path.join(codeDir,'debug.f90'))+\
                         ' %s'%(os.path.join(codeDir,'math.f90'))+\
+                        ' %s'%(os.path.join(codeDir,'DAMASK_python.f90'))+\
                         ' -L%s/lib -lfftw3'%(damaskEnv.pathInfo['fftw'])+\
                         ' %s'%lib_lapack,
                         'mv %s `readlink -f %s`' %(os.path.join(codeDir,'core.so'),os.path.join(damaskEnv.relPath('lib/damask'),'core.so')),