From 385085de7333b9ab5e0fec47b620c5487600ba55 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 4 Jan 2020 01:05:40 +0100 Subject: [PATCH 01/37] correct names for numpy (differ from math) --- python/damask/colormaps.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/python/damask/colormaps.py b/python/damask/colormaps.py index ad924325f..e4183e830 100644 --- a/python/damask/colormaps.py +++ b/python/damask/colormaps.py @@ -328,9 +328,9 @@ class Color(): Msh = np.zeros(3,'d') Msh[0] = np.sqrt(np.dot(self.color,self.color)) if (Msh[0] > 0.001): - Msh[1] = np.acos(self.color[0]/Msh[0]) + Msh[1] = np.arccos(self.color[0]/Msh[0]) if (self.color[1] != 0.0): - Msh[2] = np.atan2(self.color[2],self.color[1]) + Msh[2] = np.arctan2(self.color[2],self.color[1]) converted = Color('MSH', Msh) self.model = converted.model From bd5f96326092c33402c3b07f1132c468108312c0 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 4 Jan 2020 14:37:09 +0100 Subject: [PATCH 02/37] polishing --- CMakeLists.txt | 11 +++-------- src/IO.f90 | 17 ++++++++--------- src/grid/grid_mech_FEM.f90 | 1 - 3 files changed, 11 insertions(+), 18 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 24079457a..e44f5eab2 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -117,9 +117,7 @@ elseif (DAMASK_SOLVER STREQUAL "fem" OR DAMASK_SOLVER STREQUAL "mesh") else () message (FATAL_ERROR "Build target (DAMASK_SOLVER) is not defined") endif () - -# set linker commands (needs to be done after defining the project) -set (CMAKE_LINKER "${PETSC_LINKER}") +list(APPEND CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/cmake) if (CMAKE_BUILD_TYPE STREQUAL "") set (CMAKE_BUILD_TYPE "RELEASE") @@ -168,9 +166,6 @@ add_definitions (-DDAMASKVERSION="${DAMASK_V}") # definition of other macros add_definitions (-DPETSc) -set (DAMASK_INCLUDE_FLAGS "${DAMASK_INCLUDE_FLAGS} ${PETSC_INCLUDES}") -list(APPEND CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/cmake) - if (CMAKE_Fortran_COMPILER_ID STREQUAL "Intel") include(Compiler-Intel) elseif(CMAKE_Fortran_COMPILER_ID STREQUAL "GNU") @@ -183,14 +178,14 @@ endif () set (CMAKE_Fortran_FLAGS_${CMAKE_BUILD_TYPE} "${BUILDCMD_PRE} ${OPENMP_FLAGS} ${STANDARD_CHECK} ${OPTIMIZATION_FLAGS} ${COMPILE_FLAGS} ${PRECISION_FLAGS}") -set (CMAKE_Fortran_LINK_EXECUTABLE "${BUILDCMD_PRE} ${CMAKE_LINKER} ${OPENMP_FLAGS} ${OPTIMIZATION_FLAGS} ${LINKER_FLAGS}") +set (CMAKE_Fortran_LINK_EXECUTABLE "${BUILDCMD_PRE} ${PETSC_LINKER} ${OPENMP_FLAGS} ${OPTIMIZATION_FLAGS} ${LINKER_FLAGS}") if (CMAKE_BUILD_TYPE STREQUAL "DEBUG") set (CMAKE_Fortran_FLAGS_${CMAKE_BUILD_TYPE} "${CMAKE_Fortran_FLAGS_${CMAKE_BUILD_TYPE}} ${DEBUG_FLAGS}") set (CMAKE_Fortran_LINK_EXECUTABLE "${CMAKE_Fortran_LINK_EXECUTABLE} ${DEBUG_FLAGS}") endif () -set (CMAKE_Fortran_FLAGS_${CMAKE_BUILD_TYPE} "${CMAKE_Fortran_FLAGS_${CMAKE_BUILD_TYPE}} ${DAMASK_INCLUDE_FLAGS} ${BUILDCMD_POST}") +set (CMAKE_Fortran_FLAGS_${CMAKE_BUILD_TYPE} "${CMAKE_Fortran_FLAGS_${CMAKE_BUILD_TYPE}} ${PETSC_INCLUDES} ${BUILDCMD_POST}") set (CMAKE_Fortran_LINK_EXECUTABLE "${CMAKE_Fortran_LINK_EXECUTABLE} -o ${PETSC_EXTERNAL_LIB} ${BUILDCMD_POST}") message ("Fortran Compiler Flags:\n${CMAKE_Fortran_FLAGS_${CMAKE_BUILD_TYPE}}\n") diff --git a/src/IO.f90 b/src/IO.f90 index 6d8506506..6c6c3b7c7 100644 --- a/src/IO.f90 +++ b/src/IO.f90 @@ -384,7 +384,7 @@ function IO_stringValue(string,chunkPos,myChunk,silent) logical :: warn if (present(silent)) then - warn = silent + warn = .not. silent else warn = .false. endif @@ -414,11 +414,10 @@ real(pReal) function IO_floatValue (string,chunkPos,myChunk) valuePresent: if (myChunk > chunkPos(1) .or. myChunk < 1) then call IO_warning(201,el=myChunk,ext_msg=MYNAME//trim(string)) - else valuePresent - IO_floatValue = & - verifyFloatValue(trim(adjustl(string(chunkPos(myChunk*2):chunkPos(myChunk*2+1)))),& - VALIDCHARACTERS,MYNAME) - endif valuePresent + else valuePresent + IO_floatValue = verifyFloatValue(trim(adjustl(string(chunkPos(myChunk*2):chunkPos(myChunk*2+1)))),& + VALIDCHARACTERS,MYNAME) + endif valuePresent end function IO_floatValue @@ -466,12 +465,12 @@ real(pReal) function IO_fixedNoEFloatValue (string,ends,myChunk) pos_exp = scan(string(ends(myChunk)+1:ends(myChunk+1)),'+-',back=.true.) hasExponent: if (pos_exp > 1) then base = verifyFloatValue(trim(adjustl(string(ends(myChunk)+1:ends(myChunk)+pos_exp-1))),& - VALIDBASE,MYNAME//'(base): ') + VALIDBASE,MYNAME//'(base): ') expon = verifyIntValue(trim(adjustl(string(ends(myChunk)+pos_exp:ends(myChunk+1)))),& - VALIDEXP,MYNAME//'(exp): ') + VALIDEXP,MYNAME//'(exp): ') else hasExponent base = verifyFloatValue(trim(adjustl(string(ends(myChunk)+1:ends(myChunk+1)))),& - VALIDBASE,MYNAME//'(base): ') + VALIDBASE,MYNAME//'(base): ') expon = 0 endif hasExponent IO_fixedNoEFloatValue = base*10.0_pReal**real(expon,pReal) diff --git a/src/grid/grid_mech_FEM.f90 b/src/grid/grid_mech_FEM.f90 index b2ea20f16..713e83029 100644 --- a/src/grid/grid_mech_FEM.f90 +++ b/src/grid/grid_mech_FEM.f90 @@ -94,7 +94,6 @@ subroutine grid_mech_FEM_init 1.0_pReal,-1.0_pReal,-1.0_pReal,-1.0_pReal, & 1.0_pReal, 1.0_pReal, 1.0_pReal, 1.0_pReal], [4,8]) PetscErrorCode :: ierr - integer :: rank integer(HID_T) :: fileHandle, groupHandle character(len=pStringLen) :: fileName real(pReal), dimension(3,3,3,3) :: devNull From 3999c0b63079bed881902f47a9f269269a612dda Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 4 Jan 2020 17:15:12 +0100 Subject: [PATCH 03/37] is not used anymore (and IO_fixedXXXvalue seem to be superfluous) --- src/IO.f90 | 43 ++----------------------------------------- src/mesh_marc.f90 | 7 +++---- 2 files changed, 5 insertions(+), 45 deletions(-) diff --git a/src/IO.f90 b/src/IO.f90 index 6c6c3b7c7..2d45a8add 100644 --- a/src/IO.f90 +++ b/src/IO.f90 @@ -23,7 +23,7 @@ module IO public :: & IO_init, & IO_read_ASCII, & - IO_open_file, & + IO_open_file, & ! deprecated, use IO_read_ASCII IO_open_jobFile_binary, & IO_isBlank, & IO_getTag, & @@ -44,8 +44,7 @@ module IO IO_countDataLines #elif defined(Marc4DAMASK) IO_fixedNoEFloatValue, & - IO_fixedIntValue, & - IO_countNumericalDataLines + IO_fixedIntValue #endif #endif @@ -549,14 +548,8 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg) msg = 'could not read file:' case (103) msg = 'could not assemble input files' - case (104) - msg = '{input} recursion limit reached' - case (105) - msg = 'unknown output:' case (106) msg = 'working directory does not exist:' - case (107) - msg = 'line length exceeds limit of 256' !-------------------------------------------------------------------------------------------------- ! lattice error messages @@ -923,38 +916,6 @@ end function IO_countDataLines #endif -#ifdef Marc4DAMASK -!-------------------------------------------------------------------------------------------------- -!> @brief count lines containig data up to next *keyword -!-------------------------------------------------------------------------------------------------- -integer function IO_countNumericalDataLines(fileUnit) - - integer, intent(in) :: fileUnit !< file handle - - - integer, allocatable, dimension(:) :: chunkPos - character(len=pStringLen) :: line, & - tmp - - IO_countNumericalDataLines = 0 - line = '' - - do while (trim(line) /= IO_EOF) - line = IO_read(fileUnit) - chunkPos = IO_stringPos(line) - tmp = IO_lc(IO_stringValue(line,chunkPos,1)) - if (verify(trim(tmp),'0123456789') == 0) then ! numerical values - IO_countNumericalDataLines = IO_countNumericalDataLines + 1 - else - exit - endif - enddo - backspace(fileUnit) - -end function IO_countNumericalDataLines -#endif - - !-------------------------------------------------------------------------------------------------- !> @brief count items in consecutive lines depending on lines !> @details Marc: ints concatenated by "c" as last char or range of values a "to" b diff --git a/src/mesh_marc.f90 b/src/mesh_marc.f90 index 1d7fd09c3..00fd76326 100644 --- a/src/mesh_marc.f90 +++ b/src/mesh_marc.f90 @@ -488,8 +488,7 @@ subroutine inputRead_mapNodes(fileContent) chunkPos = IO_stringPos(fileContent(l)) if( IO_lc(IO_stringValue(fileContent(l),chunkPos,1)) == 'coordinates' ) then do i = 1,size(mesh_mapFEtoCPnode,2) - mesh_mapFEtoCPnode(1,i) = IO_fixedIntValue (fileContent(l+1+i),[0,10],1) - mesh_mapFEtoCPnode(2,i) = i + mesh_mapFEtoCPnode(1:2,i) = [IO_fixedIntValue (fileContent(l+1+i),[0,10],1),i] ! ToDo: use IO_intValue enddo exit endif @@ -520,9 +519,9 @@ subroutine inputRead_elemNodes(nodes, & chunkPos = IO_stringPos(fileContent(l)) if( IO_lc(IO_stringValue(fileContent(l),chunkPos,1)) == 'coordinates' ) then do i=1,nNode - m = mesh_FEasCP('node',IO_fixedIntValue(fileContent(l+1+i),node_ends,1)) + m = mesh_FEasCP('node',IO_fixedIntValue(fileContent(l+1+i),node_ends,1)) !ToDo: use IO_intValue do j = 1,3 - nodes(j,m) = mesh_unitlength * IO_fixedNoEFloatValue(fileContent(l+1+i),node_ends,j+1) + nodes(j,m) = mesh_unitlength * IO_fixedNoEFloatValue(fileContent(l+1+i),node_ends,j+1) !ToDo: use IO_floatValue enddo enddo exit From 6b6ad5235535d8dbf9dcb4154f854269d49c51d8 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 4 Jan 2020 18:25:59 +0100 Subject: [PATCH 04/37] use variable string as return (no need for trim) --- src/DAMASK_interface.f90 | 37 ++++++++++++++++++++----------------- src/system_routines.f90 | 14 ++++++++++---- 2 files changed, 30 insertions(+), 21 deletions(-) diff --git a/src/DAMASK_interface.f90 b/src/DAMASK_interface.f90 index 4cf155ac1..27a0084f5 100644 --- a/src/DAMASK_interface.f90 +++ b/src/DAMASK_interface.f90 @@ -269,10 +269,10 @@ subroutine DAMASK_interface_init write(6,'(a,a)') ' Working dir argument: ', trim(workingDirArg) write(6,'(a,a)') ' Geometry argument: ', trim(geometryArg) write(6,'(a,a)') ' Load case argument: ', trim(loadcaseArg) - write(6,'(a,a)') ' Working directory: ', trim(getCWD()) + write(6,'(a,a)') ' Working directory: ', getCWD() write(6,'(a,a)') ' Geometry file: ', trim(geometryFile) write(6,'(a,a)') ' Loadcase file: ', trim(loadCaseFile) - write(6,'(a,a)') ' Solver job name: ', trim(getSolverJobName()) + write(6,'(a,a)') ' Solver job name: ', getSolverJobName() if (interface_restartInc > 0) & write(6,'(a,i6.6)') ' Restart from increment: ', interface_restartInc @@ -308,7 +308,7 @@ subroutine setWorkingDirectory(workingDirectoryArg) workingDirectory = trim(rectifyPath(workingDirectory)) error = setCWD(trim(workingDirectory)) if(error) then - write(6,'(/,a)') ' ERROR: Working directory "'//trim(workingDirectory)//'" does not exist' + write(6,'(/,a)') ' ERROR: Invalid Working directory: '//trim(workingDirectory) call quit(1) endif @@ -318,8 +318,9 @@ end subroutine setWorkingDirectory !-------------------------------------------------------------------------------------------------- !> @brief solver job name (no extension) as combination of geometry and load case name !-------------------------------------------------------------------------------------------------- -character(len=1024) function getSolverJobName() +function getSolverJobName() + character(len=:), allocatable :: getSolverJobName integer :: posExt,posSep posExt = scan(geometryFile,'.',back=.true.) @@ -330,7 +331,7 @@ character(len=1024) function getSolverJobName() posExt = scan(loadCaseFile,'.',back=.true.) posSep = scan(loadCaseFile,'/',back=.true.) - getSolverJobName = trim(getSolverJobName)//'_'//loadCaseFile(posSep+1:posExt-1) + getSolverJobName = getSolverJobName//'_'//loadCaseFile(posSep+1:posExt-1) end function getSolverJobName @@ -338,15 +339,16 @@ end function getSolverJobName !-------------------------------------------------------------------------------------------------- !> @brief basename of geometry file with extension from command line arguments !-------------------------------------------------------------------------------------------------- -character(len=1024) function getGeometryFile(geometryParameter) +function getGeometryFile(geometryParameter) - character(len=1024), intent(in) :: geometryParameter - logical :: file_exists - external :: quit + character(len=:), allocatable :: getGeometryFile + character(len=*), intent(in) :: geometryParameter + logical :: file_exists + external :: quit getGeometryFile = trim(geometryParameter) - if (scan(getGeometryFile,'/') /= 1) getGeometryFile = trim(getCWD())//'/'//trim(getGeometryFile) - getGeometryFile = makeRelativePath(trim(getCWD()), getGeometryFile) + if (scan(getGeometryFile,'/') /= 1) getGeometryFile = getCWD()//'/'//trim(getGeometryFile) + getGeometryFile = makeRelativePath(getCWD(), getGeometryFile) inquire(file=trim(getGeometryFile), exist=file_exists) if (.not. file_exists) then @@ -360,15 +362,16 @@ end function getGeometryFile !-------------------------------------------------------------------------------------------------- !> @brief relative path of loadcase from command line arguments !-------------------------------------------------------------------------------------------------- -character(len=1024) function getLoadCaseFile(loadCaseParameter) +function getLoadCaseFile(loadCaseParameter) - character(len=1024), intent(in) :: loadCaseParameter - logical :: file_exists - external :: quit + character(len=:), allocatable :: getLoadCaseFile + character(len=*), intent(in) :: loadCaseParameter + logical :: file_exists + external :: quit getLoadCaseFile = trim(loadCaseParameter) - if (scan(getLoadCaseFile,'/') /= 1) getLoadCaseFile = trim(getCWD())//'/'//trim(getLoadCaseFile) - getLoadCaseFile = makeRelativePath(trim(getCWD()), getLoadCaseFile) + if (scan(getLoadCaseFile,'/') /= 1) getLoadCaseFile = getCWD()//'/'//trim(getLoadCaseFile) + getLoadCaseFile = makeRelativePath(getCWD(), getLoadCaseFile) inquire(file=trim(getLoadCaseFile), exist=file_exists) if (.not. file_exists) then diff --git a/src/system_routines.f90 b/src/system_routines.f90 index 0611c96db..932eefeb6 100644 --- a/src/system_routines.f90 +++ b/src/system_routines.f90 @@ -93,21 +93,24 @@ end function isDirectory !-------------------------------------------------------------------------------------------------- !> @brief gets the current working directory !-------------------------------------------------------------------------------------------------- -character(len=1024) function getCWD() +function getCWD() character(kind=C_CHAR), dimension(1024) :: charArray ! C string is an array + character(len=:), allocatable :: getCWD integer(C_INT) :: stat integer :: i call getCurrentWorkDir_C(charArray,stat) + if (stat /= 0_C_INT) then getCWD = 'Error occured when getting currend working directory' else - getCWD = repeat('',len(getCWD)) + allocate(character(len=1024)::getCWD) arrayToString: do i=1,len(getCWD) if (charArray(i) /= C_NULL_CHAR) then getCWD(i:i)=charArray(i) else + getCWD = getCWD(:i-1) exit endif enddo arrayToString @@ -119,21 +122,24 @@ end function getCWD !-------------------------------------------------------------------------------------------------- !> @brief gets the current host name !-------------------------------------------------------------------------------------------------- -character(len=1024) function getHostName() +function getHostName() character(kind=C_CHAR), dimension(1024) :: charArray ! C string is an array + character(len=:), allocatable :: getHostName integer(C_INT) :: stat integer :: i call getHostName_C(charArray,stat) + if (stat /= 0_C_INT) then getHostName = 'Error occured when getting host name' else - getHostName = repeat('',len(getHostName)) + allocate(character(len=1024)::getHostName) arrayToString: do i=1,len(getHostName) if (charArray(i) /= C_NULL_CHAR) then getHostName(i:i)=charArray(i) else + getHostName = getHostName(:i-1) exit endif enddo arrayToString From bd6f2a6b5c0b26adbce7129df25d854160f46834 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 4 Jan 2020 19:01:36 +0100 Subject: [PATCH 05/37] consistent string length --- src/IO.f90 | 8 ++++---- src/grid/grid_mech_FEM.f90 | 2 +- src/grid/grid_mech_spectral_basic.f90 | 2 +- src/grid/grid_mech_spectral_polarisation.f90 | 2 +- src/grid/spectral_utilities.f90 | 2 +- 5 files changed, 8 insertions(+), 8 deletions(-) diff --git a/src/IO.f90 b/src/IO.f90 index 2d45a8add..1e39daa1e 100644 --- a/src/IO.f90 +++ b/src/IO.f90 @@ -528,8 +528,8 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg) character(len=*), optional, intent(in) :: ext_msg external :: quit - character(len=1024) :: msg - character(len=1024) :: formatString + character(len=pStringLen) :: msg + character(len=pStringLen) :: formatString select case (error_ID) @@ -769,8 +769,8 @@ subroutine IO_warning(warning_ID,el,ip,g,ext_msg) integer, optional, intent(in) :: el,ip,g character(len=*), optional, intent(in) :: ext_msg - character(len=1024) :: msg - character(len=1024) :: formatString + character(len=pStringLen) :: msg + character(len=pStringLen) :: formatString select case (warning_ID) case (1) diff --git a/src/grid/grid_mech_FEM.f90 b/src/grid/grid_mech_FEM.f90 index 713e83029..05fc3520a 100644 --- a/src/grid/grid_mech_FEM.f90 +++ b/src/grid/grid_mech_FEM.f90 @@ -53,7 +53,7 @@ module grid_mech_FEM F_aim_lastInc = math_I3, & !< previous average deformation gradient P_av = 0.0_pReal !< average 1st Piola--Kirchhoff stress - character(len=1024), private :: incInfo !< time and increment information + character(len=pStringLen), private :: incInfo !< time and increment information real(pReal), private, dimension(3,3,3,3) :: & C_volAvg = 0.0_pReal, & !< current volume average stiffness diff --git a/src/grid/grid_mech_spectral_basic.f90 b/src/grid/grid_mech_spectral_basic.f90 index 9fcd3d563..af8cbc377 100644 --- a/src/grid/grid_mech_spectral_basic.f90 +++ b/src/grid/grid_mech_spectral_basic.f90 @@ -55,7 +55,7 @@ module grid_mech_spectral_basic F_aim_lastInc = math_I3, & !< previous average deformation gradient P_av = 0.0_pReal !< average 1st Piola--Kirchhoff stress - character(len=1024), private :: incInfo !< time and increment information + character(len=pStringLen), private :: incInfo !< time and increment information real(pReal), private, dimension(3,3,3,3) :: & C_volAvg = 0.0_pReal, & !< current volume average stiffness C_volAvgLastInc = 0.0_pReal, & !< previous volume average stiffness diff --git a/src/grid/grid_mech_spectral_polarisation.f90 b/src/grid/grid_mech_spectral_polarisation.f90 index 4986a3457..59ab84869 100644 --- a/src/grid/grid_mech_spectral_polarisation.f90 +++ b/src/grid/grid_mech_spectral_polarisation.f90 @@ -59,7 +59,7 @@ module grid_mech_spectral_polarisation F_av = 0.0_pReal, & !< average incompatible def grad field P_av = 0.0_pReal !< average 1st Piola--Kirchhoff stress - character(len=1024), private :: incInfo !< time and increment information + character(len=pStringLen), private :: incInfo !< time and increment information real(pReal), private, dimension(3,3,3,3) :: & C_volAvg = 0.0_pReal, & !< current volume average stiffness C_volAvgLastInc = 0.0_pReal, & !< previous volume average stiffness diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index 04a24e877..68c372b26 100644 --- a/src/grid/spectral_utilities.f90 +++ b/src/grid/spectral_utilities.f90 @@ -700,7 +700,7 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C) c_reduced, & !< reduced stiffness (depending on number of stress BC) sTimesC !< temp variable to check inversion logical :: errmatinv - character(len=1024):: formatString + character(len=pStringLen):: formatString mask_stressVector = reshape(transpose(mask_stress), [9]) size_reduced = count(mask_stressVector) From 09f42a399136df4747c807544fc624bd5dfcd539 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 8 Jan 2020 15:34:21 +0100 Subject: [PATCH 06/37] ang files might have more columns --- python/damask/table.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/python/damask/table.py b/python/damask/table.py index 28ce3efe0..ef8a84276 100644 --- a/python/damask/table.py +++ b/python/damask/table.py @@ -138,9 +138,12 @@ class Table(): break data = np.loadtxt(content) + for c in range(data.shape[1]-10): + shapes['user_defined{}'.format(c+1)] = (1,) return Table(data,shapes,comments) + @property def labels(self): return list(self.shapes.keys()) From 70e23fea9365c74efefcc58c0199f984286595c6 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 10 Jan 2020 01:28:32 +0100 Subject: [PATCH 07/37] polishing --- src/kinematics_cleavage_opening.f90 | 2 +- src/kinematics_slipplane_opening.f90 | 29 ++++++++++------------------ src/kinematics_thermal_expansion.f90 | 2 +- 3 files changed, 12 insertions(+), 21 deletions(-) diff --git a/src/kinematics_cleavage_opening.f90 b/src/kinematics_cleavage_opening.f90 index eec8a1986..6b060a8d9 100644 --- a/src/kinematics_cleavage_opening.f90 +++ b/src/kinematics_cleavage_opening.f90 @@ -65,7 +65,7 @@ subroutine kinematics_cleavage_opening_init integer :: maxNinstance,p,instance - write(6,'(/,a)') ' <<<+- kinematics_'//KINEMATICS_cleavage_opening_LABEL//' init -+>>>' + write(6,'(/,a)') ' <<<+- kinematics_'//KINEMATICS_cleavage_opening_LABEL//' init -+>>>'; flush(6) maxNinstance = count(phase_kinematics == KINEMATICS_cleavage_opening_ID) if (maxNinstance == 0) return diff --git a/src/kinematics_slipplane_opening.f90 b/src/kinematics_slipplane_opening.f90 index c0f198985..018a5b4b5 100644 --- a/src/kinematics_slipplane_opening.f90 +++ b/src/kinematics_slipplane_opening.f90 @@ -51,7 +51,7 @@ subroutine kinematics_slipplane_opening_init integer :: maxNinstance,p,instance - write(6,'(/,a)') ' <<<+- kinematics_'//KINEMATICS_slipplane_opening_LABEL//' init -+>>>' + write(6,'(/,a)') ' <<<+- kinematics_'//KINEMATICS_slipplane_opening_LABEL//' init -+>>>'; flush(6) maxNinstance = count(phase_kinematics == KINEMATICS_slipplane_opening_ID) if (maxNinstance == 0) return @@ -144,40 +144,31 @@ subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, traction_crit = prm%critLoad(i)* damage(homog)%p(damageOffset) ! degrading critical load carrying capacity by damage - udotd = sign(1.0_pReal,traction_d)* & - prm%sdot0* & - (abs(traction_d)/traction_crit - & - abs(traction_d)/prm%critLoad(i))**prm%n + udotd = sign(1.0_pReal,traction_d)* prm%sdot0* ( abs(traction_d)/traction_crit & + - abs(traction_d)/prm%critLoad(i))**prm%n if (abs(udotd) > tol_math_check) then Ld = Ld + udotd*projection_d dudotd_dt = udotd*prm%n/traction_d forall (k=1:3,l=1:3,m=1:3,n=1:3) & - dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + & - dudotd_dt*projection_d(k,l)*projection_d(m,n) + dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + dudotd_dt*projection_d(k,l)*projection_d(m,n) endif - udott = sign(1.0_pReal,traction_t)* & - prm%sdot0* & - (abs(traction_t)/traction_crit - & - abs(traction_t)/prm%critLoad(i))**prm%n + udott = sign(1.0_pReal,traction_t)* prm%sdot0* ( abs(traction_t)/traction_crit & + - abs(traction_t)/prm%critLoad(i))**prm%n if (abs(udott) > tol_math_check) then Ld = Ld + udott*projection_t dudott_dt = udott*prm%n/traction_t forall (k=1:3,l=1:3,m=1:3,n=1:3) & - dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + & - dudott_dt*projection_t(k,l)*projection_t(m,n) + dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + dudott_dt*projection_t(k,l)*projection_t(m,n) endif - udotn = & - prm%sdot0* & - (max(0.0_pReal,traction_n)/traction_crit - & - max(0.0_pReal,traction_n)/prm%critLoad(i))**prm%n + udotn = prm%sdot0* ( max(0.0_pReal,traction_n)/traction_crit & + - max(0.0_pReal,traction_n)/prm%critLoad(i))**prm%n if (abs(udotn) > tol_math_check) then Ld = Ld + udotn*projection_n dudotn_dt = udotn*prm%n/traction_n forall (k=1:3,l=1:3,m=1:3,n=1:3) & - dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + & - dudotn_dt*projection_n(k,l)*projection_n(m,n) + dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + dudotn_dt*projection_n(k,l)*projection_n(m,n) endif enddo diff --git a/src/kinematics_thermal_expansion.f90 b/src/kinematics_thermal_expansion.f90 index 814d604ed..7f7994959 100644 --- a/src/kinematics_thermal_expansion.f90 +++ b/src/kinematics_thermal_expansion.f90 @@ -42,7 +42,7 @@ subroutine kinematics_thermal_expansion_init real(pReal), dimension(:), allocatable :: & temp - write(6,'(/,a)') ' <<<+- kinematics_'//KINEMATICS_thermal_expansion_LABEL//' init -+>>>' + write(6,'(/,a)') ' <<<+- kinematics_'//KINEMATICS_thermal_expansion_LABEL//' init -+>>>'; flush(6) Ninstance = count(phase_kinematics == KINEMATICS_thermal_expansion_ID) From 7d2012f492021243e3d8c2bd8f1f51369150c075 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 10 Jan 2020 01:29:35 +0100 Subject: [PATCH 08/37] no need to exclude small values no danger of division by zero --- src/kinematics_slipplane_opening.f90 | 30 +++++++++++----------------- 1 file changed, 12 insertions(+), 18 deletions(-) diff --git a/src/kinematics_slipplane_opening.f90 b/src/kinematics_slipplane_opening.f90 index 018a5b4b5..a736338f9 100644 --- a/src/kinematics_slipplane_opening.f90 +++ b/src/kinematics_slipplane_opening.f90 @@ -146,30 +146,24 @@ subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, udotd = sign(1.0_pReal,traction_d)* prm%sdot0* ( abs(traction_d)/traction_crit & - abs(traction_d)/prm%critLoad(i))**prm%n - if (abs(udotd) > tol_math_check) then - Ld = Ld + udotd*projection_d - dudotd_dt = udotd*prm%n/traction_d - forall (k=1:3,l=1:3,m=1:3,n=1:3) & - dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + dudotd_dt*projection_d(k,l)*projection_d(m,n) - endif + Ld = Ld + udotd*projection_d + dudotd_dt = udotd*prm%n/traction_d + forall (k=1:3,l=1:3,m=1:3,n=1:3) & + dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + dudotd_dt*projection_d(k,l)*projection_d(m,n) udott = sign(1.0_pReal,traction_t)* prm%sdot0* ( abs(traction_t)/traction_crit & - abs(traction_t)/prm%critLoad(i))**prm%n - if (abs(udott) > tol_math_check) then - Ld = Ld + udott*projection_t - dudott_dt = udott*prm%n/traction_t - forall (k=1:3,l=1:3,m=1:3,n=1:3) & - dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + dudott_dt*projection_t(k,l)*projection_t(m,n) - endif + Ld = Ld + udott*projection_t + dudott_dt = udott*prm%n/traction_t + forall (k=1:3,l=1:3,m=1:3,n=1:3) & + dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + dudott_dt*projection_t(k,l)*projection_t(m,n) udotn = prm%sdot0* ( max(0.0_pReal,traction_n)/traction_crit & - max(0.0_pReal,traction_n)/prm%critLoad(i))**prm%n - if (abs(udotn) > tol_math_check) then - Ld = Ld + udotn*projection_n - dudotn_dt = udotn*prm%n/traction_n - forall (k=1:3,l=1:3,m=1:3,n=1:3) & - dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + dudotn_dt*projection_n(k,l)*projection_n(m,n) - endif + Ld = Ld + udotn*projection_n + dudotn_dt = udotn*prm%n/traction_n + forall (k=1:3,l=1:3,m=1:3,n=1:3) & + dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + dudotn_dt*projection_n(k,l)*projection_n(m,n) enddo end associate From 0771411cd82a66ff7504748a731f1035bf4bf8f3 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 10 Jan 2020 01:33:03 +0100 Subject: [PATCH 09/37] group similar operations --- src/kinematics_slipplane_opening.f90 | 29 ++++++++++++++-------------- 1 file changed, 14 insertions(+), 15 deletions(-) diff --git a/src/kinematics_slipplane_opening.f90 b/src/kinematics_slipplane_opening.f90 index a736338f9..2c94448bd 100644 --- a/src/kinematics_slipplane_opening.f90 +++ b/src/kinematics_slipplane_opening.f90 @@ -134,36 +134,35 @@ subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, dLd_dTstar = 0.0_pReal do i = 1, prm%totalNslip - projection_d = math_outer(prm%slip_direction(1:3,i),prm%slip_normal(1:3,i)) + projection_d = math_outer(prm%slip_direction(1:3,i), prm%slip_normal(1:3,i)) projection_t = math_outer(prm%slip_transverse(1:3,i),prm%slip_normal(1:3,i)) - projection_n = math_outer(prm%slip_normal(1:3,i),prm%slip_normal(1:3,i)) + projection_n = math_outer(prm%slip_normal(1:3,i), prm%slip_normal(1:3,i)) traction_d = math_mul33xx33(S,projection_d) traction_t = math_mul33xx33(S,projection_t) traction_n = math_mul33xx33(S,projection_n) - traction_crit = prm%critLoad(i)* damage(homog)%p(damageOffset) ! degrading critical load carrying capacity by damage + traction_crit = prm%critLoad(i)* damage(homog)%p(damageOffset) ! degrading critical load carrying capacity by damage udotd = sign(1.0_pReal,traction_d)* prm%sdot0* ( abs(traction_d)/traction_crit & - abs(traction_d)/prm%critLoad(i))**prm%n - Ld = Ld + udotd*projection_d - dudotd_dt = udotd*prm%n/traction_d - forall (k=1:3,l=1:3,m=1:3,n=1:3) & - dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + dudotd_dt*projection_d(k,l)*projection_d(m,n) - udott = sign(1.0_pReal,traction_t)* prm%sdot0* ( abs(traction_t)/traction_crit & - abs(traction_t)/prm%critLoad(i))**prm%n - Ld = Ld + udott*projection_t - dudott_dt = udott*prm%n/traction_t - forall (k=1:3,l=1:3,m=1:3,n=1:3) & - dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + dudott_dt*projection_t(k,l)*projection_t(m,n) - udotn = prm%sdot0* ( max(0.0_pReal,traction_n)/traction_crit & - max(0.0_pReal,traction_n)/prm%critLoad(i))**prm%n - Ld = Ld + udotn*projection_n + + dudotd_dt = udotd*prm%n/traction_d + dudott_dt = udott*prm%n/traction_t dudotn_dt = udotn*prm%n/traction_n + forall (k=1:3,l=1:3,m=1:3,n=1:3) & - dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + dudotn_dt*projection_n(k,l)*projection_n(m,n) + dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + dudotd_dt*projection_d(k,l)*projection_d(m,n) & + + dudott_dt*projection_t(k,l)*projection_t(m,n) & + + dudotn_dt*projection_n(k,l)*projection_n(m,n) + + Ld = Ld + udotd*projection_d & + + udott*projection_t & + + udotn*projection_n enddo end associate From 5d2802f6c98eb952b6d5af69d610fc11234d7acd Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 10 Jan 2020 01:36:48 +0100 Subject: [PATCH 10/37] no patches available at the moment --- installation/patch/README.md | 8 -------- 1 file changed, 8 deletions(-) diff --git a/installation/patch/README.md b/installation/patch/README.md index 0b8251510..95f377691 100644 --- a/installation/patch/README.md +++ b/installation/patch/README.md @@ -9,14 +9,6 @@ cd DAMASK_ROOT patch -p1 < installation/patch/nameOfPatch ``` -## Available patches - - * **disable_HDF5** disables all HDF5 output. - HDF5 output is an experimental feature. Also, some routines not present in HDF5 1.8.x are removed to allow compilation of DAMASK with HDF5 < 1.10.x - - * **disable_old_output** disables all non-HDF5 output. - Saves some memory when using only HDF5 output - ## Create patch commit your changes From 8f43f054375b4689fbf19389baf343c8efc6a262 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 10 Jan 2020 01:45:00 +0100 Subject: [PATCH 11/37] stronger encapsulation --- src/CPFEM.f90 | 2 +- src/CPFEM2.f90 | 2 +- src/results.f90 | 12 ++++++++++++ 3 files changed, 14 insertions(+), 2 deletions(-) diff --git a/src/CPFEM.f90 b/src/CPFEM.f90 index 3a7f35633..d74155703 100644 --- a/src/CPFEM.f90 +++ b/src/CPFEM.f90 @@ -377,7 +377,7 @@ subroutine CPFEM_results(inc,time) call constitutive_results call crystallite_results call homogenization_results - call results_removeLink('current') ! ToDo: put this into closeJobFile + call results_finalizeIncrement call results_closeJobFile end subroutine CPFEM_results diff --git a/src/CPFEM2.f90 b/src/CPFEM2.f90 index 6b4479b89..bc3a424e3 100644 --- a/src/CPFEM2.f90 +++ b/src/CPFEM2.f90 @@ -201,7 +201,7 @@ subroutine CPFEM_results(inc,time) call crystallite_results call homogenization_results call discretization_results - call results_removeLink('current') ! ToDo: put this into closeJobFile? + call results_finalizeIncrement call results_closeJobFile end subroutine CPFEM_results diff --git a/src/results.f90 b/src/results.f90 index 73c5d7dbb..ca2a868fb 100644 --- a/src/results.f90 +++ b/src/results.f90 @@ -47,6 +47,7 @@ module results results_openJobFile, & results_closeJobFile, & results_addIncrement, & + results_finalizeIncrement, & results_addGroup, & results_openGroup, & results_closeGroup, & @@ -119,6 +120,17 @@ subroutine results_addIncrement(inc,time) end subroutine results_addIncrement +!-------------------------------------------------------------------------------------------------- +!> @brief finalize increment +!> @details remove soft link +!-------------------------------------------------------------------------------------------------- +subroutine results_finalizeIncrement + + call results_removeLink('current') + +end subroutine results_finalizeIncrement + + !-------------------------------------------------------------------------------------------------- !> @brief open a group from the results file !-------------------------------------------------------------------------------------------------- From 87c7a5d5a32fbd810f8d3800f0ed7f43cb12ba80 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 10 Jan 2020 02:11:19 +0100 Subject: [PATCH 12/37] polishing --- src/math.f90 | 30 +++++++++++++++++------------- 1 file changed, 17 insertions(+), 13 deletions(-) diff --git a/src/math.f90 b/src/math.f90 index d6d87d30a..62c22b6d8 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -398,19 +398,19 @@ pure function math_exp33(A,n) real(pReal), dimension(3,3), intent(in) :: A real(pReal), dimension(3,3) :: B, math_exp33 real(pReal) :: invFac - integer :: order - - B = math_I3 ! init - invFac = 1.0_pReal ! 0! - math_exp33 = B ! A^0 = eye2 + integer :: n_ if (present(n)) then - order = n + n_ = n else - order = 5 + n_ = 5 endif - - do i = 1, order + + invFac = 1.0_pReal ! 0! + B = math_I3 + math_exp33 = math_I3 ! A^0 = I + + do i = 1, n_ invFac = invFac/real(i,pReal) ! invfac = 1/(i!) B = matmul(B,A) math_exp33 = math_exp33 + invFac*B ! exp = SUM (A^i)/(i!) @@ -882,16 +882,20 @@ real(pReal) function math_sampleGaussVar(meanvalue, stddev, width) real(pReal), intent(in), optional :: width ! width of considered values as multiples of standard deviation real(pReal), dimension(2) :: rnd ! random numbers real(pReal) :: scatter, & ! normalized scatter around meanvalue - myWidth + width_ if (abs(stddev) < tol_math_check) then math_sampleGaussVar = meanvalue else - myWidth = merge(width,3.0_pReal,present(width)) ! use +-3*sigma as default value for scatter if not given - + if (present(width)) then + width_ = width + else + width_ = 3.0_pReal ! use +-3*sigma as default scatter + endif + do call random_number(rnd) - scatter = myWidth * (2.0_pReal * rnd(1) - 1.0_pReal) + scatter = width_ * (2.0_pReal * rnd(1) - 1.0_pReal) if (rnd(2) <= exp(-0.5_pReal * scatter ** 2.0_pReal)) exit ! test if scattered value is drawn enddo From 115a2552f8f0b2fd3b5eda8b983ae880bacdfc44 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 10 Jan 2020 03:19:39 +0100 Subject: [PATCH 13/37] 4 newer versions are out --- src/element.f90 | 22 ++++++++-------------- 1 file changed, 8 insertions(+), 14 deletions(-) diff --git a/src/element.f90 b/src/element.f90 index 02f5fb762..3a1e3f5a3 100644 --- a/src/element.f90 +++ b/src/element.f90 @@ -39,7 +39,7 @@ module element integer, parameter, private :: & NELEMTYPE = 13 - integer, dimension(NelemType), parameter, private :: NNODE = & + integer, dimension(NELEMTYPE), parameter, private :: NNODE = & [ & 3, & ! 2D 3node 1ip 6, & ! 2D 6node 3ip @@ -57,7 +57,7 @@ module element 20 & ! 3D 20node 27ip ] !< number of nodes that constitute a specific type of element - integer, dimension(NelemType), parameter, public :: GEOMTYPE = & + integer, dimension(NELEMTYPE), parameter, public :: GEOMTYPE = & [ & 1, & 2, & @@ -74,8 +74,7 @@ module element 10 & ] !< geometry type of particular element type - !integer, dimension(maxval(geomType)), parameter, private :: NCELLNODE = & ! Intel 16.0 complains - integer, dimension(10), parameter, private :: NCELLNODE = & + integer, dimension(maxval(GEOMTYPE)), parameter, private :: NCELLNODE = & [ & 3, & 7, & @@ -89,8 +88,7 @@ module element 64 & ] !< number of cell nodes in a specific geometry type - !integer, dimension(maxval(geomType)), parameter, private :: NIP = & ! Intel 16.0 complains - integer, dimension(10), parameter, private :: NIP = & + integer, dimension(maxval(GEOMTYPE)), parameter, private :: NIP = & [ & 1, & 3, & @@ -104,8 +102,7 @@ module element 27 & ] !< number of IPs in a specific geometry type - !integer, dimension(maxval(geomType)), parameter, private :: CELLTYPE = & ! Intel 16.0 complains - integer, dimension(10), parameter, private :: CELLTYPE = & + integer, dimension(maxval(GEOMTYPE)), parameter, private :: CELLTYPE = & [ & 1, & ! 2D 3node 2, & ! 2D 4node @@ -119,8 +116,7 @@ module element 4 & ! 3D 8node ] !< cell type that is used by each geometry type - !integer, dimension(maxval(cellType)), parameter, private :: nIPNeighbor = & ! Intel 16.0 complains - integer, dimension(4), parameter, private :: NIPNEIGHBOR = & + integer, dimension(maxval(CELLTYPE)), parameter, private :: NIPNEIGHBOR = & [ & 3, & ! 2D 3node 4, & ! 2D 4node @@ -128,8 +124,7 @@ module element 6 & ! 3D 8node ] !< number of ip neighbors / cell faces in a specific cell type - !integer, dimension(maxval(cellType)), parameter, private :: NCELLNODESPERCELLFACE = & ! Intel 16.0 complains - integer, dimension(4), parameter, private :: NCELLNODEPERCELLFACE = & + integer, dimension(maxval(CELLTYPE)), parameter, private :: NCELLNODEPERCELLFACE = & [ & 2, & ! 2D 3node 2, & ! 2D 4node @@ -137,8 +132,7 @@ module element 4 & ! 3D 8node ] !< number of cell nodes in a specific cell type - !integer, dimension(maxval(CELLTYPE)), parameter, private :: NCELLNODEPERCELL = & ! Intel 16.0 complains - integer, dimension(4), parameter, private :: NCELLNODEPERCELL = & + integer, dimension(maxval(CELLTYPE)), parameter, private :: NCELLNODEPERCELL = & [ & 3, & ! 2D 3node 4, & ! 2D 4node From 3a08a8bbe27ec3b6b74bbda5516ec97747f5b6ed Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Fri, 10 Jan 2020 12:07:30 -0500 Subject: [PATCH 14/37] always using intrinsic init when assigning quaternions as output variables --- src/quaternions.f90 | 68 +++++++++++++++++---------------------------- 1 file changed, 26 insertions(+), 42 deletions(-) diff --git a/src/quaternions.f90 b/src/quaternions.f90 index 8efb985ed..e55d3804e 100644 --- a/src/quaternions.f90 +++ b/src/quaternions.f90 @@ -142,32 +142,29 @@ type(quaternion) pure function init__(array) real(pReal), intent(in), dimension(4) :: array - init__%w=array(1) - init__%x=array(2) - init__%y=array(3) - init__%z=array(4) + init__%w = array(1) + init__%x = array(2) + init__%y = array(3) + init__%z = array(4) end function init__ !--------------------------------------------------------------------------------------------------- -!> assing a quaternion +!> assigning a quaternion !--------------------------------------------------------------------------------------------------- elemental pure subroutine assign_quat__(self,other) type(quaternion), intent(out) :: self type(quaternion), intent(in) :: other - self%w = other%w - self%x = other%x - self%y = other%y - self%z = other%z - + self = [other%w,other%x,other%y,other%z] + end subroutine assign_quat__ !--------------------------------------------------------------------------------------------------- -!> assing a 4-vector +!> assigning a 4-vector !--------------------------------------------------------------------------------------------------- pure subroutine assign_vec__(self,other) @@ -189,11 +186,8 @@ type(quaternion) elemental pure function add__(self,other) class(quaternion), intent(in) :: self,other - add__%w = self%w + other%w - add__%x = self%x + other%x - add__%y = self%y + other%y - add__%z = self%z + other%z - + add__ = [self%w + other%w,self%x + other%x,self%y + other%y,self%z + other%z] + end function add__ @@ -204,11 +198,8 @@ type(quaternion) elemental pure function pos__(self) class(quaternion), intent(in) :: self - pos__%w = self%w - pos__%x = self%x - pos__%y = self%y - pos__%z = self%z - + pos__ = [self%w,self%x,self%y,self%z] + end function pos__ @@ -219,26 +210,20 @@ type(quaternion) elemental pure function sub__(self,other) class(quaternion), intent(in) :: self,other - sub__%w = self%w - other%w - sub__%x = self%x - other%x - sub__%y = self%y - other%y - sub__%z = self%z - other%z - + sub__ = [self%w - other%w,self%x - other%x,self%y - other%y,self%z - other%z] + end function sub__ !--------------------------------------------------------------------------------------------------- -!> unary positive operator +!> unary negative operator !--------------------------------------------------------------------------------------------------- type(quaternion) elemental pure function neg__(self) class(quaternion), intent(in) :: self - neg__%w = -self%w - neg__%x = -self%x - neg__%y = -self%y - neg__%z = -self%z - + neg__ = [-self%w,-self%x,-self%y,-self%z] + end function neg__ @@ -258,18 +243,15 @@ end function mul_quat__ !--------------------------------------------------------------------------------------------------- -!> multiplication of quaternions with scalar +!> multiplication of quaternion with scalar !--------------------------------------------------------------------------------------------------- type(quaternion) elemental pure function mul_scal__(self,scal) class(quaternion), intent(in) :: self real(pReal), intent(in) :: scal - mul_scal__%w = self%w*scal - mul_scal__%x = self%x*scal - mul_scal__%y = self%y*scal - mul_scal__%z = self%z*scal - + mul_scal__ = [self%w,self%x,self%y,self%z]*scal + end function mul_scal__ @@ -418,7 +400,7 @@ type(quaternion) elemental pure function conjg__(a) class(quaternion), intent(in) :: a - conjg__ = quaternion([a%w, -a%x, -a%y, -a%z]) + conjg__ = [a%w, -a%x, -a%y, -a%z] end function conjg__ @@ -430,7 +412,7 @@ type(quaternion) elemental pure function quat_homomorphed(self) class(quaternion), intent(in) :: self - quat_homomorphed = quaternion(-[self%w,self%x,self%y,self%z]) + quat_homomorphed = -[self%w,self%x,self%y,self%z] end function quat_homomorphed @@ -444,6 +426,7 @@ pure function asArray(self) class(quaternion), intent(in) :: self asArray = [self%w,self%x,self%y,self%z] + if (self%w < 0) asArray = -asArray end function asArray @@ -484,6 +467,7 @@ subroutine unitTest type(quaternion) :: q, q_2 call random_number(qu) + if (qu(1) < 0.0_pReal) qu = -qu q = qu q_2 = q + q @@ -492,10 +476,10 @@ subroutine unitTest q_2 = q - q if(any(dNeq0(q_2%asArray()))) call IO_error(401,ext_msg='sub__') - q_2 = q * 5.0_preal + q_2 = q * 5.0_pReal if(any(dNeq(q_2%asArray(),5.0_pReal*qu))) call IO_error(401,ext_msg='mul__') - q_2 = q / 0.5_preal + q_2 = q / 0.5_pReal if(any(dNeq(q_2%asArray(),2.0_pReal*qu))) call IO_error(401,ext_msg='div__') q_2 = q From 79cafebffe5ae92c029ab8af987c8401ec3ec8f1 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 11 Jan 2020 03:08:39 +0100 Subject: [PATCH 15/37] following https://www.python.org/dev/peps/pep-0257/ --- src/quaternions.f90 | 106 ++++++++++++++++++++++++-------------------- 1 file changed, 57 insertions(+), 49 deletions(-) diff --git a/src/quaternions.f90 b/src/quaternions.f90 index e55d3804e..3df230e31 100644 --- a/src/quaternions.f90 +++ b/src/quaternions.f90 @@ -31,7 +31,8 @@ !> @author Marc De Graef, Carnegie Mellon University !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @brief general quaternion math, not limited to unit quaternions -!> @details w is the real part, (x, y, z) are the imaginary parts. +!> @details w is the real part, (x, y, z) are the imaginary parts. +!> @details https://users.aalto.fi/~ssarkka/pub/quat.pdf !--------------------------------------------------------------------------------------------------- module quaternions use prec @@ -117,6 +118,14 @@ module quaternions interface log module procedure log__ end interface log + + interface real + module procedure real__ + end interface real + + interface aimag + module procedure aimag__ + end interface aimag private :: & unitTest @@ -125,18 +134,18 @@ contains !-------------------------------------------------------------------------------------------------- -!> @brief doing self test +!> @brief do self test !-------------------------------------------------------------------------------------------------- subroutine quaternions_init - write(6,'(/,a)') ' <<<+- quaternions init -+>>>' + write(6,'(/,a)') ' <<<+- quaternions init -+>>>'; flush(6) call unitTest end subroutine quaternions_init !--------------------------------------------------------------------------------------------------- -!> constructor for a quaternion from a 4-vector +!> construct a quaternion from a 4-vector !--------------------------------------------------------------------------------------------------- type(quaternion) pure function init__(array) @@ -151,7 +160,7 @@ end function init__ !--------------------------------------------------------------------------------------------------- -!> assigning a quaternion +!> assign a quaternion !--------------------------------------------------------------------------------------------------- elemental pure subroutine assign_quat__(self,other) @@ -164,7 +173,7 @@ end subroutine assign_quat__ !--------------------------------------------------------------------------------------------------- -!> assigning a 4-vector +!> assign a 4-vector !--------------------------------------------------------------------------------------------------- pure subroutine assign_vec__(self,other) @@ -180,7 +189,7 @@ end subroutine assign_vec__ !--------------------------------------------------------------------------------------------------- -!> addition of two quaternions +!> add a quaternion !--------------------------------------------------------------------------------------------------- type(quaternion) elemental pure function add__(self,other) @@ -192,7 +201,7 @@ end function add__ !--------------------------------------------------------------------------------------------------- -!> unary positive operator +!> return (unary positive operator) !--------------------------------------------------------------------------------------------------- type(quaternion) elemental pure function pos__(self) @@ -204,7 +213,7 @@ end function pos__ !--------------------------------------------------------------------------------------------------- -!> subtraction of two quaternions +!> subtract a quaternion !--------------------------------------------------------------------------------------------------- type(quaternion) elemental pure function sub__(self,other) @@ -216,7 +225,7 @@ end function sub__ !--------------------------------------------------------------------------------------------------- -!> unary negative operator +!> negate (unary negative operator) !--------------------------------------------------------------------------------------------------- type(quaternion) elemental pure function neg__(self) @@ -228,7 +237,7 @@ end function neg__ !--------------------------------------------------------------------------------------------------- -!> multiplication of two quaternions +!> multiply with a quaternion !--------------------------------------------------------------------------------------------------- type(quaternion) elemental pure function mul_quat__(self,other) @@ -243,7 +252,7 @@ end function mul_quat__ !--------------------------------------------------------------------------------------------------- -!> multiplication of quaternion with scalar +!> multiply with a scalar !--------------------------------------------------------------------------------------------------- type(quaternion) elemental pure function mul_scal__(self,scal) @@ -256,7 +265,7 @@ end function mul_scal__ !--------------------------------------------------------------------------------------------------- -!> division of two quaternions +!> divide by a quaternion !--------------------------------------------------------------------------------------------------- type(quaternion) elemental pure function div_quat__(self,other) @@ -268,7 +277,7 @@ end function div_quat__ !--------------------------------------------------------------------------------------------------- -!> divisiont of quaternions by scalar +!> divide by a scalar !--------------------------------------------------------------------------------------------------- type(quaternion) elemental pure function div_scal__(self,scal) @@ -281,7 +290,7 @@ end function div_scal__ !--------------------------------------------------------------------------------------------------- -!> equality of two quaternions +!> test equality !--------------------------------------------------------------------------------------------------- logical elemental pure function eq__(self,other) @@ -294,7 +303,7 @@ end function eq__ !--------------------------------------------------------------------------------------------------- -!> inequality of two quaternions +!> test inequality !--------------------------------------------------------------------------------------------------- logical elemental pure function neq__(self,other) @@ -306,20 +315,7 @@ end function neq__ !--------------------------------------------------------------------------------------------------- -!> quaternion to the power of a scalar -!--------------------------------------------------------------------------------------------------- -type(quaternion) elemental pure function pow_scal__(self,expon) - - class(quaternion), intent(in) :: self - real(pReal), intent(in) :: expon - - pow_scal__ = exp(log(self)*expon) - -end function pow_scal__ - - -!--------------------------------------------------------------------------------------------------- -!> quaternion to the power of a quaternion +!> raise to the power of a quaternion !--------------------------------------------------------------------------------------------------- type(quaternion) elemental pure function pow_quat__(self,expon) @@ -332,7 +328,20 @@ end function pow_quat__ !--------------------------------------------------------------------------------------------------- -!> exponential of a quaternion +!> raise to the power of a scalar +!--------------------------------------------------------------------------------------------------- +type(quaternion) elemental pure function pow_scal__(self,expon) + + class(quaternion), intent(in) :: self + real(pReal), intent(in) :: expon + + pow_scal__ = exp(log(self)*expon) + +end function pow_scal__ + + +!--------------------------------------------------------------------------------------------------- +!> take exponential !> ToDo: Lacks any check for invalid operations !--------------------------------------------------------------------------------------------------- type(quaternion) elemental pure function exp__(self) @@ -340,7 +349,7 @@ type(quaternion) elemental pure function exp__(self) class(quaternion), intent(in) :: self real(pReal) :: absImag - absImag = norm2([self%x, self%y, self%z]) + absImag = norm2(aimag(self)) exp__ = exp(self%w) * [ cos(absImag), & self%x/absImag * sin(absImag), & @@ -351,7 +360,7 @@ end function exp__ !--------------------------------------------------------------------------------------------------- -!> logarithm of a quaternion +!> take logarithm !> ToDo: Lacks any check for invalid operations !--------------------------------------------------------------------------------------------------- type(quaternion) elemental pure function log__(self) @@ -359,7 +368,7 @@ type(quaternion) elemental pure function log__(self) class(quaternion), intent(in) :: self real(pReal) :: absImag - absImag = norm2([self%x, self%y, self%z]) + absImag = norm2(aimag(self)) log__ = [log(abs(self)), & self%x/absImag * acos(self%w/abs(self)), & @@ -370,7 +379,7 @@ end function log__ !--------------------------------------------------------------------------------------------------- -!> norm of a quaternion +!> return norm !--------------------------------------------------------------------------------------------------- real(pReal) elemental pure function abs__(a) @@ -382,7 +391,7 @@ end function abs__ !--------------------------------------------------------------------------------------------------- -!> dot product of two quaternions +!> calculate dot product !--------------------------------------------------------------------------------------------------- real(pReal) elemental pure function dot_product__(a,b) @@ -394,7 +403,7 @@ end function dot_product__ !--------------------------------------------------------------------------------------------------- -!> conjugate complex of a quaternion +!> take conjugate complex !--------------------------------------------------------------------------------------------------- type(quaternion) elemental pure function conjg__(a) @@ -406,7 +415,7 @@ end function conjg__ !--------------------------------------------------------------------------------------------------- -!> homomorphed quaternion of a quaternion +!> homomorph !--------------------------------------------------------------------------------------------------- type(quaternion) elemental pure function quat_homomorphed(self) @@ -418,7 +427,7 @@ end function quat_homomorphed !--------------------------------------------------------------------------------------------------- -!> quaternion as plain array +!> return as plain array !--------------------------------------------------------------------------------------------------- pure function asArray(self) @@ -432,7 +441,7 @@ end function asArray !--------------------------------------------------------------------------------------------------- -!> real part of a quaternion +!> real part (scalar) !--------------------------------------------------------------------------------------------------- pure function real__(self) @@ -445,7 +454,7 @@ end function real__ !--------------------------------------------------------------------------------------------------- -!> imaginary part of a quaternion +!> imaginary part (3-vector) !--------------------------------------------------------------------------------------------------- pure function aimag__(self) @@ -463,37 +472,36 @@ end function aimag__ subroutine unitTest real(pReal), dimension(4) :: qu - type(quaternion) :: q, q_2 call random_number(qu) if (qu(1) < 0.0_pReal) qu = -qu q = qu - + q_2 = q + q if(any(dNeq(q_2%asArray(),2.0_pReal*qu))) call IO_error(401,ext_msg='add__') - + q_2 = q - q if(any(dNeq0(q_2%asArray()))) call IO_error(401,ext_msg='sub__') - + q_2 = q * 5.0_pReal if(any(dNeq(q_2%asArray(),5.0_pReal*qu))) call IO_error(401,ext_msg='mul__') - + q_2 = q / 0.5_pReal if(any(dNeq(q_2%asArray(),2.0_pReal*qu))) call IO_error(401,ext_msg='div__') - + q_2 = q if(q_2 /= q) call IO_error(401,ext_msg='eq__') if(any(dNeq(q%asArray(),qu))) call IO_error(401,ext_msg='eq__') if(dNeq(q%real(), qu(1))) call IO_error(401,ext_msg='real()') if(any(dNeq(q%aimag(), qu(2:4)))) call IO_error(401,ext_msg='aimag()') - + q_2 = q%homomorphed() if(q /= q_2* (-1.0_pReal)) call IO_error(401,ext_msg='homomorphed') if(dNeq(q_2%real(), qu(1)* (-1.0_pReal))) call IO_error(401,ext_msg='homomorphed/real') if(any(dNeq(q_2%aimag(),qu(2:4)*(-1.0_pReal)))) call IO_error(401,ext_msg='homomorphed/aimag') - + q_2 = conjg(q) if(dNeq(q_2%real(), q%real())) call IO_error(401,ext_msg='conjg/real') if(any(dNeq(q_2%aimag(),q%aimag()*(-1.0_pReal)))) call IO_error(401,ext_msg='conjg/aimag') From aefd401e8cdda0ad4164f09728595f4c2467d7b1 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 11 Jan 2020 03:11:45 +0100 Subject: [PATCH 16/37] this is a quaternion class it is meant to represent any quaternion, not only unit quaternions/rotations that follow a specific convention. Need to check in rotations.f90 where the homomorph should happen --- src/quaternions.f90 | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/quaternions.f90 b/src/quaternions.f90 index 3df230e31..9a14fa211 100644 --- a/src/quaternions.f90 +++ b/src/quaternions.f90 @@ -435,7 +435,6 @@ pure function asArray(self) class(quaternion), intent(in) :: self asArray = [self%w,self%x,self%y,self%z] - if (self%w < 0) asArray = -asArray end function asArray @@ -475,7 +474,7 @@ subroutine unitTest type(quaternion) :: q, q_2 call random_number(qu) - if (qu(1) < 0.0_pReal) qu = -qu + qu = (qu-0.5_pReal) * 2.0_pReal q = qu q_2 = q + q From c7180c329571a952b37bca68953910c82588ef79 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 11 Jan 2020 03:50:17 +0100 Subject: [PATCH 17/37] some more tests for quaternion operations --- src/quaternions.f90 | 17 +++++++++++++++-- 1 file changed, 15 insertions(+), 2 deletions(-) diff --git a/src/quaternions.f90 b/src/quaternions.f90 index 9a14fa211..252135eec 100644 --- a/src/quaternions.f90 +++ b/src/quaternions.f90 @@ -475,7 +475,10 @@ subroutine unitTest call random_number(qu) qu = (qu-0.5_pReal) * 2.0_pReal - q = qu + q = quaternion(qu) + + q_2= qu + if(any(dNeq(q%asArray(),q_2%asArray()))) call IO_error(401,ext_msg='assign_vec__') q_2 = q + q if(any(dNeq(q_2%asArray(),2.0_pReal*qu))) call IO_error(401,ext_msg='add__') @@ -489,8 +492,13 @@ subroutine unitTest q_2 = q / 0.5_pReal if(any(dNeq(q_2%asArray(),2.0_pReal*qu))) call IO_error(401,ext_msg='div__') + q_2 = q * 0.3_pReal + if(dNeq0(abs(q)) .and. q_2 == q) call IO_error(401,ext_msg='eq__') + q_2 = q - if(q_2 /= q) call IO_error(401,ext_msg='eq__') + if(q_2 /= q) call IO_error(401,ext_msg='neq__') + + if(dNeq(abs(q),norm2(qu))) call IO_error(401,ext_msg='abs__') if(any(dNeq(q%asArray(),qu))) call IO_error(401,ext_msg='eq__') if(dNeq(q%real(), qu(1))) call IO_error(401,ext_msg='real()') @@ -504,6 +512,11 @@ subroutine unitTest q_2 = conjg(q) if(dNeq(q_2%real(), q%real())) call IO_error(401,ext_msg='conjg/real') if(any(dNeq(q_2%aimag(),q%aimag()*(-1.0_pReal)))) call IO_error(401,ext_msg='conjg/aimag') + + if (norm2(aimag(q)) * abs(real(q)) > 0.0_pReal) then + if (dNeq0(abs(q-exp(log(q))),1.0e-12_pReal)) call IO_error(401,ext_msg='exp/log') + if (dNeq0(abs(q-log(exp(q))),1.0e-12_pReal)) call IO_error(401,ext_msg='log/exp') + endif end subroutine unitTest From 115716b8c21df19868bc5f52c53c0ba2ce38724f Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 11 Jan 2020 03:58:12 +0100 Subject: [PATCH 18/37] polishing/use existing functions --- src/quaternions.f90 | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/src/quaternions.f90 b/src/quaternions.f90 index 252135eec..a49f53007 100644 --- a/src/quaternions.f90 +++ b/src/quaternions.f90 @@ -195,7 +195,8 @@ type(quaternion) elemental pure function add__(self,other) class(quaternion), intent(in) :: self,other - add__ = [self%w + other%w,self%x + other%x,self%y + other%y,self%z + other%z] + add__ = [ self%w, self%x, self%y ,self%z] & + + [other%w, other%x, other%y,other%z] end function add__ @@ -207,7 +208,7 @@ type(quaternion) elemental pure function pos__(self) class(quaternion), intent(in) :: self - pos__ = [self%w,self%x,self%y,self%z] + pos__ = self * (+1.0_pReal) end function pos__ @@ -219,7 +220,8 @@ type(quaternion) elemental pure function sub__(self,other) class(quaternion), intent(in) :: self,other - sub__ = [self%w - other%w,self%x - other%x,self%y - other%y,self%z - other%z] + sub__ = [ self%w, self%x, self%y ,self%z] & + - [other%w, other%x, other%y,other%z] end function sub__ @@ -231,7 +233,7 @@ type(quaternion) elemental pure function neg__(self) class(quaternion), intent(in) :: self - neg__ = [-self%w,-self%x,-self%y,-self%z] + neg__ = self * (-1.0_pReal) end function neg__ @@ -333,7 +335,7 @@ end function pow_quat__ type(quaternion) elemental pure function pow_scal__(self,expon) class(quaternion), intent(in) :: self - real(pReal), intent(in) :: expon + real(pReal), intent(in) :: expon pow_scal__ = exp(log(self)*expon) @@ -421,7 +423,7 @@ type(quaternion) elemental pure function quat_homomorphed(self) class(quaternion), intent(in) :: self - quat_homomorphed = -[self%w,self%x,self%y,self%z] + quat_homomorphed = - self end function quat_homomorphed From de95ca590608c150f081e8fa4d9f247feb188118 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 11 Jan 2020 04:15:51 +0100 Subject: [PATCH 19/37] inverse of a quaternion --- src/quaternions.f90 | 38 ++++++++++++++++++++++++++++---------- 1 file changed, 28 insertions(+), 10 deletions(-) diff --git a/src/quaternions.f90 b/src/quaternions.f90 index a49f53007..0ce7765e6 100644 --- a/src/quaternions.f90 +++ b/src/quaternions.f90 @@ -82,12 +82,13 @@ module quaternions procedure, public :: conjg__ procedure, public :: exp__ procedure, public :: log__ - - procedure, public :: homomorphed => quat_homomorphed - procedure, public :: asArray procedure, public :: real => real__ procedure, public :: aimag => aimag__ + procedure, public :: homomorphed + procedure, public :: asArray + procedure, public :: inverse + end type interface assignment (=) @@ -137,7 +138,6 @@ contains !> @brief do self test !-------------------------------------------------------------------------------------------------- subroutine quaternions_init - write(6,'(/,a)') ' <<<+- quaternions init -+>>>'; flush(6) call unitTest @@ -419,13 +419,13 @@ end function conjg__ !--------------------------------------------------------------------------------------------------- !> homomorph !--------------------------------------------------------------------------------------------------- -type(quaternion) elemental pure function quat_homomorphed(self) +type(quaternion) elemental pure function homomorphed(self) class(quaternion), intent(in) :: self - quat_homomorphed = - self + homomorphed = - self -end function quat_homomorphed +end function homomorphed !--------------------------------------------------------------------------------------------------- @@ -467,6 +467,18 @@ pure function aimag__(self) end function aimag__ +!--------------------------------------------------------------------------------------------------- +!> inverse +!--------------------------------------------------------------------------------------------------- +type(quaternion) elemental pure function inverse(self) + + class(quaternion), intent(in) :: self + + inverse = conjg(self)/abs(self)**2.0_pReal + +end function inverse + + !-------------------------------------------------------------------------------------------------- !> @brief check correctness of (some) quaternions functions !-------------------------------------------------------------------------------------------------- @@ -478,7 +490,7 @@ subroutine unitTest call random_number(qu) qu = (qu-0.5_pReal) * 2.0_pReal q = quaternion(qu) - + q_2= qu if(any(dNeq(q%asArray(),q_2%asArray()))) call IO_error(401,ext_msg='assign_vec__') @@ -515,9 +527,15 @@ subroutine unitTest if(dNeq(q_2%real(), q%real())) call IO_error(401,ext_msg='conjg/real') if(any(dNeq(q_2%aimag(),q%aimag()*(-1.0_pReal)))) call IO_error(401,ext_msg='conjg/aimag') + if(abs(q) > 0.0_pReal) then + q_2 = q * q%inverse() + if( dNeq(real(q_2), 1.0_pReal,1.0e-15_pReal)) call IO_error(401,ext_msg='inverse/real') + if(any(dNeq0(aimag(q_2), 1.0e-15_pReal))) call IO_error(401,ext_msg='inverse/aimag') + endif + if (norm2(aimag(q)) * abs(real(q)) > 0.0_pReal) then - if (dNeq0(abs(q-exp(log(q))),1.0e-12_pReal)) call IO_error(401,ext_msg='exp/log') - if (dNeq0(abs(q-log(exp(q))),1.0e-12_pReal)) call IO_error(401,ext_msg='log/exp') + if (dNeq0(abs(q-exp(log(q))),1.0e-13_pReal)) call IO_error(401,ext_msg='exp/log') + if (dNeq0(abs(q-log(exp(q))),1.0e-13_pReal)) call IO_error(401,ext_msg='log/exp') endif end subroutine unitTest From f02851959765844654bd43265cfd34b1c5cc6a6a Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 11 Jan 2020 04:44:30 +0100 Subject: [PATCH 20/37] some facts from wikipedia as tests --- src/quaternions.f90 | 50 ++++++++++++++++++++++++++------------------- 1 file changed, 29 insertions(+), 21 deletions(-) diff --git a/src/quaternions.f90 b/src/quaternions.f90 index 0ce7765e6..d474f3994 100644 --- a/src/quaternions.f90 +++ b/src/quaternions.f90 @@ -79,9 +79,9 @@ module quaternions procedure, public :: abs__ procedure, public :: dot_product__ - procedure, public :: conjg__ procedure, public :: exp__ procedure, public :: log__ + procedure, public :: conjg => conjg__ procedure, public :: real => real__ procedure, public :: aimag => aimag__ @@ -138,6 +138,7 @@ contains !> @brief do self test !-------------------------------------------------------------------------------------------------- subroutine quaternions_init + write(6,'(/,a)') ' <<<+- quaternions init -+>>>'; flush(6) call unitTest @@ -259,7 +260,7 @@ end function mul_quat__ type(quaternion) elemental pure function mul_scal__(self,scal) class(quaternion), intent(in) :: self - real(pReal), intent(in) :: scal + real(pReal), intent(in) :: scal mul_scal__ = [self%w,self%x,self%y,self%z]*scal @@ -284,7 +285,7 @@ end function div_quat__ type(quaternion) elemental pure function div_scal__(self,scal) class(quaternion), intent(in) :: self - real(pReal), intent(in) :: scal + real(pReal), intent(in) :: scal div_scal__ = [self%w,self%x,self%y,self%z]/scal @@ -492,50 +493,57 @@ subroutine unitTest q = quaternion(qu) q_2= qu - if(any(dNeq(q%asArray(),q_2%asArray()))) call IO_error(401,ext_msg='assign_vec__') + if(any(dNeq(q%asArray(),q_2%asArray()))) call IO_error(401,ext_msg='assign_vec__') q_2 = q + q - if(any(dNeq(q_2%asArray(),2.0_pReal*qu))) call IO_error(401,ext_msg='add__') + if(any(dNeq(q_2%asArray(),2.0_pReal*qu))) call IO_error(401,ext_msg='add__') q_2 = q - q - if(any(dNeq0(q_2%asArray()))) call IO_error(401,ext_msg='sub__') + if(any(dNeq0(q_2%asArray()))) call IO_error(401,ext_msg='sub__') q_2 = q * 5.0_pReal - if(any(dNeq(q_2%asArray(),5.0_pReal*qu))) call IO_error(401,ext_msg='mul__') + if(any(dNeq(q_2%asArray(),5.0_pReal*qu))) call IO_error(401,ext_msg='mul__') q_2 = q / 0.5_pReal - if(any(dNeq(q_2%asArray(),2.0_pReal*qu))) call IO_error(401,ext_msg='div__') + if(any(dNeq(q_2%asArray(),2.0_pReal*qu))) call IO_error(401,ext_msg='div__') q_2 = q * 0.3_pReal - if(dNeq0(abs(q)) .and. q_2 == q) call IO_error(401,ext_msg='eq__') + if(dNeq0(abs(q)) .and. q_2 == q) call IO_error(401,ext_msg='eq__') q_2 = q - if(q_2 /= q) call IO_error(401,ext_msg='neq__') + if(q_2 /= q) call IO_error(401,ext_msg='neq__') - if(dNeq(abs(q),norm2(qu))) call IO_error(401,ext_msg='abs__') + if(dNeq(abs(q),norm2(qu))) call IO_error(401,ext_msg='abs__') + if(dNeq(abs(q)**2.0_pReal, real(q*q%conjg()))) call IO_error(401,ext_msg='abs__/*conjg') - if(any(dNeq(q%asArray(),qu))) call IO_error(401,ext_msg='eq__') - if(dNeq(q%real(), qu(1))) call IO_error(401,ext_msg='real()') - if(any(dNeq(q%aimag(), qu(2:4)))) call IO_error(401,ext_msg='aimag()') + if(any(dNeq(q%asArray(),qu))) call IO_error(401,ext_msg='eq__') + if(dNeq(q%real(), qu(1))) call IO_error(401,ext_msg='real()') + if(any(dNeq(q%aimag(), qu(2:4)))) call IO_error(401,ext_msg='aimag()') q_2 = q%homomorphed() - if(q /= q_2* (-1.0_pReal)) call IO_error(401,ext_msg='homomorphed') - if(dNeq(q_2%real(), qu(1)* (-1.0_pReal))) call IO_error(401,ext_msg='homomorphed/real') - if(any(dNeq(q_2%aimag(),qu(2:4)*(-1.0_pReal)))) call IO_error(401,ext_msg='homomorphed/aimag') + if(q /= q_2* (-1.0_pReal)) call IO_error(401,ext_msg='homomorphed') + if(dNeq(q_2%real(), qu(1)* (-1.0_pReal))) call IO_error(401,ext_msg='homomorphed/real') + if(any(dNeq(q_2%aimag(),qu(2:4)*(-1.0_pReal)))) call IO_error(401,ext_msg='homomorphed/aimag') q_2 = conjg(q) - if(dNeq(q_2%real(), q%real())) call IO_error(401,ext_msg='conjg/real') - if(any(dNeq(q_2%aimag(),q%aimag()*(-1.0_pReal)))) call IO_error(401,ext_msg='conjg/aimag') + if(dNeq(abs(q),abs(q_2))) call IO_error(401,ext_msg='conjg/abs') + if(q /= conjg(q_2)) call IO_error(401,ext_msg='conjg/involution') + if(dNeq(q_2%real(), q%real())) call IO_error(401,ext_msg='conjg/real') + if(any(dNeq(q_2%aimag(),q%aimag()*(-1.0_pReal)))) call IO_error(401,ext_msg='conjg/aimag') if(abs(q) > 0.0_pReal) then q_2 = q * q%inverse() if( dNeq(real(q_2), 1.0_pReal,1.0e-15_pReal)) call IO_error(401,ext_msg='inverse/real') if(any(dNeq0(aimag(q_2), 1.0e-15_pReal))) call IO_error(401,ext_msg='inverse/aimag') + + q_2 = q/abs(q) + q_2 = conjg(q_2) - inverse(q_2) + if(any(dNeq0(q_2%asArray(),1.0e-15_pReal))) call IO_error(401,ext_msg='inverse/conjg') endif if (norm2(aimag(q)) * abs(real(q)) > 0.0_pReal) then - if (dNeq0(abs(q-exp(log(q))),1.0e-13_pReal)) call IO_error(401,ext_msg='exp/log') - if (dNeq0(abs(q-log(exp(q))),1.0e-13_pReal)) call IO_error(401,ext_msg='log/exp') + if (dNeq0(abs(q-exp(log(q))),1.0e-13_pReal)) call IO_error(401,ext_msg='exp/log') + if (dNeq0(abs(q-log(exp(q))),1.0e-13_pReal)) call IO_error(401,ext_msg='log/exp') endif end subroutine unitTest From 3a6819f548897d20a04a1df5498b764ee79d5413 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 11 Jan 2020 05:14:17 +0100 Subject: [PATCH 21/37] check for invalid operations --- src/quaternions.f90 | 24 +++++++++++++----------- 1 file changed, 13 insertions(+), 11 deletions(-) diff --git a/src/quaternions.f90 b/src/quaternions.f90 index d474f3994..bdd44b93b 100644 --- a/src/quaternions.f90 +++ b/src/quaternions.f90 @@ -345,7 +345,6 @@ end function pow_scal__ !--------------------------------------------------------------------------------------------------- !> take exponential -!> ToDo: Lacks any check for invalid operations !--------------------------------------------------------------------------------------------------- type(quaternion) elemental pure function exp__(self) @@ -354,17 +353,18 @@ type(quaternion) elemental pure function exp__(self) absImag = norm2(aimag(self)) - exp__ = exp(self%w) * [ cos(absImag), & - self%x/absImag * sin(absImag), & - self%y/absImag * sin(absImag), & - self%z/absImag * sin(absImag)] + exp__ = merge(exp(self%w) * [ cos(absImag), & + self%x/absImag * sin(absImag), & + self%y/absImag * sin(absImag), & + self%z/absImag * sin(absImag)], & + IEEE_value(1.0_pReal,IEEE_SIGNALING_NAN), & + dNeq0(absImag)) end function exp__ !--------------------------------------------------------------------------------------------------- !> take logarithm -!> ToDo: Lacks any check for invalid operations !--------------------------------------------------------------------------------------------------- type(quaternion) elemental pure function log__(self) @@ -373,10 +373,12 @@ type(quaternion) elemental pure function log__(self) absImag = norm2(aimag(self)) - log__ = [log(abs(self)), & - self%x/absImag * acos(self%w/abs(self)), & - self%y/absImag * acos(self%w/abs(self)), & - self%z/absImag * acos(self%w/abs(self))] + log__ = merge([log(abs(self)), & + self%x/absImag * acos(self%w/abs(self)), & + self%y/absImag * acos(self%w/abs(self)), & + self%z/absImag * acos(self%w/abs(self))], & + IEEE_value(1.0_pReal,IEEE_SIGNALING_NAN), & + dNeq0(absImag)) end function log__ @@ -541,7 +543,7 @@ subroutine unitTest if(any(dNeq0(q_2%asArray(),1.0e-15_pReal))) call IO_error(401,ext_msg='inverse/conjg') endif - if (norm2(aimag(q)) * abs(real(q)) > 0.0_pReal) then + if (norm2(aimag(q)) > 0.0_pReal) then if (dNeq0(abs(q-exp(log(q))),1.0e-13_pReal)) call IO_error(401,ext_msg='exp/log') if (dNeq0(abs(q-log(exp(q))),1.0e-13_pReal)) call IO_error(401,ext_msg='log/exp') endif From 842666cc20a4aafd1608f22e70ab1b152c1f9396 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 11 Jan 2020 11:23:53 +0100 Subject: [PATCH 22/37] no overlap with Marc's code --- src/quaternions.f90 | 33 ++------------------------------- 1 file changed, 2 insertions(+), 31 deletions(-) diff --git a/src/quaternions.f90 b/src/quaternions.f90 index bdd44b93b..886c7d071 100644 --- a/src/quaternions.f90 +++ b/src/quaternions.f90 @@ -1,38 +1,9 @@ -! ################################################################### -! Copyright (c) 2013-2015, Marc De Graef/Carnegie Mellon University -! Modified 2017-2019, Martin Diehl/Max-Planck-Institut für Eisenforschung GmbH -! All rights reserved. -! -! Redistribution and use in source and binary forms, with or without modification, are -! permitted provided that the following conditions are met: -! -! - Redistributions of source code must retain the above copyright notice, this list -! of conditions and the following disclaimer. -! - Redistributions in binary form must reproduce the above copyright notice, this -! list of conditions and the following disclaimer in the documentation and/or -! other materials provided with the distribution. -! - Neither the names of Marc De Graef, Carnegie Mellon University nor the names -! of its contributors may be used to endorse or promote products derived from -! this software without specific prior written permission. -! -! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" -! AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE -! IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE -! ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE -! LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL -! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR -! SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER -! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, -! OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE -! USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -! ################################################################### - !--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH +!> @author Philip Eisenlohr, Michigan State University !> @brief general quaternion math, not limited to unit quaternions !> @details w is the real part, (x, y, z) are the imaginary parts. -!> @details https://users.aalto.fi/~ssarkka/pub/quat.pdf +!> @details https://en.wikipedia.org/wiki/Quaternion !--------------------------------------------------------------------------------------------------- module quaternions use prec From e762cb4dfd5352ac18b2d7f53c926dacac1ad779 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 11 Jan 2020 12:36:35 +0100 Subject: [PATCH 23/37] issue with gfortran < 9 the false branch of merge seems to be evaluated which results in a signaling NaN --- src/quaternions.f90 | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/quaternions.f90 b/src/quaternions.f90 index 886c7d071..a67f345a0 100644 --- a/src/quaternions.f90 +++ b/src/quaternions.f90 @@ -513,12 +513,14 @@ subroutine unitTest q_2 = conjg(q_2) - inverse(q_2) if(any(dNeq0(q_2%asArray(),1.0e-15_pReal))) call IO_error(401,ext_msg='inverse/conjg') endif - + +#if !(defined(__GFORTRAN__) && __GNUC__ < 9) if (norm2(aimag(q)) > 0.0_pReal) then if (dNeq0(abs(q-exp(log(q))),1.0e-13_pReal)) call IO_error(401,ext_msg='exp/log') if (dNeq0(abs(q-log(exp(q))),1.0e-13_pReal)) call IO_error(401,ext_msg='log/exp') endif - +#endif + end subroutine unitTest From ac112d2d36f783b942186c5d6fa003bb62f194d3 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 11 Jan 2020 13:55:56 +0100 Subject: [PATCH 24/37] tolerance needed for optimized code --- src/quaternions.f90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/quaternions.f90 b/src/quaternions.f90 index a67f345a0..0ca404880 100644 --- a/src/quaternions.f90 +++ b/src/quaternions.f90 @@ -487,7 +487,8 @@ subroutine unitTest if(q_2 /= q) call IO_error(401,ext_msg='neq__') if(dNeq(abs(q),norm2(qu))) call IO_error(401,ext_msg='abs__') - if(dNeq(abs(q)**2.0_pReal, real(q*q%conjg()))) call IO_error(401,ext_msg='abs__/*conjg') + if(dNeq(abs(q)**2.0_pReal, real(q*q%conjg()),1.0e-14_pReal)) & + call IO_error(401,ext_msg='abs__/*conjg') if(any(dNeq(q%asArray(),qu))) call IO_error(401,ext_msg='eq__') if(dNeq(q%real(), qu(1))) call IO_error(401,ext_msg='real()') From 300f1b7015cf644e4c175aea43483d651493bff8 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Sat, 11 Jan 2020 11:36:22 -0500 Subject: [PATCH 25/37] added options to return "natural" versions of asQ, asRodrig, and asAxisAngle --- python/damask/orientation.py | 40 +++++++++++++++++++++++------------- 1 file changed, 26 insertions(+), 14 deletions(-) diff --git a/python/damask/orientation.py b/python/damask/orientation.py index a63444155..a86ba9331 100644 --- a/python/damask/orientation.py +++ b/python/damask/orientation.py @@ -170,9 +170,18 @@ class Rotation: ################################################################################################ # convert to different orientation representations (numpy arrays) - def asQuaternion(self): - """Unit quaternion: (q, p_1, p_2, p_3).""" - return self.quaternion.asArray() + def asQuaternion(self, + quaternion = False): + """ + Unit quaternion [q, p_1, p_2, p_3] unless quaternion == True: damask.quaternion object. + + Parameters + ---------- + quaternion : bool, optional + return quaternion as DAMASK object. + + """ + return self.quaternion if quaternion else self.quaternion.asArray() def asEulers(self, degrees = False): @@ -190,33 +199,36 @@ class Rotation: return eu def asAxisAngle(self, - degrees = False): + degrees = False, + pair = False): """ - Axis angle pair: ([n_1, n_2, n_3], ω). + Axis angle representation [n_1, n_2, n_3, ω] unless pair == True: ([n_1, n_2, n_3], ω). Parameters ---------- degrees : bool, optional return rotation angle in degrees. + pair : bool, optional + return tuple of axis and angle. """ ax = qu2ax(self.quaternion.asArray()) if degrees: ax[3] = np.degrees(ax[3]) - return ax + return (ax[:3],np.degrees(ax[3])) if pair else ax def asMatrix(self): """Rotation matrix.""" return qu2om(self.quaternion.asArray()) def asRodrigues(self, - vector=False): + vector = False): """ - Rodrigues-Frank vector: ([n_1, n_2, n_3], tan(ω/2)). + Rodrigues-Frank vector representation [n_1, n_2, n_3, tan(ω/2)] unless vector == True: [n_1, n_2, n_3] * tan(ω/2). Parameters ---------- vector : bool, optional - return as array of length 3, i.e. scale the unit vector giving the rotation axis. + return as actual Rodrigues--Frank vector, i.e. rotation axis scaled by tan(ω/2). """ ro = qu2ro(self.quaternion.asArray()) @@ -252,8 +264,8 @@ class Rotation: acceptHomomorph = False, P = -1): - qu = quaternion if isinstance(quaternion, np.ndarray) and quaternion.dtype == np.dtype(float) \ - else np.array(quaternion,dtype=float) + qu = quaternion if isinstance(quaternion, np.ndarray) and quaternion.dtype == np.dtype(float) \ + else np.array(quaternion,dtype=float) if P > 0: qu[1:4] *= -1 # convert from P=1 to P=-1 if qu[0] < 0.0: if acceptHomomorph: @@ -1193,9 +1205,9 @@ class Orientation: ref = orientations[0] for o in orientations: closest.append(o.equivalentOrientations( - ref.disorientation(o, - SST = False, # select (o[ther]'s) sym orientation - symmetries = True)[2]).rotation) # with lowest misorientation + ref.disorientation(o, + SST = False, # select (o[ther]'s) sym orientation + symmetries = True)[2]).rotation) # with lowest misorientation return Orientation(Rotation.fromAverage(closest,weights),ref.lattice) From d4535dadb4ea4cf6973fc5f9964bd49fbf33dde3 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 11 Jan 2020 20:33:29 +0100 Subject: [PATCH 26/37] use American english --- processing/post/addCompatibilityMismatch.py | 4 ++-- src/mesh_grid.f90 | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/processing/post/addCompatibilityMismatch.py b/processing/post/addCompatibilityMismatch.py index 8aff13c53..29e874614 100755 --- a/processing/post/addCompatibilityMismatch.py +++ b/processing/post/addCompatibilityMismatch.py @@ -182,8 +182,8 @@ for name in filenames: nodes = damask.grid_filters.node_coord(size,F) if options.shape: - centres = damask.grid_filters.cell_coord(size,F) - shapeMismatch = shapeMismatch( size,table.get(options.defgrad).reshape(grid[2],grid[1],grid[0],3,3),nodes,centres) + centers = damask.grid_filters.cell_coord(size,F) + shapeMismatch = shapeMismatch( size,table.get(options.defgrad).reshape(grid[2],grid[1],grid[0],3,3),nodes,centers) table.add('shapeMismatch(({}))'.format(options.defgrad), shapeMismatch.reshape((-1,1)), scriptID+' '+' '.join(sys.argv[1:])) diff --git a/src/mesh_grid.f90 b/src/mesh_grid.f90 index 3d839e1c9..62ce67397 100644 --- a/src/mesh_grid.f90 +++ b/src/mesh_grid.f90 @@ -296,7 +296,7 @@ end subroutine readGeom !--------------------------------------------------------------------------------------------------- -!> @brief Calculate undeformed position of IPs/cell centres (pretend to be an element) +!> @brief Calculate undeformed position of IPs/cell centers (pretend to be an element) !--------------------------------------------------------------------------------------------------- function IPcoordinates0(grid,geomSize,grid3Offset) From 8f99dd85ac5ec48ba6152e152d234228e491e212 Mon Sep 17 00:00:00 2001 From: Test User Date: Sat, 11 Jan 2020 22:16:08 +0100 Subject: [PATCH 27/37] [skip ci] updated version information after successful test of v2.0.3-1308-g612e8e34 --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index 287da9b11..cf29053c3 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.3-1294-g034367fa +v2.0.3-1308-g612e8e34 From d9afdaed165984d5a03d376ae37494c4bd61d900 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 11 Jan 2020 23:43:03 +0100 Subject: [PATCH 28/37] using updated tests (no crystallite anymore) --- PRIVATE | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PRIVATE b/PRIVATE index cda69f9a5..3eeffba15 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit cda69f9a59fe64223439a2c725e1a78cf22b28aa +Subproject commit 3eeffba1513c0495032d0f33522e3aaa1ed77465 From 97ddd56540d82bdad9c12a1ba6c166423928b5fb Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 11 Jan 2020 23:43:31 +0100 Subject: [PATCH 29/37] avoid strange deformation gradients otherwise, the test fails in a few cases (determinants were in the range 1e-5 to 1e16) --- python/tests/test_mechanics.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/tests/test_mechanics.py b/python/tests/test_mechanics.py index 1ea1f2bba..9e1d9bc0c 100644 --- a/python/tests/test_mechanics.py +++ b/python/tests/test_mechanics.py @@ -113,7 +113,7 @@ class TestMechanics: def test_strain_tensor_rotation_equivalence(self): """Ensure that left and right strain differ only by a rotation.""" - F = np.random.random((self.n,3,3)) + F = np.broadcast_to(np.eye(3),[self.n,3,3]) + (np.random.random((self.n,3,3))*0.5 - 0.25) m = np.random.random()*5.0-2.5 assert np.allclose(np.linalg.det(mechanics.strain_tensor(F,'U',m)), np.linalg.det(mechanics.strain_tensor(F,'V',m))) From ddd8027b8a354c0ad7535ac10bdaa96e41599cf3 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 12 Jan 2020 00:10:42 +0100 Subject: [PATCH 30/37] autodetect string length --- src/FEsolving.f90 | 4 ++-- src/IO.f90 | 8 ++++---- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/FEsolving.f90 b/src/FEsolving.f90 index 38788a065..8c9683266 100644 --- a/src/FEsolving.f90 +++ b/src/FEsolving.f90 @@ -46,10 +46,10 @@ subroutine FE_init call IO_open_inputFile(FILEUNIT) rewind(FILEUNIT) do - read (FILEUNIT,'(a256)',END=100) line + read (FILEUNIT,'(A)',END=100) line chunkPos = IO_stringPos(line) if(IO_lc(IO_stringValue(line,chunkPos,1)) == 'solver') then - read (FILEUNIT,'(a256)',END=100) line ! next line + read (FILEUNIT,'(A)',END=100) line ! next line chunkPos = IO_stringPos(line) symmetricSolver = (IO_intValue(line,chunkPos,2) /= 1) endif diff --git a/src/IO.f90 b/src/IO.f90 index 1e39daa1e..0436e9b19 100644 --- a/src/IO.f90 +++ b/src/IO.f90 @@ -244,7 +244,7 @@ subroutine IO_open_inputFile(fileUnit) do - read(unit2,'(A256)',END=220) line + read(unit2,'(A)',END=220) line chunkPos = IO_stringPos(line) if (IO_lc(IO_StringValue(line,chunkPos,1))=='*include') then @@ -1000,7 +1000,7 @@ function IO_continuousIntValues(fileUnit,maxN,lookupName,lookupMap,lookupMaxN) #if defined(Marc4DAMASK) do - read(fileUnit,'(A256)',end=100) line + read(fileUnit,'(A)',end=100) line chunkPos = IO_stringPos(line) if (chunkPos(1) < 1) then ! empty line exit @@ -1041,14 +1041,14 @@ function IO_continuousIntValues(fileUnit,maxN,lookupName,lookupMap,lookupMaxN) !-------------------------------------------------------------------------------------------------- ! check if the element values in the elset are auto generated backspace(fileUnit) - read(fileUnit,'(A256)',end=100) line + read(fileUnit,'(A)',end=100) line chunkPos = IO_stringPos(line) do i = 1,chunkPos(1) if (IO_lc(IO_stringValue(line,chunkPos,i)) == 'generate') rangeGeneration = .true. enddo do l = 1,c - read(fileUnit,'(A256)',end=100) line + read(fileUnit,'(A)',end=100) line chunkPos = IO_stringPos(line) if (verify(IO_stringValue(line,chunkPos,1),'0123456789') > 0) then ! a non-int, i.e. set names follow on this line do i = 1,chunkPos(1) ! loop over set names in line From 1bb94f03b8f15387918439deae6c79ee57ebda1b Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 12 Jan 2020 00:14:35 +0100 Subject: [PATCH 31/37] polishing (prospector was complaining) --- processing/post/addStrainTensors.py | 3 +-- processing/post/addTable.py | 1 + processing/post/averageDown.py | 1 + processing/post/blowUp.py | 1 + processing/post/growTable.py | 3 +-- python/damask/table.py | 9 ++++++--- 6 files changed, 11 insertions(+), 7 deletions(-) diff --git a/processing/post/addStrainTensors.py b/processing/post/addStrainTensors.py index 756761e46..77015a91f 100755 --- a/processing/post/addStrainTensors.py +++ b/processing/post/addStrainTensors.py @@ -2,10 +2,9 @@ import os import sys +from io import StringIO from optparse import OptionParser -import numpy as np - import damask diff --git a/processing/post/addTable.py b/processing/post/addTable.py index c7d34f1e8..944214e69 100755 --- a/processing/post/addTable.py +++ b/processing/post/addTable.py @@ -2,6 +2,7 @@ import os import sys +from io import StringIO from optparse import OptionParser import damask diff --git a/processing/post/averageDown.py b/processing/post/averageDown.py index 817000dfa..cbd1f9637 100755 --- a/processing/post/averageDown.py +++ b/processing/post/averageDown.py @@ -2,6 +2,7 @@ import os import sys +from io import StringIO from optparse import OptionParser import numpy as np diff --git a/processing/post/blowUp.py b/processing/post/blowUp.py index ec42e2365..9cb6347ab 100755 --- a/processing/post/blowUp.py +++ b/processing/post/blowUp.py @@ -2,6 +2,7 @@ import os import sys +from io import StringIO from optparse import OptionParser from scipy import ndimage diff --git a/processing/post/growTable.py b/processing/post/growTable.py index 4fd647b1e..1dbfa8423 100755 --- a/processing/post/growTable.py +++ b/processing/post/growTable.py @@ -2,10 +2,9 @@ import os import sys +from io import StringIO from optparse import OptionParser -import numpy as np - import damask scriptName = os.path.splitext(os.path.basename(__file__))[0] diff --git a/python/damask/table.py b/python/damask/table.py index ef8a84276..343ba638d 100644 --- a/python/damask/table.py +++ b/python/damask/table.py @@ -274,7 +274,9 @@ class Table(): def append(self,other): """ - Append other table vertically (similar to numpy.vstack). Requires matching shapes and order. + Append other table vertically (similar to numpy.vstack). + + Requires matching labels/shapes and order. Parameters ---------- @@ -290,8 +292,9 @@ class Table(): def join(self,other): """ - Append other table horizontally (similar to numpy.hstack). Requires matching number of rows - and no common lables + Append other table horizontally (similar to numpy.hstack). + + Requires matching number of rows and no common labels. Parameters ---------- From 5bc1c98da7f520b218da6014f22d64442ad3e2e6 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 12 Jan 2020 00:49:03 +0100 Subject: [PATCH 32/37] use 0-based indexing for worldrank --- src/grid/grid_damage_spectral.f90 | 6 +++--- src/grid/grid_mech_FEM.f90 | 6 +++--- src/grid/grid_mech_spectral_polarisation.f90 | 6 +++--- src/grid/grid_thermal_spectral.f90 | 6 +++--- 4 files changed, 12 insertions(+), 12 deletions(-) diff --git a/src/grid/grid_damage_spectral.f90 b/src/grid/grid_damage_spectral.f90 index e2ce28c0e..97b20a6db 100644 --- a/src/grid/grid_damage_spectral.f90 +++ b/src/grid/grid_damage_spectral.f90 @@ -54,7 +54,7 @@ contains !-------------------------------------------------------------------------------------------------- subroutine grid_damage_spectral_init - PetscInt, dimension(worldsize) :: localK + PetscInt, dimension(0:worldsize-1) :: localK integer :: i, j, k, cell DM :: damage_grid Vec :: uBound, lBound @@ -78,8 +78,8 @@ subroutine grid_damage_spectral_init ! initialize solver specific parts of PETSc call SNESCreate(PETSC_COMM_WORLD,damage_snes,ierr); CHKERRQ(ierr) call SNESSetOptionsPrefix(damage_snes,'damage_',ierr);CHKERRQ(ierr) - localK = 0 - localK(worldrank+1) = grid3 + localK = 0 + localK(worldrank) = grid3 call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,PETSC_COMM_WORLD,ierr) call DMDACreate3D(PETSC_COMM_WORLD, & DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & ! cut off stencil at boundary diff --git a/src/grid/grid_mech_FEM.f90 b/src/grid/grid_mech_FEM.f90 index 05fc3520a..b51ebb56a 100644 --- a/src/grid/grid_mech_FEM.f90 +++ b/src/grid/grid_mech_FEM.f90 @@ -81,7 +81,7 @@ contains subroutine grid_mech_FEM_init real(pReal) :: HGCoeff = 0.0e-2_pReal - PetscInt, dimension(worldsize) :: localK + PetscInt, dimension(0:worldsize-1) :: localK real(pReal), dimension(3,3) :: & temp33_Real = 0.0_pReal real(pReal), dimension(4,8) :: & @@ -120,8 +120,8 @@ subroutine grid_mech_FEM_init ! initialize solver specific parts of PETSc call SNESCreate(PETSC_COMM_WORLD,mech_snes,ierr); CHKERRQ(ierr) call SNESSetOptionsPrefix(mech_snes,'mech_',ierr);CHKERRQ(ierr) - localK = 0 - localK(worldrank+1) = grid3 + localK = 0 + localK(worldrank) = grid3 call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,PETSC_COMM_WORLD,ierr) call DMDACreate3d(PETSC_COMM_WORLD, & DM_BOUNDARY_PERIODIC, DM_BOUNDARY_PERIODIC, DM_BOUNDARY_PERIODIC, & diff --git a/src/grid/grid_mech_spectral_polarisation.f90 b/src/grid/grid_mech_spectral_polarisation.f90 index 59ab84869..bdc65a8c5 100644 --- a/src/grid/grid_mech_spectral_polarisation.f90 +++ b/src/grid/grid_mech_spectral_polarisation.f90 @@ -100,7 +100,7 @@ subroutine grid_mech_spectral_polarisation_init FandF_tau, & ! overall pointer to solution data F, & ! specific (sub)pointer F_tau ! specific (sub)pointer - PetscInt, dimension(worldsize) :: localK + PetscInt, dimension(0:worldsize-1) :: localK integer(HID_T) :: fileHandle, groupHandle integer :: fileUnit character(len=pStringLen) :: fileName @@ -130,8 +130,8 @@ subroutine grid_mech_spectral_polarisation_init ! initialize solver specific parts of PETSc call SNESCreate(PETSC_COMM_WORLD,snes,ierr); CHKERRQ(ierr) call SNESSetOptionsPrefix(snes,'mech_',ierr);CHKERRQ(ierr) - localK = 0 - localK(worldrank+1) = grid3 + localK = 0 + localK(worldrank) = grid3 call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,PETSC_COMM_WORLD,ierr) call DMDACreate3d(PETSC_COMM_WORLD, & DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & ! cut off stencil at boundary diff --git a/src/grid/grid_thermal_spectral.f90 b/src/grid/grid_thermal_spectral.f90 index c381d837d..2d69f075e 100644 --- a/src/grid/grid_thermal_spectral.f90 +++ b/src/grid/grid_thermal_spectral.f90 @@ -55,7 +55,7 @@ contains !-------------------------------------------------------------------------------------------------- subroutine grid_thermal_spectral_init - PetscInt, dimension(worldsize) :: localK + PetscInt, dimension(0:worldsize-1) :: localK integer :: i, j, k, cell DM :: thermal_grid PetscScalar, dimension(:,:,:), pointer :: x_scal @@ -77,8 +77,8 @@ subroutine grid_thermal_spectral_init ! initialize solver specific parts of PETSc call SNESCreate(PETSC_COMM_WORLD,thermal_snes,ierr); CHKERRQ(ierr) call SNESSetOptionsPrefix(thermal_snes,'thermal_',ierr);CHKERRQ(ierr) - localK = 0 - localK(worldrank+1) = grid3 + localK = 0 + localK(worldrank) = grid3 call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,PETSC_COMM_WORLD,ierr) call DMDACreate3D(PETSC_COMM_WORLD, & DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & ! cut off stencil at boundary From 13150291966c96d8f02aeef7bbb267e5a2513688 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 12 Jan 2020 02:06:53 +0100 Subject: [PATCH 33/37] do not clutter comments --- python/damask/table.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/python/damask/table.py b/python/damask/table.py index 343ba638d..c40baa23c 100644 --- a/python/damask/table.py +++ b/python/damask/table.py @@ -22,7 +22,7 @@ class Table(): Additional, human-readable information. """ - self.comments = ['table.py v {}'.format(version)] if not comments else [c for c in comments] + self.comments = [c for c in comments] self.data = pd.DataFrame(data=data) self.shapes = shapes self.__label_condensed() @@ -79,8 +79,7 @@ class Table(): else: raise Exception - comments = ['table.py:from_ASCII v {}'.format(version)] - comments+= [f.readline()[:-1] for i in range(1,header)] + comments = [f.readline()[:-1] for i in range(1,header)] labels = f.readline().split() shapes = {} From 19e88df57131af9a82bce7d8da995e36bc21cb74 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 12 Jan 2020 07:53:41 +0100 Subject: [PATCH 34/37] polishing --- python/damask/table.py | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/python/damask/table.py b/python/damask/table.py index c40baa23c..a3c0a8af0 100644 --- a/python/damask/table.py +++ b/python/damask/table.py @@ -3,7 +3,6 @@ import re import pandas as pd import numpy as np -from . import version class Table(): """Store spreadsheet-like data.""" @@ -22,7 +21,7 @@ class Table(): Additional, human-readable information. """ - self.comments = [c for c in comments] + self.comments = [] if comments is None else [c for c in comments] self.data = pd.DataFrame(data=data) self.shapes = shapes self.__label_condensed() @@ -77,8 +76,8 @@ class Table(): if keyword == 'header': header = int(header) else: - raise Exception - + raise TypeError + comments = [f.readline()[:-1] for i in range(1,header)] labels = f.readline().split() @@ -138,7 +137,7 @@ class Table(): data = np.loadtxt(content) for c in range(data.shape[1]-10): - shapes['user_defined{}'.format(c+1)] = (1,) + shapes['n/a_{}'.format(c+1)] = (1,) return Table(data,shapes,comments) From 70b73762edc4b3c2e7a6831724c86e8fbe169a2b Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 12 Jan 2020 19:21:16 +0100 Subject: [PATCH 35/37] avoid warning due to change in default parameter --- python/damask/dadf5.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/damask/dadf5.py b/python/damask/dadf5.py index 73e785b10..5c9ebe08e 100644 --- a/python/damask/dadf5.py +++ b/python/damask/dadf5.py @@ -880,7 +880,7 @@ class DADF5(): else: nodes = vtk.vtkPoints() - with h5py.File(self.fname) as f: + with h5py.File(self.fname,'r') as f: nodes.SetData(numpy_support.numpy_to_vtk(f['/geometry/x_n'][()],deep=True)) vtk_geom = vtk.vtkUnstructuredGrid() From 26442ee7838b801caf2238a607c02e0315720868 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 12 Jan 2020 19:35:45 +0100 Subject: [PATCH 36/37] bugfixes missing import/outdated test --- PRIVATE | 2 +- python/damask/table.py | 1 + 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/PRIVATE b/PRIVATE index 3eeffba15..99b076706 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 3eeffba1513c0495032d0f33522e3aaa1ed77465 +Subproject commit 99b076706a186ec7deb051ae181c2d9697c799fc diff --git a/python/damask/table.py b/python/damask/table.py index a3c0a8af0..b3dfc2433 100644 --- a/python/damask/table.py +++ b/python/damask/table.py @@ -3,6 +3,7 @@ import re import pandas as pd import numpy as np +from . import version class Table(): """Store spreadsheet-like data.""" From 96b2b6f602053f0e1d2fc0e9669f291ed9ebcd28 Mon Sep 17 00:00:00 2001 From: Test User Date: Sun, 12 Jan 2020 20:41:01 +0100 Subject: [PATCH 37/37] [skip ci] updated version information after successful test of v2.0.3-1354-gef5a8a3a --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index cf29053c3..7712dccc1 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.3-1308-g612e8e34 +v2.0.3-1354-gef5a8a3a