diff --git a/code/DAMASK_python.f90 b/code/DAMASK_python.f90
new file mode 100644
index 000000000..1017f7a01
--- /dev/null
+++ b/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 .
+!
+!##################################################################################################
+!* $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
diff --git a/code/DAMASK_python_interface.f90 b/code/DAMASK_python_interface.f90
new file mode 100644
index 000000000..9979c7bed
--- /dev/null
+++ b/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 .
+!
+!--------------------------------------------------------------------------------------------------
+!* $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 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(640) 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(640) 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
diff --git a/code/DAMASK_spectral.f90 b/code/DAMASK_spectral.f90
index 2343b730c..01157c43e 100644
--- a/code/DAMASK_spectral.f90
+++ b/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)
diff --git a/code/DAMASK_spectral_AL.f90 b/code/DAMASK_spectral_AL.f90
index c2608ae73..24677987a 100644
--- a/code/DAMASK_spectral_AL.f90
+++ b/code/DAMASK_spectral_AL.f90
@@ -17,7 +17,7 @@
! along with DAMASK. If not, see .
!
!##################################################################################################
-!* $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)', ' $Id: DAMASK_spectral.f90 1321 2012-02-15 18:58:38Z MPIE\u.diehl $'
+ print '(a)', ' <<<+- DAMASK_spectral_AL init -+>>>'
+ 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)', 'The spectral method boundary value problem solver for'
+ print '(a)', 'DAMASK spectral_AL:'
+ 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
- 1.0_pReal,F_real(1:res(1),1:res(2),1:res(3),1:3,1:3),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
!--------------------------------------------------------------------------------------------------
! 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)
- write (*,'(a,/,3(3(es14.7,1x)/))',advance='no') 'P(F) =',&
- math_transpose33(P_av)/1.e6_pReal
+ F_star_av = F_star_av *wgt
+ write (*,'(a,/,3(3(f12.7,1x)/))',advance='no') 'F* =',&
+ 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 p', err_p
+ write(6,'(a,es14.7,1x)') 'error F', err_f
+ 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: &
+ pInt
+
implicit none
-
integer(pInt), intent(in) :: stop_id
- ! if (stop_id == 0_pInt) stop 0_pInt ! normal termination
- stop 'abnormal termination of DAMASK_spectral'
+ 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
diff --git a/code/math.f90 b/code/math.f90
index c9514c668..dbc4e123c 100644
--- a/code/math.f90
+++ b/code/math.f90
@@ -35,63 +35,70 @@ 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 = &
- 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/))
+ real(pReal), dimension(3,3), parameter, public :: &
+ 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])
! *** Mandel notation ***
- integer(pInt), dimension (2,6), parameter :: 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/))
+ integer(pInt), dimension (2,6), parameter, private :: &
+ 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])
- real(pReal), dimension(6), parameter :: nrmMandel = &
- (/1.0_pReal,1.0_pReal,1.0_pReal, 1.414213562373095_pReal, 1.414213562373095_pReal, 1.414213562373095_pReal/)
- real(pReal), dimension(6), parameter :: invnrmMandel = &
- (/1.0_pReal,1.0_pReal,1.0_pReal,0.7071067811865476_pReal,0.7071067811865476_pReal,0.7071067811865476_pReal/)
+ real(pReal), dimension(6), parameter, private :: &
+ nrmMandel = [&
+ 1.0_pReal, 1.0_pReal, 1.0_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 = &
- 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/))
+ integer(pInt), dimension (2,6), parameter, private :: &
+ 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])
- real(pReal), dimension(6), parameter :: nrmVoigt = &
- (/1.0_pReal,1.0_pReal,1.0_pReal, 1.0_pReal, 1.0_pReal, 1.0_pReal/)
- real(pReal), dimension(6), parameter :: invnrmVoigt = &
- (/1.0_pReal,1.0_pReal,1.0_pReal, 1.0_pReal, 1.0_pReal, 1.0_pReal/)
+ real(pReal), dimension(6), parameter, private :: &
+ nrmVoigt = 1.0_pReal, &
+ invnrmVoigt = 1.0_pReal
! *** Plain notation ***
- integer(pInt), dimension (2,9), parameter :: mapPlain = &
- reshape((/&
- 1_pInt,1_pInt, &
- 1_pInt,2_pInt, &
- 1_pInt,3_pInt, &
- 2_pInt,1_pInt, &
- 2_pInt,2_pInt, &
- 2_pInt,3_pInt, &
- 3_pInt,1_pInt, &
- 3_pInt,2_pInt, &
- 3_pInt,3_pInt &
- /),(/2,9/))
+ integer(pInt), dimension (2,9), parameter, private :: &
+ mapPlain = reshape([&
+ 1_pInt,1_pInt, &
+ 1_pInt,2_pInt, &
+ 1_pInt,3_pInt, &
+ 2_pInt,1_pInt, &
+ 2_pInt,2_pInt, &
+ 2_pInt,3_pInt, &
+ 3_pInt,1_pInt, &
+ 3_pInt,2_pInt, &
+ 3_pInt,3_pInt &
+ ],[2,9])
! Symmetry operations as quaternions
! 24 for cubic, 12 for hexagonal = 36
-integer(pInt), dimension(2), parameter :: math_NsymOperations = (/24_pInt,12_pInt/)
-real(pReal), dimension(4,36), parameter :: math_symOperations = &
- reshape((/&
+integer(pInt), dimension(2), parameter, private :: &
+ math_NsymOperations = [24_pInt,12_pInt]
+
+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, &
@@ -3316,7 +3312,7 @@ subroutine deformed_fft(res,geomdim,defgrad_av,scaling,defgrad,coords)
integer(pInt) :: i, j, k, res1_red
integer(pInt), dimension(3) :: k_s
real(pReal), dimension(3) :: step, offset_coords, integrator
-
+
integrator = geomdim / 2.0_pReal / pi ! see notes where it is used
if (iand(debug_what(debug_math),debug_levelBasic) /= 0_pInt) then
@@ -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
@@ -3355,7 +3350,7 @@ subroutine deformed_fft(res,geomdim,defgrad_av,scaling,defgrad,coords)
enddo; enddo; enddo
call fftw_execute_dft_r2c(fftw_forth, defgrad_real, defgrad_fourier)
-
+
!remove highest frequency in each direction
if(res(1)>1_pInt) &
defgrad_fourier( res(1)/2_pInt+1_pInt,1:res(2) ,1:res(3) ,&
@@ -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)
@@ -3384,7 +3378,7 @@ subroutine deformed_fft(res,geomdim,defgrad_av,scaling,defgrad,coords)
if(k/=1_pInt) coords_fourier(i,j,k,1:3) = coords_fourier(i,j,k,1:3)&
+ defgrad_fourier(i,j,k,1:3,3)*cmplx(0.0_pReal,integrator(3)/real(k_s(3),pReal),pReal)
enddo; enddo; enddo
-
+
call fftw_execute_dft_c2r(fftw_back,coords_fourier,coords_real)
coords_real = coords_real/real(res(1)*res(2)*res(3),pReal)
diff --git a/processing/pre/OIMang_hex2cub.py b/processing/pre/OIMang_hex2cub.py
new file mode 100755
index 000000000..30045f8ef
--- /dev/null
+++ b/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))
diff --git a/processing/pre/ang_hex2cub.py b/processing/pre/ang_hex2cub.py
deleted file mode 100644
index d700b8148..000000000
--- a/processing/pre/ang_hex2cub.py
+++ /dev/null
@@ -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
diff --git a/processing/pre/spectral_ang2geom.py b/processing/pre/spectral_ang2geom.py
new file mode 100755
index 000000000..e2b765b93
--- /dev/null
+++ b/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\n#---\n'+
+ '\n[SX]\ntype isostrain\nNgrains 1\n\n'+
+ '#---\n\n#---\n\n')
+ tempPart2= '#---\n\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])