Merge branch 'polishing' into 'development'

Polishing

See merge request damask/DAMASK!601
This commit is contained in:
Nikhil Prabhu 2022-06-16 09:31:02 +00:00
commit ac626db9df
23 changed files with 304 additions and 289 deletions

@ -1 +1 @@
Subproject commit 226d7b627e44247b800bce0d9eb7bef1aac6f537 Subproject commit b14f78e96a8e2986aaf6845b98ea77fec92bc997

View File

@ -0,0 +1,15 @@
type: dislotwin
references:
- N. Jia et al.,
Acta Materialia 60(3):1099-1115, 2012,
https://doi.org/10.1016/j.actamat.2011.10.047
- N. Jia et al.,
Acta Materialia 60:3415-3434, 2012,
https://doi.org/10.1016/j.actamat.2012.03.005
gamma_0_sb: 0.0001
tau_sb: 180.0e6 # tau_hat_sb
Q_sb: 4.0e-19 # Q_0
p_sb: 1.15
q_sb: 1.0

View File

@ -156,15 +156,15 @@ subroutine CLI_init
if (CLI_restartInc < 0 .or. stat /=0) then if (CLI_restartInc < 0 .or. stat /=0) then
print'(/,a)', ' ERROR: Could not parse restart increment: '//trim(arg) print'(/,a)', ' ERROR: Could not parse restart increment: '//trim(arg)
call quit(1) call quit(1)
endif end if
end select end select
if (err /= 0) call quit(1) if (err /= 0) call quit(1)
enddo end do
if (len_trim(loadcaseArg) == 0 .or. len_trim(geometryArg) == 0) then if (len_trim(loadcaseArg) == 0 .or. len_trim(geometryArg) == 0) then
print'(/,a)', ' ERROR: Please specify geometry AND load case (-h for help)' print'(/,a)', ' ERROR: Please specify geometry AND load case (-h for help)'
call quit(1) call quit(1)
endif end if
if (len_trim(workingDirArg) > 0) call setWorkingDirectory(trim(workingDirArg)) if (len_trim(workingDirArg) > 0) call setWorkingDirectory(trim(workingDirArg))
CLI_geomFile = getGeometryFile(geometryArg) CLI_geomFile = getGeometryFile(geometryArg)
@ -205,14 +205,14 @@ subroutine setWorkingDirectory(workingDirectoryArg)
else absolutePath else absolutePath
workingDirectory = getCWD() workingDirectory = getCWD()
workingDirectory = trim(workingDirectory)//'/'//workingDirectoryArg workingDirectory = trim(workingDirectory)//'/'//workingDirectoryArg
endif absolutePath end if absolutePath
workingDirectory = trim(rectifyPath(workingDirectory)) workingDirectory = trim(rectifyPath(workingDirectory))
error = setCWD(trim(workingDirectory)) error = setCWD(trim(workingDirectory))
if(error) then if(error) then
print*, 'ERROR: Invalid Working directory: '//trim(workingDirectory) print*, 'ERROR: Invalid Working directory: '//trim(workingDirectory)
call quit(1) call quit(1)
endif end if
end subroutine setWorkingDirectory end subroutine setWorkingDirectory
@ -256,7 +256,7 @@ function getGeometryFile(geometryParameter)
if (.not. file_exists) then if (.not. file_exists) then
print*, 'ERROR: Geometry file does not exists: '//trim(getGeometryFile) print*, 'ERROR: Geometry file does not exists: '//trim(getGeometryFile)
call quit(1) call quit(1)
endif end if
end function getGeometryFile end function getGeometryFile
@ -279,7 +279,7 @@ function getLoadCaseFile(loadCaseParameter)
if (.not. file_exists) then if (.not. file_exists) then
print*, 'ERROR: Load case file does not exists: '//trim(getLoadCaseFile) print*, 'ERROR: Load case file does not exists: '//trim(getLoadCaseFile)
call quit(1) call quit(1)
endif end if
end function getLoadCaseFile end function getLoadCaseFile
@ -300,14 +300,14 @@ function rectifyPath(path)
l = len_trim(rectifyPath) l = len_trim(rectifyPath)
do i = l,3,-1 do i = l,3,-1
if (rectifyPath(i-2:i) == '/./') rectifyPath(i-1:l) = rectifyPath(i+1:l)//' ' if (rectifyPath(i-2:i) == '/./') rectifyPath(i-1:l) = rectifyPath(i+1:l)//' '
enddo end do
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! remove // from path ! remove // from path
l = len_trim(rectifyPath) l = len_trim(rectifyPath)
do i = l,2,-1 do i = l,2,-1
if (rectifyPath(i-1:i) == '//') rectifyPath(i-1:l) = rectifyPath(i:l)//' ' if (rectifyPath(i-1:i) == '//') rectifyPath(i-1:l) = rectifyPath(i:l)//' '
enddo end do
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! remove ../ and corresponding directory from rectifyPath ! remove ../ and corresponding directory from rectifyPath
@ -321,9 +321,9 @@ function rectifyPath(path)
k = len_trim(rectifyPath) k = len_trim(rectifyPath)
rectifyPath(j+1:k-1) = rectifyPath(j+2:k) rectifyPath(j+1:k-1) = rectifyPath(j+2:k)
rectifyPath(k:k) = ' ' rectifyPath(k:k) = ' '
endif end if
i = j+index(rectifyPath(j+1:l),'../') i = j+index(rectifyPath(j+1:l),'../')
enddo end do
if(len_trim(rectifyPath) == 0) rectifyPath = '/' if(len_trim(rectifyPath) == 0) rectifyPath = '/'
rectifyPath = trim(rectifyPath) rectifyPath = trim(rectifyPath)
@ -349,10 +349,10 @@ function makeRelativePath(a,b)
do i = 1, min(len_trim(a_cleaned),len_trim(rectifyPath(b_cleaned))) do i = 1, min(len_trim(a_cleaned),len_trim(rectifyPath(b_cleaned)))
if (a_cleaned(i:i) /= b_cleaned(i:i)) exit if (a_cleaned(i:i) /= b_cleaned(i:i)) exit
if (a_cleaned(i:i) == '/') posLastCommonSlash = i if (a_cleaned(i:i) == '/') posLastCommonSlash = i
enddo end do
do i = posLastCommonSlash+1,len_trim(a_cleaned) do i = posLastCommonSlash+1,len_trim(a_cleaned)
if (a_cleaned(i:i) == '/') remainingSlashes = remainingSlashes + 1 if (a_cleaned(i:i) == '/') remainingSlashes = remainingSlashes + 1
enddo end do
makeRelativePath = repeat('..'//'/',remainingSlashes)//b_cleaned(posLastCommonSlash+1:len_trim(b_cleaned)) makeRelativePath = repeat('..'//'/',remainingSlashes)//b_cleaned(posLastCommonSlash+1:len_trim(b_cleaned))

View File

@ -509,7 +509,7 @@ subroutine HDF5_addAttribute_str_array(loc_id,attrLabel,attrValue,path)
do i=1,size(attrValue) do i=1,size(attrValue)
attrValue_(i) = attrValue(i)//C_NULL_CHAR attrValue_(i) = attrValue(i)//C_NULL_CHAR
ptr(i) = c_loc(attrValue_(i)) ptr(i) = c_loc(attrValue_(i))
enddo end do
call H5Screate_simple_f(1,shape(attrValue_,kind=HSIZE_T),space_id,hdferr,shape(attrValue_,kind=HSIZE_T)) call H5Screate_simple_f(1,shape(attrValue_,kind=HSIZE_T),space_id,hdferr,shape(attrValue_,kind=HSIZE_T))
if(hdferr < 0) error stop 'HDF5 error' if(hdferr < 0) error stop 'HDF5 error'

View File

@ -92,15 +92,15 @@ function IO_readlines(fileName) result(fileContent)
if (.not. warned) then if (.not. warned) then
call IO_warning(207,trim(fileName),label1='line',ID1=l) call IO_warning(207,trim(fileName),label1='line',ID1=l)
warned = .true. warned = .true.
endif end if
else else
line = rawData(startPos:endpos) line = rawData(startPos:endpos)
endif end if
startPos = endPos + 2 ! jump to next line start startPos = endPos + 2 ! jump to next line start
fileContent(l) = trim(line)//'' fileContent(l) = trim(line)//''
l = l + 1 l = l + 1
enddo end do
end function IO_readlines end function IO_readlines
@ -129,7 +129,7 @@ function IO_read(fileName) result(fileContent)
if (fileLength==0) then if (fileLength==0) then
close(fileUnit) close(fileUnit)
return return
endif end if
read(fileUnit,iostat=myStat) fileContent read(fileUnit,iostat=myStat) fileContent
if (myStat /= 0) call IO_error(102,trim(fileName)) if (myStat /= 0) call IO_error(102,trim(fileName))
@ -183,8 +183,8 @@ pure function IO_stringPos(string)
endOfString: if (right < left) then endOfString: if (right < left) then
IO_stringPos(IO_stringPos(1)*2+1) = len_trim(string) IO_stringPos(IO_stringPos(1)*2+1) = len_trim(string)
exit exit
endif endOfString end if endOfString
enddo end do
end function IO_stringPos end function IO_stringPos
@ -204,7 +204,7 @@ function IO_stringValue(string,chunkPos,myChunk)
call IO_error(110,'IO_stringValue: "'//trim(string)//'"',label1='chunk',ID1=myChunk) call IO_error(110,'IO_stringValue: "'//trim(string)//'"',label1='chunk',ID1=myChunk)
else validChunk else validChunk
IO_stringValue = string(chunkPos(myChunk*2):chunkPos(myChunk*2+1)) IO_stringValue = string(chunkPos(myChunk*2):chunkPos(myChunk*2+1))
endif validChunk end if validChunk
end function IO_stringValue end function IO_stringValue
@ -257,8 +257,8 @@ pure function IO_lc(string)
IO_lc(i:i) = LOWER(n:n) IO_lc(i:i) = LOWER(n:n)
else else
IO_lc(i:i) = string(i:i) IO_lc(i:i) = string(i:i)
endif end if
enddo end do
end function IO_lc end function IO_lc
@ -280,7 +280,7 @@ function IO_rmComment(line)
IO_rmComment = trim(line) IO_rmComment = trim(line)
else else
IO_rmComment = trim(line(:split-1)) IO_rmComment = trim(line(:split-1))
endif end if
end function IO_rmComment end function IO_rmComment
@ -302,7 +302,7 @@ integer function IO_stringAsInt(string)
else valid else valid
IO_stringAsInt = 0 IO_stringAsInt = 0
call IO_error(111,string) call IO_error(111,string)
endif valid end if valid
end function IO_stringAsInt end function IO_stringAsInt
@ -324,7 +324,7 @@ real(pReal) function IO_stringAsFloat(string)
else valid else valid
IO_stringAsFloat = 0.0_pReal IO_stringAsFloat = 0.0_pReal
call IO_error(112,string) call IO_error(112,string)
endif valid end if valid
end function IO_stringAsFloat end function IO_stringAsFloat
@ -344,7 +344,7 @@ logical function IO_stringAsBool(string)
else else
IO_stringAsBool = .false. IO_stringAsBool = .false.
call IO_error(113,string) call IO_error(113,string)
endif end if
end function IO_stringAsBool end function IO_stringAsBool
@ -598,7 +598,7 @@ pure function CRLF2LF(string)
CRLF2LF(c-n:c-n) = string(c:c) CRLF2LF(c-n:c-n) = string(c:c)
if (c == len_trim(string)) exit if (c == len_trim(string)) exit
if (string(c:c+1) == CR//LF) n = n + 1 if (string(c:c+1) == CR//LF) n = n + 1
enddo end do
CRLF2LF = CRLF2LF(:c-n) CRLF2LF = CRLF2LF(:c-n)
@ -639,17 +639,17 @@ subroutine panel(paneltype,ID,msg,ext_msg,label1,ID1,label2,ID2)
write(formatString,'(a,i3.3,a,i3.3,a)') '(1x,a4,a',max(1,len_trim(ext_msg)),',',& write(formatString,'(a,i3.3,a,i3.3,a)') '(1x,a4,a',max(1,len_trim(ext_msg)),',',&
max(1,panelwidth+3-len_trim(ext_msg)-4),'x,a)' max(1,panelwidth+3-len_trim(ext_msg)-4),'x,a)'
write(IO_STDERR,formatString) '│ ',trim(ext_msg), '│' write(IO_STDERR,formatString) '│ ',trim(ext_msg), '│'
endif end if
if (present(label1)) then if (present(label1)) then
write(formatString,'(a,i3.3,a,i3.3,a)') '(1x,a7,a',max(1,len_trim(label1)),',i9,',& write(formatString,'(a,i3.3,a,i3.3,a)') '(1x,a7,a',max(1,len_trim(label1)),',i9,',&
max(1,panelwidth+3-len_trim(label1)-9-7),'x,a)' max(1,panelwidth+3-len_trim(label1)-9-7),'x,a)'
write(IO_STDERR,formatString) '│ at ',trim(label1),ID1, '│' write(IO_STDERR,formatString) '│ at ',trim(label1),ID1, '│'
endif end if
if (present(label2)) then if (present(label2)) then
write(formatString,'(a,i3.3,a,i3.3,a)') '(1x,a7,a',max(1,len_trim(label2)),',i9,',& write(formatString,'(a,i3.3,a,i3.3,a)') '(1x,a7,a',max(1,len_trim(label2)),',i9,',&
max(1,panelwidth+3-len_trim(label2)-9-7),'x,a)' max(1,panelwidth+3-len_trim(label2)-9-7),'x,a)'
write(IO_STDERR,formatString) '│ at ',trim(label2),ID2, '│' write(IO_STDERR,formatString) '│ at ',trim(label2),ID2, '│'
endif end if
write(formatString,'(a,i2.2,a)') '(a,',max(1,panelwidth),'x,a)' write(formatString,'(a,i2.2,a)') '(a,',max(1,panelwidth),'x,a)'
write(IO_STDERR,formatString) ' │', '│' write(IO_STDERR,formatString) ' │', '│'
write(IO_STDERR,'(a)') ' └'//DIVIDER//'┘' write(IO_STDERR,'(a)') ' └'//DIVIDER//'┘'

View File

@ -103,7 +103,7 @@ recursive function parse_flow(YAML_flow) result(node)
class is (tDict) class is (tDict)
call node%set(key,myVal) call node%set(key,myVal)
end select end select
enddo end do
elseif (flow_string(1:1) == '[') then ! start of a list elseif (flow_string(1:1) == '[') then ! start of a list
e = 1 e = 1
allocate(tList::node) allocate(tList::node)
@ -116,7 +116,7 @@ recursive function parse_flow(YAML_flow) result(node)
class is (tList) class is (tList)
call node%append(myVal) call node%append(myVal)
end select end select
enddo end do
else ! scalar value else ! scalar value
allocate(tScalar::node) allocate(tScalar::node)
select type (node) select type (node)
@ -156,7 +156,7 @@ integer function find_end(str,e_char)
N_sq = N_sq - merge(1,0,str(i:i) == ']') N_sq = N_sq - merge(1,0,str(i:i) == ']')
N_cu = N_cu - merge(1,0,str(i:i) == '}') N_cu = N_cu - merge(1,0,str(i:i) == '}')
i = i + 1 i = i + 1
enddo end do
find_end = i find_end = i
end function find_end end function find_end
@ -333,7 +333,7 @@ subroutine skip_empty_lines(blck,s_blck)
do while(empty .and. len_trim(blck(s_blck:)) /= 0) do while(empty .and. len_trim(blck(s_blck:)) /= 0)
empty = len_trim(IO_rmComment(blck(s_blck:s_blck + index(blck(s_blck:),IO_EOL) - 2))) == 0 empty = len_trim(IO_rmComment(blck(s_blck:s_blck + index(blck(s_blck:),IO_EOL) - 2))) == 0
if(empty) s_blck = s_blck + index(blck(s_blck:),IO_EOL) if(empty) s_blck = s_blck + index(blck(s_blck:),IO_EOL)
enddo end do
end subroutine skip_empty_lines end subroutine skip_empty_lines
@ -387,7 +387,7 @@ logical function flow_is_closed(str,e_char)
N_cu = N_cu + merge(1,0,line(i:i) == '{') N_cu = N_cu + merge(1,0,line(i:i) == '{')
N_sq = N_sq - merge(1,0,line(i:i) == ']') N_sq = N_sq - merge(1,0,line(i:i) == ']')
N_cu = N_cu - merge(1,0,line(i:i) == '}') N_cu = N_cu - merge(1,0,line(i:i) == '}')
enddo end do
end function flow_is_closed end function flow_is_closed
@ -410,7 +410,7 @@ subroutine remove_line_break(blck,s_blck,e_char,flow_line)
flow_line = flow_line//IO_rmComment(blck(s_blck:s_blck + index(blck(s_blck:),IO_EOL) - 2))//' ' flow_line = flow_line//IO_rmComment(blck(s_blck:s_blck + index(blck(s_blck:),IO_EOL) - 2))//' '
line_end = flow_is_closed(flow_line,e_char) line_end = flow_is_closed(flow_line,e_char)
s_blck = s_blck + index(blck(s_blck:),IO_EOL) s_blck = s_blck + index(blck(s_blck:),IO_EOL)
enddo end do
end subroutine remove_line_break end subroutine remove_line_break
@ -439,7 +439,7 @@ subroutine list_item_inline(blck,s_blck,inline,offset)
inline = inline//' '//trim(adjustl(IO_rmComment(blck(s_blck:s_blck + index(blck(s_blck:),IO_EOL) - 2)))) inline = inline//' '//trim(adjustl(IO_rmComment(blck(s_blck:s_blck + index(blck(s_blck:),IO_EOL) - 2))))
s_blck = s_blck + index(blck(s_blck:),IO_EOL) s_blck = s_blck + index(blck(s_blck:),IO_EOL)
indent_next = indentDepth(blck(s_blck:)) indent_next = indentDepth(blck(s_blck:))
enddo end do
if(scan(inline,",") > 0) inline = '"'//inline//'"' if(scan(inline,",") > 0) inline = '"'//inline//'"'
@ -481,7 +481,7 @@ recursive subroutine line_isFlow(flow,s_flow,line)
flow(s_flow:s_flow+1) = ', ' flow(s_flow:s_flow+1) = ', '
s_flow = s_flow +2 s_flow = s_flow +2
s = s + find_end(line(s+1:),']') s = s + find_end(line(s+1:),']')
enddo end do
s_flow = s_flow - 1 s_flow = s_flow - 1
if (flow(s_flow-1:s_flow-1) == ',') s_flow = s_flow - 1 if (flow(s_flow-1:s_flow-1) == ',') s_flow = s_flow - 1
flow(s_flow:s_flow) = ']' flow(s_flow:s_flow) = ']'
@ -498,7 +498,7 @@ recursive subroutine line_isFlow(flow,s_flow,line)
flow(s_flow:s_flow+1) = ', ' flow(s_flow:s_flow+1) = ', '
s_flow = s_flow +2 s_flow = s_flow +2
s = s + find_end(line(s+1:),'}') s = s + find_end(line(s+1:),'}')
enddo end do
s_flow = s_flow -1 s_flow = s_flow -1
if(flow(s_flow-1:s_flow-1) == ',') s_flow = s_flow -1 if(flow(s_flow-1:s_flow-1) == ',') s_flow = s_flow -1
flow(s_flow:s_flow) = '}' flow(s_flow:s_flow) = '}'
@ -646,7 +646,7 @@ recursive subroutine lst(blck,flow,s_blck,s_flow,offset)
s_flow = s_flow + 2 s_flow = s_flow + 2
end if end if
enddo end do
s_flow = s_flow - 1 s_flow = s_flow - 1
if (flow(s_flow-1:s_flow-1) == ',') s_flow = s_flow - 1 if (flow(s_flow-1:s_flow-1) == ',') s_flow = s_flow - 1
@ -733,7 +733,7 @@ recursive subroutine dct(blck,flow,s_blck,s_flow,offset)
flow(s_flow:s_flow) = ' ' flow(s_flow:s_flow) = ' '
s_flow = s_flow + 1 s_flow = s_flow + 1
offset = 0 offset = 0
enddo end do
s_flow = s_flow - 1 s_flow = s_flow - 1
if (flow(s_flow-1:s_flow-1) == ',') s_flow = s_flow - 1 if (flow(s_flow-1:s_flow-1) == ',') s_flow = s_flow - 1

View File

@ -411,7 +411,7 @@ function tNode_get_byIndex(self,i) result(node)
do j = 2,i do j = 2,i
item => item%next item => item%next
enddo end do
node => item%node node => item%node
end function tNode_get_byIndex end function tNode_get_byIndex
@ -681,7 +681,7 @@ function tNode_contains(self,k) result(exists)
exists = .true. exists = .true.
return return
end if end if
enddo end do
class is(tList) class is(tList)
list => self%asList() list => self%asList()
do j=1, list%length do j=1, list%length
@ -689,7 +689,7 @@ function tNode_contains(self,k) result(exists)
exists = .true. exists = .true.
return return
end if end if
enddo end do
class default class default
call IO_error(706,ext_msg='Expected list or dict') call IO_error(706,ext_msg='Expected list or dict')
end select end select
@ -731,7 +731,7 @@ function tNode_get_byKey(self,k,defaultVal) result(node)
end if end if
item => item%next item => item%next
j = j + 1 j = j + 1
enddo end do
if (.not. found) then if (.not. found) then
call IO_error(143,ext_msg=k) call IO_error(143,ext_msg=k)
@ -1333,7 +1333,7 @@ function tList_as1dString(self)
scalar => item%node%asScalar() scalar => item%node%asScalar()
tList_as1dString(i) = scalar%asString() tList_as1dString(i) = scalar%asString()
item => item%next item => item%next
enddo end do
end function tList_as1dString end function tList_as1dString
@ -1384,7 +1384,7 @@ subroutine tDict_set(self,key,node)
searchExisting: do while (associated(item%next)) searchExisting: do while (associated(item%next))
if (item%key == key) exit if (item%key == key) exit
item => item%next item => item%next
enddo searchExisting end do searchExisting
if (item%key /= key) then if (item%key /= key) then
allocate(item%next) allocate(item%next)
item => item%next item => item%next

View File

@ -68,7 +68,7 @@ subroutine discretization_init(materialAt,&
discretization_sharedNodesBegin = sharedNodesBegin discretization_sharedNodesBegin = sharedNodesBegin
else else
discretization_sharedNodesBegin = size(discretization_NodeCoords0,2) discretization_sharedNodesBegin = size(discretization_NodeCoords0,2)
endif end if
end subroutine discretization_init end subroutine discretization_init

View File

@ -365,7 +365,7 @@ subroutine homogenization_results
call thermal_results(ho,group) call thermal_results(ho,group)
end if end if
enddo end do
end subroutine homogenization_results end subroutine homogenization_results
@ -383,7 +383,7 @@ subroutine homogenization_forward
homogState (ho)%state0 = homogState (ho)%state homogState (ho)%state0 = homogState (ho)%state
if(damageState_h(ho)%sizeState > 0) & if(damageState_h(ho)%sizeState > 0) &
damageState_h(ho)%state0 = damageState_h(ho)%state damageState_h(ho)%state0 = damageState_h(ho)%state
enddo end do
end subroutine homogenization_forward end subroutine homogenization_forward
@ -408,7 +408,7 @@ subroutine homogenization_restartWrite(fileHandle)
call HDF5_closeGroup(groupHandle(2)) call HDF5_closeGroup(groupHandle(2))
enddo end do
call HDF5_closeGroup(groupHandle(1)) call HDF5_closeGroup(groupHandle(1))
@ -435,7 +435,7 @@ subroutine homogenization_restartRead(fileHandle)
call HDF5_closeGroup(groupHandle(2)) call HDF5_closeGroup(groupHandle(2))
enddo end do
call HDF5_closeGroup(groupHandle(1)) call HDF5_closeGroup(groupHandle(1))
@ -476,7 +476,7 @@ subroutine parseHomogenization
case default case default
call IO_error(500,ext_msg=homogThermal%get_asString('type')) call IO_error(500,ext_msg=homogThermal%get_asString('type'))
end select end select
endif end if
if (homog%contains('damage')) then if (homog%contains('damage')) then
homogDamage => homog%get('damage') homogDamage => homog%get('damage')
@ -486,8 +486,8 @@ subroutine parseHomogenization
case default case default
call IO_error(500,ext_msg=homogDamage%get_asString('type')) call IO_error(500,ext_msg=homogDamage%get_asString('type'))
end select end select
endif end if
enddo end do
end subroutine parseHomogenization end subroutine parseHomogenization

View File

@ -65,9 +65,9 @@ module subroutine damage_init()
allocate(damageState_h(ho)%state (1,Nmembers), source=1.0_pReal) allocate(damageState_h(ho)%state (1,Nmembers), source=1.0_pReal)
else else
prm%output = emptyStringArray prm%output = emptyStringArray
endif end if
end associate end associate
enddo end do
call pass_init() call pass_init()
@ -174,7 +174,7 @@ module subroutine damage_results(ho,group)
call results_writeDataset(damagestate_h(ho)%state(1,:),group,prm%output(o),& call results_writeDataset(damagestate_h(ho)%state(1,:),group,prm%output(o),&
'damage indicator','-') 'damage indicator','-')
end select end select
enddo outputsLoop end do outputsLoop
end associate end associate
end subroutine damage_results end subroutine damage_results

View File

@ -561,7 +561,7 @@ module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy)
*cosh(prm%c_alpha*nDefNorm) & *cosh(prm%c_alpha*nDefNorm) &
*0.5_pReal*nVect(l)*nDef(i,k)/nDefNorm*math_LeviCivita(k,l,j) & *0.5_pReal*nVect(l)*nDef(i,k)/nDefNorm*math_LeviCivita(k,l,j) &
*tanh(nDefNorm/num%xSmoo) *tanh(nDefNorm/num%xSmoo)
end do; end do;enddo; end do end do; end do;end do; end do
end do interfaceLoop end do interfaceLoop

View File

@ -484,8 +484,8 @@ function lattice_characteristicShear_Twin(Ntwin,lattice,CoverA) result(character
case default case default
call IO_error(137,ext_msg='lattice_characteristicShear_Twin: '//trim(lattice)) call IO_error(137,ext_msg='lattice_characteristicShear_Twin: '//trim(lattice))
end select end select
enddo mySystems end do mySystems
enddo myFamilies end do myFamilies
end function lattice_characteristicShear_Twin end function lattice_characteristicShear_Twin
@ -523,7 +523,7 @@ function lattice_C66_twin(Ntwin,C66,lattice,CoverA)
do i = 1, sum(Ntwin) do i = 1, sum(Ntwin)
call R%fromAxisAngle([coordinateSystem(1:3,2,i),PI],P=1) ! ToDo: Why always 180 deg? call R%fromAxisAngle([coordinateSystem(1:3,2,i),PI],P=1) ! ToDo: Why always 180 deg?
lattice_C66_twin(1:6,1:6,i) = R%rotStiffness(C66) lattice_C66_twin(1:6,1:6,i) = R%rotStiffness(C66)
enddo end do
end function lattice_C66_twin end function lattice_C66_twin
@ -572,19 +572,19 @@ function lattice_C66_trans(Ntrans,C_parent66,lattice_target, &
C_target_unrotated66 = C_parent66 C_target_unrotated66 = C_parent66
else else
call IO_error(137,ext_msg='lattice_C66_trans : '//trim(lattice_target)) call IO_error(137,ext_msg='lattice_C66_trans : '//trim(lattice_target))
endif end if
do i = 1,6 do i = 1,6
if (abs(C_target_unrotated66(i,i))<tol_math_check) & if (abs(C_target_unrotated66(i,i))<tol_math_check) &
call IO_error(135,'matrix diagonal in transformation',label1='entry',ID1=i) call IO_error(135,'matrix diagonal in transformation',label1='entry',ID1=i)
enddo end do
call buildTransformationSystem(Q,S,Ntrans,cOverA_trans,a_cF,a_cI) call buildTransformationSystem(Q,S,Ntrans,cOverA_trans,a_cF,a_cI)
do i = 1,sum(Ntrans) do i = 1,sum(Ntrans)
call R%fromMatrix(Q(1:3,1:3,i)) call R%fromMatrix(Q(1:3,1:3,i))
lattice_C66_trans(1:6,1:6,i) = R%rotStiffness(C_target_unrotated66) lattice_C66_trans(1:6,1:6,i) = R%rotStiffness(C_target_unrotated66)
enddo end do
end function lattice_C66_trans end function lattice_C66_trans
@ -632,7 +632,7 @@ function lattice_nonSchmidMatrix(Nslip,nonSchmidCoefficients,sense) result(nonSc
math_cross(normal, direction)) math_cross(normal, direction))
if (size(nonSchmidCoefficients)>5) nonSchmidMatrix(1:3,1:3,i) = nonSchmidMatrix(1:3,1:3,i) & if (size(nonSchmidCoefficients)>5) nonSchmidMatrix(1:3,1:3,i) = nonSchmidMatrix(1:3,1:3,i) &
+ nonSchmidCoefficients(6) * math_outer(direction, direction) + nonSchmidCoefficients(6) * math_outer(direction, direction)
enddo end do
end function lattice_nonSchmidMatrix end function lattice_nonSchmidMatrix
@ -1432,7 +1432,7 @@ function lattice_SchmidMatrix_slip(Nslip,lattice,cOverA) result(SchmidMatrix)
SchmidMatrix(1:3,1:3,i) = math_outer(coordinateSystem(1:3,1,i),coordinateSystem(1:3,2,i)) SchmidMatrix(1:3,1:3,i) = math_outer(coordinateSystem(1:3,1,i),coordinateSystem(1:3,2,i))
if (abs(math_trace33(SchmidMatrix(1:3,1:3,i))) > tol_math_check) & if (abs(math_trace33(SchmidMatrix(1:3,1:3,i))) > tol_math_check) &
error stop 'dilatational Schmid matrix for slip' error stop 'dilatational Schmid matrix for slip'
enddo end do
end function lattice_SchmidMatrix_slip end function lattice_SchmidMatrix_slip
@ -1479,7 +1479,7 @@ function lattice_SchmidMatrix_twin(Ntwin,lattice,cOverA) result(SchmidMatrix)
SchmidMatrix(1:3,1:3,i) = math_outer(coordinateSystem(1:3,1,i),coordinateSystem(1:3,2,i)) SchmidMatrix(1:3,1:3,i) = math_outer(coordinateSystem(1:3,1,i),coordinateSystem(1:3,2,i))
if (abs(math_trace33(SchmidMatrix(1:3,1:3,i))) > tol_math_check) & if (abs(math_trace33(SchmidMatrix(1:3,1:3,i))) > tol_math_check) &
error stop 'dilatational Schmid matrix for twin' error stop 'dilatational Schmid matrix for twin'
enddo end do
end function lattice_SchmidMatrix_twin end function lattice_SchmidMatrix_twin
@ -1552,7 +1552,7 @@ function lattice_SchmidMatrix_cleavage(Ncleavage,lattice,cOverA) result(SchmidMa
SchmidMatrix(1:3,1:3,1,i) = math_outer(coordinateSystem(1:3,1,i),coordinateSystem(1:3,2,i)) SchmidMatrix(1:3,1:3,1,i) = math_outer(coordinateSystem(1:3,1,i),coordinateSystem(1:3,2,i))
SchmidMatrix(1:3,1:3,2,i) = math_outer(coordinateSystem(1:3,3,i),coordinateSystem(1:3,2,i)) SchmidMatrix(1:3,1:3,2,i) = math_outer(coordinateSystem(1:3,3,i),coordinateSystem(1:3,2,i))
SchmidMatrix(1:3,1:3,3,i) = math_outer(coordinateSystem(1:3,2,i),coordinateSystem(1:3,2,i)) SchmidMatrix(1:3,1:3,3,i) = math_outer(coordinateSystem(1:3,2,i),coordinateSystem(1:3,2,i))
enddo end do
end function lattice_SchmidMatrix_cleavage end function lattice_SchmidMatrix_cleavage
@ -1719,8 +1719,8 @@ pure function lattice_symmetrize_C66(C66,lattice) result(C66_sym)
do i = 1, 6 do i = 1, 6
do j = i+1, 6 do j = i+1, 6
C66_sym(j,i) = C66_sym(i,j) C66_sym(j,i) = C66_sym(i,j)
enddo end do
enddo end do
end function lattice_symmetrize_C66 end function lattice_symmetrize_C66
@ -1782,7 +1782,7 @@ function slipProjection_transverse(Nslip,lattice,cOverA) result(projection)
do i=1, sum(Nslip); do j=1, sum(Nslip) do i=1, sum(Nslip); do j=1, sum(Nslip)
projection(i,j) = abs(math_inner(n(:,i),t(:,j))) projection(i,j) = abs(math_inner(n(:,i),t(:,j)))
enddo; enddo end do; end do
end function slipProjection_transverse end function slipProjection_transverse
@ -1806,7 +1806,7 @@ function slipProjection_direction(Nslip,lattice,cOverA) result(projection)
do i=1, sum(Nslip); do j=1, sum(Nslip) do i=1, sum(Nslip); do j=1, sum(Nslip)
projection(i,j) = abs(math_inner(n(:,i),d(:,j))) projection(i,j) = abs(math_inner(n(:,i),d(:,j)))
enddo; enddo end do; end do
end function slipProjection_direction end function slipProjection_direction
@ -1890,8 +1890,8 @@ function buildInteraction(reacting_used,acting_used,reacting_max,acting_max,valu
buildInteraction(l,k) = values(matrix(i,j)) buildInteraction(l,k) = values(matrix(i,j))
enddo; enddo end do; end do
enddo; enddo end do; end do
end function buildInteraction end function buildInteraction
@ -1957,8 +1957,8 @@ function buildCoordinateSystem(active,potential,system,lattice,cOverA)
buildCoordinateSystem(1:3,3,a) = math_cross(direction/norm2(direction),& buildCoordinateSystem(1:3,3,a) = math_cross(direction/norm2(direction),&
normal /norm2(normal)) normal /norm2(normal))
enddo activeSystems end do activeSystems
enddo activeFamilies end do activeFamilies
end function buildCoordinateSystem end function buildCoordinateSystem
@ -2066,7 +2066,7 @@ subroutine buildTransformationSystem(Q,S,Ntrans,cOverA,a_cF,a_cI)
U = (a_cI/a_cF) * (math_outer(x,x) + (math_outer(y,y)+math_outer(z,z)) * sqrt(2.0_pReal)) U = (a_cI/a_cF) * (math_outer(x,x) + (math_outer(y,y)+math_outer(z,z)) * sqrt(2.0_pReal))
Q(1:3,1:3,i) = matmul(R%asMatrix(),B%asMatrix()) Q(1:3,1:3,i) = matmul(R%asMatrix(),B%asMatrix())
S(1:3,1:3,i) = matmul(R%asMatrix(),U) - MATH_I3 S(1:3,1:3,i) = matmul(R%asMatrix(),U) - MATH_I3
enddo end do
else if (present(cOverA)) then else if (present(cOverA)) then
ss = MATH_I3 ss = MATH_I3
sd = MATH_I3 sd = MATH_I3
@ -2125,7 +2125,7 @@ function getlabels(active,potential,system) result(labels)
write(label(i+1:i+2),'(I2.1)') int(system(j,p)) write(label(i+1:i+2),'(I2.1)') int(system(j,p))
label(i+3:i+3) = ' ' label(i+3:i+3) = ' '
i = i + 3 i = i + 3
enddo direction end do direction
label(i:i) = ']' label(i:i) = ']'
i = i +1 i = i +1
@ -2134,13 +2134,13 @@ function getlabels(active,potential,system) result(labels)
write(label(i+1:i+2),'(I2.1)') int(system(j,p)) write(label(i+1:i+2),'(I2.1)') int(system(j,p))
label(i+3:i+3) = ' ' label(i+3:i+3) = ' '
i = i + 3 i = i + 3
enddo normal end do normal
label(i:i) = ')' label(i:i) = ')'
labels(a) = label labels(a) = label
enddo activeSystems end do activeSystems
enddo activeFamilies end do activeFamilies
end function getlabels end function getlabels
@ -2170,7 +2170,7 @@ pure function lattice_equivalent_nu(C,assumption) result(nu)
/ (S(1,1)+S(2,2)+S(3,3) +2.0_pReal*(S(1,2)+S(2,3)+S(1,3))) / (S(1,1)+S(2,2)+S(3,3) +2.0_pReal*(S(1,2)+S(2,3)+S(1,3)))
else else
error stop 'invalid assumption' error stop 'invalid assumption'
endif end if
mu = lattice_equivalent_mu(C,assumption) mu = lattice_equivalent_mu(C,assumption)
nu = (1.5_pReal*K-mu)/(3.0_pReal*K+mu) nu = (1.5_pReal*K-mu)/(3.0_pReal*K+mu)
@ -2202,7 +2202,7 @@ pure function lattice_equivalent_mu(C,assumption) result(mu)
/ (4.0_pReal*(S(1,1)+S(2,2)+S(3,3)) -4.0_pReal*(S(1,2)+S(2,3)+S(1,3)) +3.0_pReal*(S(4,4)+S(5,5)+S(6,6))) / (4.0_pReal*(S(1,1)+S(2,2)+S(3,3)) -4.0_pReal*(S(1,2)+S(2,3)+S(1,3)) +3.0_pReal*(S(4,4)+S(5,5)+S(6,6)))
else else
error stop 'invalid assumption' error stop 'invalid assumption'
endif end if
end function lattice_equivalent_mu end function lattice_equivalent_mu
@ -2266,7 +2266,7 @@ subroutine selfTest
if (any(dNeq(T(1,1),[T_hP(1,1),T_hP(2,2)]))) error stop 'Symmetry33_11-22/hP' if (any(dNeq(T(1,1),[T_hP(1,1),T_hP(2,2)]))) error stop 'Symmetry33_11-22/hP'
if (any(dNeq(T(1,1),[T_tI(1,1),T_tI(2,2)]))) error stop 'Symmetry33_11-22/tI' if (any(dNeq(T(1,1),[T_tI(1,1),T_tI(2,2)]))) error stop 'Symmetry33_11-22/tI'
enddo end do
call random_number(C) call random_number(C)
C(1,1) = C(1,1) + C(1,2) + 0.1_pReal C(1,1) = C(1,1) + C(1,2) + 0.1_pReal

View File

@ -135,25 +135,25 @@ pure recursive subroutine math_sort(a, istart, iend, sortDim)
s = istart s = istart
else else
s = lbound(a,2) s = lbound(a,2)
endif end if
if (present(iend)) then if (present(iend)) then
e = iend e = iend
else else
e = ubound(a,2) e = ubound(a,2)
endif end if
if (present(sortDim)) then if (present(sortDim)) then
d = sortDim d = sortDim
else else
d = 1 d = 1
endif end if
if (s < e) then if (s < e) then
call qsort_partition(a,ipivot, s,e, d) call qsort_partition(a,ipivot, s,e, d)
call math_sort(a, s, ipivot-1, d) call math_sort(a, s, ipivot-1, d)
call math_sort(a, ipivot+1, e, d) call math_sort(a, ipivot+1, e, d)
endif end if
contains contains
@ -175,11 +175,11 @@ pure recursive subroutine math_sort(a, istart, iend, sortDim)
! find the first element on the right side less than or equal to the pivot point ! find the first element on the right side less than or equal to the pivot point
do j = iend, istart, -1 do j = iend, istart, -1
if (a(sort,j) <= a(sort,istart)) exit if (a(sort,j) <= a(sort,istart)) exit
enddo end do
! find the first element on the left side greater than the pivot point ! find the first element on the left side greater than the pivot point
do i = istart, iend do i = istart, iend
if (a(sort,i) > a(sort,istart)) exit if (a(sort,i) > a(sort,istart)) exit
enddo end do
cross: if (i >= j) then ! exchange left value with pivot and return with the partition index cross: if (i >= j) then ! exchange left value with pivot and return with the partition index
tmp = a(:,istart) tmp = a(:,istart)
a(:,istart) = a(:,j) a(:,istart) = a(:,j)
@ -190,8 +190,8 @@ pure recursive subroutine math_sort(a, istart, iend, sortDim)
tmp = a(:,i) tmp = a(:,i)
a(:,i) = a(:,j) a(:,i) = a(:,j)
a(:,j) = tmp a(:,j) = tmp
endif cross end if cross
enddo end do
end subroutine qsort_partition end subroutine qsort_partition
@ -216,7 +216,7 @@ pure function math_expand(what,how)
do i = 1, size(how) do i = 1, size(how)
math_expand(sum(how(1:i-1))+1:sum(how(1:i))) = what(mod(i-1,size(what))+1) math_expand(sum(how(1:i-1))+1:sum(how(1:i))) = what(mod(i-1,size(what))+1)
enddo end do
end function math_expand end function math_expand
@ -251,7 +251,7 @@ pure function math_eye(d)
math_eye = 0.0_pReal math_eye = 0.0_pReal
do i=1,d do i=1,d
math_eye(i,i) = 1.0_pReal math_eye(i,i) = 1.0_pReal
enddo end do
end function math_eye end function math_eye
@ -270,7 +270,7 @@ pure function math_identity4th()
#ifndef __INTEL_COMPILER #ifndef __INTEL_COMPILER
do concurrent(i=1:3, j=1:3, k=1:3, l=1:3) do concurrent(i=1:3, j=1:3, k=1:3, l=1:3)
math_identity4th(i,j,k,l) = 0.5_pReal*(math_I3(i,k)*math_I3(j,l)+math_I3(i,l)*math_I3(j,k)) math_identity4th(i,j,k,l) = 0.5_pReal*(math_I3(i,k)*math_I3(j,l)+math_I3(i,l)*math_I3(j,k))
enddo end do
#else #else
forall(i=1:3, j=1:3, k=1:3, l=1:3) & forall(i=1:3, j=1:3, k=1:3, l=1:3) &
math_identity4th(i,j,k,l) = 0.5_pReal*(math_I3(i,k)*math_I3(j,l)+math_I3(i,l)*math_I3(j,k)) math_identity4th(i,j,k,l) = 0.5_pReal*(math_I3(i,k)*math_I3(j,l)+math_I3(i,l)*math_I3(j,k))
@ -298,7 +298,7 @@ real(pReal) pure function math_LeviCivita(i,j,k)
math_LeviCivita = -1.0_pReal math_LeviCivita = -1.0_pReal
else else
math_LeviCivita = 0.0_pReal math_LeviCivita = 0.0_pReal
endif end if
end function math_LeviCivita end function math_LeviCivita
@ -348,7 +348,7 @@ pure function math_outer(A,B)
#ifndef __INTEL_COMPILER #ifndef __INTEL_COMPILER
do concurrent(i=1:size(A,1), j=1:size(B,1)) do concurrent(i=1:size(A,1), j=1:size(B,1))
math_outer(i,j) = A(i)*B(j) math_outer(i,j) = A(i)*B(j)
enddo end do
#else #else
forall(i=1:size(A,1), j=1:size(B,1)) math_outer(i,j) = A(i)*B(j) forall(i=1:size(A,1), j=1:size(B,1)) math_outer(i,j) = A(i)*B(j)
#endif #endif
@ -398,7 +398,7 @@ pure function math_mul3333xx33(A,B)
#ifndef __INTEL_COMPILER #ifndef __INTEL_COMPILER
do concurrent(i=1:3, j=1:3) do concurrent(i=1:3, j=1:3)
math_mul3333xx33(i,j) = sum(A(i,j,1:3,1:3)*B(1:3,1:3)) math_mul3333xx33(i,j) = sum(A(i,j,1:3,1:3)*B(1:3,1:3))
enddo end do
#else #else
forall (i=1:3, j=1:3) math_mul3333xx33(i,j) = sum(A(i,j,1:3,1:3)*B(1:3,1:3)) forall (i=1:3, j=1:3) math_mul3333xx33(i,j) = sum(A(i,j,1:3,1:3)*B(1:3,1:3))
#endif #endif
@ -421,7 +421,7 @@ pure function math_mul3333xx3333(A,B)
#ifndef __INTEL_COMPILER #ifndef __INTEL_COMPILER
do concurrent(i=1:3, j=1:3, k=1:3, l=1:3) do concurrent(i=1:3, j=1:3, k=1:3, l=1:3)
math_mul3333xx3333(i,j,k,l) = sum(A(i,j,1:3,1:3)*B(1:3,1:3,k,l)) math_mul3333xx3333(i,j,k,l) = sum(A(i,j,1:3,1:3)*B(1:3,1:3,k,l))
enddo end do
#else #else
forall(i=1:3, j=1:3, k=1:3, l=1:3) math_mul3333xx3333(i,j,k,l) = sum(A(i,j,1:3,1:3)*B(1:3,1:3,k,l)) forall(i=1:3, j=1:3, k=1:3, l=1:3) math_mul3333xx3333(i,j,k,l) = sum(A(i,j,1:3,1:3)*B(1:3,1:3,k,l))
#endif #endif
@ -446,7 +446,7 @@ pure function math_exp33(A,n)
n_ = n n_ = n
else else
n_ = 5 n_ = 5
endif end if
invFac = 1.0_pReal ! 0! invFac = 1.0_pReal ! 0!
B = math_I3 B = math_I3
@ -456,7 +456,7 @@ pure function math_exp33(A,n)
invFac = invFac/real(i,pReal) ! invfac = 1/(i!) invFac = invFac/real(i,pReal) ! invfac = 1/(i!)
B = matmul(B,A) B = matmul(B,A)
math_exp33 = math_exp33 + invFac*B ! exp = SUM (A^i)/(i!) math_exp33 = math_exp33 + invFac*B ! exp = SUM (A^i)/(i!)
enddo end do
end function math_exp33 end function math_exp33
@ -514,7 +514,7 @@ pure subroutine math_invert33(InvA, DetA, error, A)
InvA = InvA/DetA InvA = InvA/DetA
error = .false. error = .false.
endif end if
end subroutine math_invert33 end subroutine math_invert33
@ -541,7 +541,7 @@ pure function math_invSym3333(A)
error stop 'matrix inversion error' error stop 'matrix inversion error'
else else
math_invSym3333 = math_66toSym3333(temp66) math_invSym3333 = math_66toSym3333(temp66)
endif end if
end function math_invSym3333 end function math_invSym3333
@ -696,7 +696,7 @@ pure function math_9to33(v9)
do i = 1, 9 do i = 1, 9
math_9to33(MAPPLAIN(1,i),MAPPLAIN(2,i)) = v9(i) math_9to33(MAPPLAIN(1,i),MAPPLAIN(2,i)) = v9(i)
enddo end do
end function math_9to33 end function math_9to33
@ -721,7 +721,7 @@ pure function math_sym33to6(m33,weighted)
w = merge(NRMMANDEL,1.0_pReal,weighted) w = merge(NRMMANDEL,1.0_pReal,weighted)
else else
w = NRMMANDEL w = NRMMANDEL
endif end if
math_sym33to6 = [(w(i)*m33(MAPNYE(1,i),MAPNYE(2,i)),i=1,6)] math_sym33to6 = [(w(i)*m33(MAPNYE(1,i),MAPNYE(2,i)),i=1,6)]
@ -748,12 +748,12 @@ pure function math_6toSym33(v6,weighted)
w = merge(INVNRMMANDEL,1.0_pReal,weighted) w = merge(INVNRMMANDEL,1.0_pReal,weighted)
else else
w = INVNRMMANDEL w = INVNRMMANDEL
endif end if
do i=1,6 do i=1,6
math_6toSym33(MAPNYE(1,i),MAPNYE(2,i)) = w(i)*v6(i) math_6toSym33(MAPNYE(1,i),MAPNYE(2,i)) = w(i)*v6(i)
math_6toSym33(MAPNYE(2,i),MAPNYE(1,i)) = w(i)*v6(i) math_6toSym33(MAPNYE(2,i),MAPNYE(1,i)) = w(i)*v6(i)
enddo end do
end function math_6toSym33 end function math_6toSym33
@ -772,7 +772,7 @@ pure function math_3333to99(m3333)
#ifndef __INTEL_COMPILER #ifndef __INTEL_COMPILER
do concurrent(i=1:9, j=1:9) do concurrent(i=1:9, j=1:9)
math_3333to99(i,j) = m3333(MAPPLAIN(1,i),MAPPLAIN(2,i),MAPPLAIN(1,j),MAPPLAIN(2,j)) math_3333to99(i,j) = m3333(MAPPLAIN(1,i),MAPPLAIN(2,i),MAPPLAIN(1,j),MAPPLAIN(2,j))
enddo end do
#else #else
forall(i=1:9, j=1:9) math_3333to99(i,j) = m3333(MAPPLAIN(1,i),MAPPLAIN(2,i),MAPPLAIN(1,j),MAPPLAIN(2,j)) forall(i=1:9, j=1:9) math_3333to99(i,j) = m3333(MAPPLAIN(1,i),MAPPLAIN(2,i),MAPPLAIN(1,j),MAPPLAIN(2,j))
#endif #endif
@ -794,7 +794,7 @@ pure function math_99to3333(m99)
#ifndef __INTEL_COMPILER #ifndef __INTEL_COMPILER
do concurrent(i=1:9, j=1:9) do concurrent(i=1:9, j=1:9)
math_99to3333(MAPPLAIN(1,i),MAPPLAIN(2,i),MAPPLAIN(1,j),MAPPLAIN(2,j)) = m99(i,j) math_99to3333(MAPPLAIN(1,i),MAPPLAIN(2,i),MAPPLAIN(1,j),MAPPLAIN(2,j)) = m99(i,j)
enddo end do
#else #else
forall(i=1:9, j=1:9) math_99to3333(MAPPLAIN(1,i),MAPPLAIN(2,i),MAPPLAIN(1,j),MAPPLAIN(2,j)) = m99(i,j) forall(i=1:9, j=1:9) math_99to3333(MAPPLAIN(1,i),MAPPLAIN(2,i),MAPPLAIN(1,j),MAPPLAIN(2,j)) = m99(i,j)
#endif #endif
@ -827,7 +827,7 @@ pure function math_sym3333to66(m3333,weighted)
#ifndef __INTEL_COMPILER #ifndef __INTEL_COMPILER
do concurrent(i=1:6, j=1:6) do concurrent(i=1:6, j=1:6)
math_sym3333to66(i,j) = w(i)*w(j)*m3333(MAPNYE(1,i),MAPNYE(2,i),MAPNYE(1,j),MAPNYE(2,j)) math_sym3333to66(i,j) = w(i)*w(j)*m3333(MAPNYE(1,i),MAPNYE(2,i),MAPNYE(1,j),MAPNYE(2,j))
enddo end do
#else #else
forall(i=1:6, j=1:6) math_sym3333to66(i,j) = w(i)*w(j)*m3333(MAPNYE(1,i),MAPNYE(2,i),MAPNYE(1,j),MAPNYE(2,j)) forall(i=1:6, j=1:6) math_sym3333to66(i,j) = w(i)*w(j)*m3333(MAPNYE(1,i),MAPNYE(2,i),MAPNYE(1,j),MAPNYE(2,j))
#endif #endif
@ -1080,8 +1080,8 @@ pure subroutine math_eigh33(w,v,m)
else fallback2 else fallback2
v(1:3,2) = v(1:3, 2) / norm v(1:3,2) = v(1:3, 2) / norm
v(1:3,3) = math_cross(v(1:3,1),v(1:3,2)) v(1:3,3) = math_cross(v(1:3,1),v(1:3,2))
endif fallback2 end if fallback2
endif fallback1 end if fallback1
end subroutine math_eigh33 end subroutine math_eigh33
@ -1120,7 +1120,7 @@ pure function math_rotationalPart(F) result(R)
lambda = sqrt(math_clip(lambda,0.0_pReal)/3.0_pReal) lambda = sqrt(math_clip(lambda,0.0_pReal)/3.0_pReal)
else else
lambda = sqrt(I_C(1)/3.0_pReal) lambda = sqrt(I_C(1)/3.0_pReal)
endif end if
I_U = [sum(lambda), lambda(1)*lambda(2)+lambda(2)*lambda(3)+lambda(3)*lambda(1), product(lambda)] I_U = [sum(lambda), lambda(1)*lambda(2)+lambda(2)*lambda(3)+lambda(3)*lambda(1), product(lambda)]
@ -1188,7 +1188,7 @@ pure function math_eigvalsh33(m)
cos((phi+2.0_pReal*TAU)/3.0_pReal) & cos((phi+2.0_pReal*TAU)/3.0_pReal) &
] & ] &
+ I(1)/3.0_pReal + I(1)/3.0_pReal
endif end if
end function math_eigvalsh33 end function math_eigvalsh33
@ -1238,7 +1238,7 @@ integer pure function math_binomial(n,k)
do i = 1, k_ do i = 1, k_
math_binomial = (math_binomial * n_)/i math_binomial = (math_binomial * n_)/i
n_ = n_ -1 n_ = n_ -1
enddo end do
end function math_binomial end function math_binomial
@ -1302,7 +1302,7 @@ real(pReal) pure elemental function math_clip(a, left, right)
if (present(right)) math_clip = min(right,math_clip) if (present(right)) math_clip = min(right,math_clip)
if (present(left) .and. present(right)) then if (present(left) .and. present(right)) then
if (left>right) error stop 'left > right' if (left>right) error stop 'left > right'
endif end if
end function math_clip end function math_clip
@ -1402,7 +1402,7 @@ subroutine selfTest()
do while(abs(math_det33(t33))<1.0e-9_pReal) do while(abs(math_det33(t33))<1.0e-9_pReal)
call random_number(t33) call random_number(t33)
enddo end do
if (any(dNeq0(matmul(t33,math_inv33(t33)) - math_eye(3),tol=1.0e-9_pReal))) & if (any(dNeq0(matmul(t33,math_inv33(t33)) - math_eye(3),tol=1.0e-9_pReal))) &
error stop 'math_inv33' error stop 'math_inv33'
@ -1418,7 +1418,7 @@ subroutine selfTest()
do while(math_det33(t33)<1.0e-2_pReal) ! O(det(F)) = 1 do while(math_det33(t33)<1.0e-2_pReal) ! O(det(F)) = 1
call random_number(t33) call random_number(t33)
enddo end do
t33_2 = math_rotationalPart(transpose(t33)) t33_2 = math_rotationalPart(transpose(t33))
t33 = math_rotationalPart(t33) t33 = math_rotationalPart(t33)
if (any(dNeq0(matmul(t33_2,t33) - math_I3,tol=1.0e-10_pReal))) & if (any(dNeq0(matmul(t33_2,t33) - math_I3,tol=1.0e-10_pReal))) &

View File

@ -97,7 +97,7 @@ subroutine parallelization_init
open(OUTPUT_UNIT,file='/dev/null',status='replace') ! close() alone will leave some temp files in cwd open(OUTPUT_UNIT,file='/dev/null',status='replace') ! close() alone will leave some temp files in cwd
else else
open(OUTPUT_UNIT,encoding='UTF-8') ! for special characters in output open(OUTPUT_UNIT,encoding='UTF-8') ! for special characters in output
endif end if
#endif #endif
print'(/,1x,a)', '<<<+- parallelization init -+>>>' print'(/,1x,a)', '<<<+- parallelization init -+>>>'
@ -142,8 +142,8 @@ subroutine parallelization_init
!$ if (OMP_NUM_THREADS < 1_pI32) then !$ if (OMP_NUM_THREADS < 1_pI32) then
!$ print'(1x,a)', 'Invalid OMP_NUM_THREADS: "'//trim(NumThreadsString)//'", using default' !$ print'(1x,a)', 'Invalid OMP_NUM_THREADS: "'//trim(NumThreadsString)//'", using default'
!$ OMP_NUM_THREADS = 4_pI32 !$ OMP_NUM_THREADS = 4_pI32
!$ endif !$ end if
!$ endif !$ end if
!$ print'(1x,a,i0)', 'OMP_NUM_THREADS: ',OMP_NUM_THREADS !$ print'(1x,a,i0)', 'OMP_NUM_THREADS: ',OMP_NUM_THREADS
!$ call omp_set_num_threads(OMP_NUM_THREADS) !$ call omp_set_num_threads(OMP_NUM_THREADS)

View File

@ -255,7 +255,7 @@ module subroutine mechanical_init(phases)
#else #else
output_mechanical(ph)%label = mech%get_as1dString('output',defaultVal=emptyStringArray) output_mechanical(ph)%label = mech%get_as1dString('output',defaultVal=emptyStringArray)
#endif #endif
enddo end do
do ce = 1, size(material_phaseID,2) do ce = 1, size(material_phaseID,2)
ma = discretization_materialAt((ce-1)/discretization_nIPs+1) ma = discretization_materialAt((ce-1)/discretization_nIPs+1)
@ -267,14 +267,14 @@ module subroutine mechanical_init(phases)
phase_mechanical_Fe(ph)%data(1:3,1:3,en) = matmul(material_V_e_0(ma)%data(1:3,1:3,co), & phase_mechanical_Fe(ph)%data(1:3,1:3,en) = matmul(material_V_e_0(ma)%data(1:3,1:3,co), &
transpose(phase_mechanical_Fp(ph)%data(1:3,1:3,en))) transpose(phase_mechanical_Fp(ph)%data(1:3,1:3,en)))
phase_mechanical_Fi(ph)%data(1:3,1:3,en) = material_O_0(ma)%data(co)%rotate(math_inv33(material_V_e_0(ma)%data(1:3,1:3,co))) phase_mechanical_Fi(ph)%data(1:3,1:3,en) = material_O_0(ma)%data(co)%rotate(math_inv33(material_V_e_0(ma)%data(1:3,1:3,co)))
enddo end do
enddo end do
do ph = 1, phases%length do ph = 1, phases%length
phase_mechanical_F0(ph)%data = phase_mechanical_F(ph)%data phase_mechanical_F0(ph)%data = phase_mechanical_F(ph)%data
phase_mechanical_Fp0(ph)%data = phase_mechanical_Fp(ph)%data phase_mechanical_Fp0(ph)%data = phase_mechanical_Fp(ph)%data
phase_mechanical_Fi0(ph)%data = phase_mechanical_Fi(ph)%data phase_mechanical_Fi0(ph)%data = phase_mechanical_Fi(ph)%data
enddo end do
call elastic_init(phases) call elastic_init(phases)
@ -284,7 +284,7 @@ module subroutine mechanical_init(phases)
call plastic_init() call plastic_init()
do ph = 1,phases%length do ph = 1,phases%length
plasticState(ph)%state0 = plasticState(ph)%state plasticState(ph)%state0 = plasticState(ph)%state
enddo end do
num_crystallite => config_numerics%get('crystallite',defaultVal=emptyDict) num_crystallite => config_numerics%get('crystallite',defaultVal=emptyDict)
@ -473,25 +473,25 @@ function integrateStress(F,subFp0,subFi0,Delta_t,ph,en) result(broken)
Lpguess = Lpguess_old & Lpguess = Lpguess_old &
+ deltaLp * stepLengthLp + deltaLp * stepLengthLp
cycle LpLoop cycle LpLoop
endif end if
calculateJacobiLp: if (mod(jacoCounterLp, num%iJacoLpresiduum) == 0) then calculateJacobiLp: if (mod(jacoCounterLp, num%iJacoLpresiduum) == 0) then
jacoCounterLp = jacoCounterLp + 1 jacoCounterLp = jacoCounterLp + 1
do o=1,3; do p=1,3 do o=1,3; do p=1,3
dFe_dLp(o,1:3,p,1:3) = - Delta_t * A(o,p)*transpose(invFi_new) ! dFe_dLp(i,j,k,l) = -Delta_t * A(i,k) invFi(l,j) dFe_dLp(o,1:3,p,1:3) = - Delta_t * A(o,p)*transpose(invFi_new) ! dFe_dLp(i,j,k,l) = -Delta_t * A(i,k) invFi(l,j)
enddo; enddo end do; end do
dRLp_dLp = math_eye(9) & dRLp_dLp = math_eye(9) &
- math_3333to99(math_mul3333xx3333(math_mul3333xx3333(dLp_dS,dS_dFe),dFe_dLp)) - math_3333to99(math_mul3333xx3333(math_mul3333xx3333(dLp_dS,dS_dFe),dFe_dLp))
temp_9 = math_33to9(residuumLp) temp_9 = math_33to9(residuumLp)
call dgesv(9,1,dRLp_dLp,9,devNull_9,temp_9,9,ierr) ! solve dRLp/dLp * delta Lp = -res for delta Lp call dgesv(9,1,dRLp_dLp,9,devNull_9,temp_9,9,ierr) ! solve dRLp/dLp * delta Lp = -res for delta Lp
if (ierr /= 0) return ! error if (ierr /= 0) return ! error
deltaLp = - math_9to33(temp_9) deltaLp = - math_9to33(temp_9)
endif calculateJacobiLp end if calculateJacobiLp
Lpguess = Lpguess & Lpguess = Lpguess &
+ deltaLp * steplengthLp + deltaLp * steplengthLp
enddo LpLoop end do LpLoop
call phase_LiAndItsTangents(Li_constitutive, dLi_dS, dLi_dFi, & call phase_LiAndItsTangents(Li_constitutive, dLi_dS, dLi_dFi, &
S, Fi_new, ph,en) S, Fi_new, ph,en)
@ -513,7 +513,7 @@ function integrateStress(F,subFp0,subFi0,Delta_t,ph,en) result(broken)
Liguess = Liguess_old & Liguess = Liguess_old &
+ deltaLi * steplengthLi + deltaLi * steplengthLi
cycle LiLoop cycle LiLoop
endif end if
calculateJacobiLi: if (mod(jacoCounterLi, num%iJacoLpresiduum) == 0) then calculateJacobiLi: if (mod(jacoCounterLi, num%iJacoLpresiduum) == 0) then
jacoCounterLi = jacoCounterLi + 1 jacoCounterLi = jacoCounterLi + 1
@ -522,10 +522,10 @@ function integrateStress(F,subFp0,subFi0,Delta_t,ph,en) result(broken)
do o=1,3; do p=1,3 do o=1,3; do p=1,3
dFe_dLi(1:3,o,1:3,p) = -Delta_t*math_I3(o,p)*temp_33 ! dFe_dLp(i,j,k,l) = -Delta_t * A(i,k) invFi(l,j) dFe_dLi(1:3,o,1:3,p) = -Delta_t*math_I3(o,p)*temp_33 ! dFe_dLp(i,j,k,l) = -Delta_t * A(i,k) invFi(l,j)
dFi_dLi(1:3,o,1:3,p) = -Delta_t*math_I3(o,p)*invFi_current dFi_dLi(1:3,o,1:3,p) = -Delta_t*math_I3(o,p)*invFi_current
enddo; enddo end do; end do
do o=1,3; do p=1,3 do o=1,3; do p=1,3
dFi_dLi(1:3,1:3,o,p) = matmul(matmul(Fi_new,dFi_dLi(1:3,1:3,o,p)),Fi_new) dFi_dLi(1:3,1:3,o,p) = matmul(matmul(Fi_new,dFi_dLi(1:3,1:3,o,p)),Fi_new)
enddo; enddo end do; end do
dRLi_dLi = math_eye(9) & dRLi_dLi = math_eye(9) &
- math_3333to99(math_mul3333xx3333(dLi_dS, math_mul3333xx3333(dS_dFe, dFe_dLi) & - math_3333to99(math_mul3333xx3333(dLi_dS, math_mul3333xx3333(dS_dFe, dFe_dLi) &
+ math_mul3333xx3333(dS_dFi, dFi_dLi))) & + math_mul3333xx3333(dS_dFi, dFi_dLi))) &
@ -534,11 +534,11 @@ function integrateStress(F,subFp0,subFi0,Delta_t,ph,en) result(broken)
call dgesv(9,1,dRLi_dLi,9,devNull_9,temp_9,9,ierr) ! solve dRLi/dLp * delta Li = -res for delta Li call dgesv(9,1,dRLi_dLi,9,devNull_9,temp_9,9,ierr) ! solve dRLi/dLp * delta Li = -res for delta Li
if (ierr /= 0) return ! error if (ierr /= 0) return ! error
deltaLi = - math_9to33(temp_9) deltaLi = - math_9to33(temp_9)
endif calculateJacobiLi end if calculateJacobiLi
Liguess = Liguess & Liguess = Liguess &
+ deltaLi * steplengthLi + deltaLi * steplengthLi
enddo LiLoop end do LiLoop
invFp_new = matmul(invFp_current,B) invFp_new = matmul(invFp_current,B)
call math_invert33(Fp_new,devNull,error,invFp_new) call math_invert33(Fp_new,devNull,error,invFp_new)
@ -613,9 +613,9 @@ function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en) result(b
if (converged(r,plasticState(ph)%state(1:sizeDotState,en),plasticState(ph)%atol(1:sizeDotState))) then if (converged(r,plasticState(ph)%state(1:sizeDotState,en),plasticState(ph)%atol(1:sizeDotState))) then
broken = plastic_deltaState(ph,en) broken = plastic_deltaState(ph,en)
exit iteration exit iteration
endif end if
enddo iteration end do iteration
contains contains
@ -638,7 +638,7 @@ function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en) result(b
damper = 0.75_pReal + 0.25_pReal * tanh(2.0_pReal + 4.0_pReal * dot_prod12 / dot_prod22) damper = 0.75_pReal + 0.25_pReal * tanh(2.0_pReal + 4.0_pReal * dot_prod12 / dot_prod22)
else else
damper = 1.0_pReal damper = 1.0_pReal
endif end if
end function damper end function damper
@ -844,7 +844,7 @@ function integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en,A,B,C,DB)
#else #else
dotState = IEEE_FMA(A(n,stage),plastic_RKdotState(1:sizeDotState,n),dotState) dotState = IEEE_FMA(A(n,stage),plastic_RKdotState(1:sizeDotState,n),dotState)
#endif #endif
enddo end do
#ifndef __INTEL_LLVM_COMPILER #ifndef __INTEL_LLVM_COMPILER
plasticState(ph)%state(1:sizeDotState,en) = subState0 + dotState*Delta_t plasticState(ph)%state(1:sizeDotState,en) = subState0 + dotState*Delta_t
@ -858,7 +858,7 @@ function integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en,A,B,C,DB)
dotState = plastic_dotState(Delta_t*C(stage), ph,en) dotState = plastic_dotState(Delta_t*C(stage), ph,en)
if (any(IEEE_is_NaN(dotState))) exit if (any(IEEE_is_NaN(dotState))) exit
enddo end do
if(broken) return if(broken) return
@ -950,7 +950,7 @@ subroutine results(group,ph)
do i = 1, size(dataset,1) do i = 1, size(dataset,1)
to_quaternion(:,i) = dataset(i)%asQuaternion() to_quaternion(:,i) = dataset(i)%asQuaternion()
enddo end do
end function to_quaternion end function to_quaternion
@ -974,7 +974,7 @@ module subroutine mechanical_forward()
phase_mechanical_Lp0(ph) = phase_mechanical_Lp(ph) phase_mechanical_Lp0(ph) = phase_mechanical_Lp(ph)
phase_mechanical_S0(ph) = phase_mechanical_S(ph) phase_mechanical_S0(ph) = phase_mechanical_S(ph)
plasticState(ph)%state0 = plasticState(ph)%state plasticState(ph)%state0 = plasticState(ph)%state
enddo end do
end subroutine mechanical_forward end subroutine mechanical_forward
@ -1037,7 +1037,7 @@ module function phase_mechanical_constitutive(Delta_t,co,ce) result(converged_)
subFp0 = phase_mechanical_Fp(ph)%data(1:3,1:3,en) subFp0 = phase_mechanical_Fp(ph)%data(1:3,1:3,en)
subFi0 = phase_mechanical_Fi(ph)%data(1:3,1:3,en) subFi0 = phase_mechanical_Fi(ph)%data(1:3,1:3,en)
subState0 = plasticState(ph)%state(:,en) subState0 = plasticState(ph)%state(:,en)
endif end if
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! cut back (reduced time and restore) ! cut back (reduced time and restore)
else else
@ -1048,10 +1048,10 @@ module function phase_mechanical_constitutive(Delta_t,co,ce) result(converged_)
if (subStep < 1.0_pReal) then ! actual (not initial) cutback if (subStep < 1.0_pReal) then ! actual (not initial) cutback
phase_mechanical_Lp(ph)%data(1:3,1:3,en) = subLp0 phase_mechanical_Lp(ph)%data(1:3,1:3,en) = subLp0
phase_mechanical_Li(ph)%data(1:3,1:3,en) = subLi0 phase_mechanical_Li(ph)%data(1:3,1:3,en) = subLi0
endif end if
plasticState(ph)%state(:,en) = subState0 plasticState(ph)%state(:,en) = subState0
todo = subStep > num%subStepMinCryst ! still on track or already done (beyond repair) todo = subStep > num%subStepMinCryst ! still on track or already done (beyond repair)
endif end if
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! prepare for integration ! prepare for integration
@ -1060,9 +1060,9 @@ module function phase_mechanical_constitutive(Delta_t,co,ce) result(converged_)
subF = subF0 & subF = subF0 &
+ subStep * (phase_mechanical_F(ph)%data(1:3,1:3,en) - phase_mechanical_F0(ph)%data(1:3,1:3,en)) + subStep * (phase_mechanical_F(ph)%data(1:3,1:3,en) - phase_mechanical_F0(ph)%data(1:3,1:3,en))
converged_ = .not. integrateState(subF0,subF,subFp0,subFi0,subState0(1:sizeDotState),subStep * Delta_t,ph,en) converged_ = .not. integrateState(subF0,subF,subFp0,subFi0,subState0(1:sizeDotState),subStep * Delta_t,ph,en)
endif end if
enddo cutbackLooping end do cutbackLooping
end function phase_mechanical_constitutive end function phase_mechanical_constitutive
@ -1086,14 +1086,14 @@ module subroutine mechanical_restore(ce,includeL)
if (includeL) then if (includeL) then
phase_mechanical_Lp(ph)%data(1:3,1:3,en) = phase_mechanical_Lp0(ph)%data(1:3,1:3,en) phase_mechanical_Lp(ph)%data(1:3,1:3,en) = phase_mechanical_Lp0(ph)%data(1:3,1:3,en)
phase_mechanical_Li(ph)%data(1:3,1:3,en) = phase_mechanical_Li0(ph)%data(1:3,1:3,en) phase_mechanical_Li(ph)%data(1:3,1:3,en) = phase_mechanical_Li0(ph)%data(1:3,1:3,en)
endif ! maybe protecting everything from overwriting makes more sense end if ! maybe protecting everything from overwriting makes more sense
phase_mechanical_Fp(ph)%data(1:3,1:3,en) = phase_mechanical_Fp0(ph)%data(1:3,1:3,en) phase_mechanical_Fp(ph)%data(1:3,1:3,en) = phase_mechanical_Fp0(ph)%data(1:3,1:3,en)
phase_mechanical_Fi(ph)%data(1:3,1:3,en) = phase_mechanical_Fi0(ph)%data(1:3,1:3,en) phase_mechanical_Fi(ph)%data(1:3,1:3,en) = phase_mechanical_Fi0(ph)%data(1:3,1:3,en)
phase_mechanical_S(ph)%data(1:3,1:3,en) = phase_mechanical_S0(ph)%data(1:3,1:3,en) phase_mechanical_S(ph)%data(1:3,1:3,en) = phase_mechanical_S0(ph)%data(1:3,1:3,en)
plasticState(ph)%state(:,en) = plasticState(ph)%State0(:,en) plasticState(ph)%state(:,en) = plasticState(ph)%State0(:,en)
enddo end do
end subroutine mechanical_restore end subroutine mechanical_restore
@ -1164,7 +1164,7 @@ module function phase_mechanical_dPdF(Delta_t,co,ce) result(dPdF)
lhs_3333(1:3,o,1:3,p) = IEEE_FMA(invFi,invFi(p,o),lhs_3333(1:3,o,1:3,p)) lhs_3333(1:3,o,1:3,p) = IEEE_FMA(invFi,invFi(p,o),lhs_3333(1:3,o,1:3,p))
rhs_3333(1:3,1:3,o,p) = IEEE_FMA(matmul(invSubFi0,dLidS(1:3,1:3,o,p)),-Delta_t,rhs_3333(1:3,1:3,o,p)) rhs_3333(1:3,1:3,o,p) = IEEE_FMA(matmul(invSubFi0,dLidS(1:3,1:3,o,p)),-Delta_t,rhs_3333(1:3,1:3,o,p))
#endif #endif
enddo; enddo end do; end do
call math_invert(temp_99,error,math_3333to99(lhs_3333)) call math_invert(temp_99,error,math_3333to99(lhs_3333))
if (error) then if (error) then
call IO_warning(600,'inversion error in analytic tangent calculation', & call IO_warning(600,'inversion error in analytic tangent calculation', &
@ -1172,9 +1172,9 @@ module function phase_mechanical_dPdF(Delta_t,co,ce) result(dPdF)
dFidS = 0.0_pReal dFidS = 0.0_pReal
else else
dFidS = math_mul3333xx3333(math_99to3333(temp_99),rhs_3333) dFidS = math_mul3333xx3333(math_99to3333(temp_99),rhs_3333)
endif end if
dLidS = math_mul3333xx3333(dLidFi,dFidS) + dLidS dLidS = math_mul3333xx3333(dLidFi,dFidS) + dLidS
endif end if
call plastic_LpAndItsTangents(devNull,dLpdS,dLpdFi, & call plastic_LpAndItsTangents(devNull,dLpdS,dLpdFi, &
phase_mechanical_S(ph)%data(1:3,1:3,en), & phase_mechanical_S(ph)%data(1:3,1:3,en), &
@ -1191,7 +1191,7 @@ module function phase_mechanical_dPdF(Delta_t,co,ce) result(dPdF)
rhs_3333(p,o,1:3,1:3) = matmul(dSdFe(p,o,1:3,1:3),temp_33_1) rhs_3333(p,o,1:3,1:3) = matmul(dSdFe(p,o,1:3,1:3),temp_33_1)
temp_3333(1:3,1:3,p,o) = matmul(matmul(temp_33_2,dLpdS(1:3,1:3,p,o)), invFi) & temp_3333(1:3,1:3,p,o) = matmul(matmul(temp_33_2,dLpdS(1:3,1:3,p,o)), invFi) &
+ matmul(temp_33_3,dLidS(1:3,1:3,p,o)) + matmul(temp_33_3,dLidS(1:3,1:3,p,o))
enddo; enddo end do; end do
#ifndef __INTEL_LLVM_COMPILER #ifndef __INTEL_LLVM_COMPILER
lhs_3333 = math_mul3333xx3333(dSdFe,temp_3333) * Delta_t & lhs_3333 = math_mul3333xx3333(dSdFe,temp_3333) * Delta_t &
+ math_mul3333xx3333(dSdFi,dFidS) + math_mul3333xx3333(dSdFi,dFidS)
@ -1206,14 +1206,14 @@ module function phase_mechanical_dPdF(Delta_t,co,ce) result(dPdF)
dSdF = rhs_3333 dSdF = rhs_3333
else else
dSdF = math_mul3333xx3333(math_99to3333(temp_99),rhs_3333) dSdF = math_mul3333xx3333(math_99to3333(temp_99),rhs_3333)
endif end if
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! calculate dFpinvdF ! calculate dFpinvdF
temp_3333 = math_mul3333xx3333(dLpdS,dSdF) temp_3333 = math_mul3333xx3333(dLpdS,dSdF)
do o=1,3; do p=1,3 do o=1,3; do p=1,3
dFpinvdF(1:3,1:3,p,o) = - matmul(invSubFp0, matmul(temp_3333(1:3,1:3,p,o),invFi)) * Delta_t dFpinvdF(1:3,1:3,p,o) = - matmul(invSubFp0, matmul(temp_3333(1:3,1:3,p,o),invFi)) * Delta_t
enddo; enddo end do; end do
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! assemble dPdF ! assemble dPdF
@ -1224,13 +1224,13 @@ module function phase_mechanical_dPdF(Delta_t,co,ce) result(dPdF)
dPdF = 0.0_pReal dPdF = 0.0_pReal
do p=1,3 do p=1,3
dPdF(p,1:3,p,1:3) = transpose(matmul(invFp,temp_33_1)) dPdF(p,1:3,p,1:3) = transpose(matmul(invFp,temp_33_1))
enddo end do
do o=1,3; do p=1,3 do o=1,3; do p=1,3
dPdF(1:3,1:3,p,o) = dPdF(1:3,1:3,p,o) & dPdF(1:3,1:3,p,o) = dPdF(1:3,1:3,p,o) &
+ matmul(matmul(phase_mechanical_F(ph)%data(1:3,1:3,en),dFpinvdF(1:3,1:3,p,o)),temp_33_1) & + matmul(matmul(phase_mechanical_F(ph)%data(1:3,1:3,en),dFpinvdF(1:3,1:3,p,o)),temp_33_1) &
+ matmul(matmul(temp_33_2,dSdF(1:3,1:3,p,o)),transpose(invFp)) & + matmul(matmul(temp_33_2,dSdF(1:3,1:3,p,o)),transpose(invFp)) &
+ matmul(temp_33_3,transpose(dFpinvdF(1:3,1:3,p,o))) + matmul(temp_33_3,transpose(dFpinvdF(1:3,1:3,p,o)))
enddo; enddo end do; end do
end function phase_mechanical_dPdF end function phase_mechanical_dPdF

View File

@ -61,7 +61,7 @@ module subroutine eigen_init(phases)
if (maxval(Nmodels) /= 0) then if (maxval(Nmodels) /= 0) then
where(thermalexpansion_init(maxval(Nmodels))) model = EIGEN_thermal_expansion_ID where(thermalexpansion_init(maxval(Nmodels))) model = EIGEN_thermal_expansion_ID
endif end if
allocate(model_damage(phases%length), source = EIGEN_UNDEFINED_ID) allocate(model_damage(phases%length), source = EIGEN_UNDEFINED_ID)

View File

@ -435,7 +435,7 @@ function plastic_active(plastic_label) result(active_plastic)
mech => phase%get('mechanical') mech => phase%get('mechanical')
pl => mech%get('plastic',defaultVal = emptyDict) pl => mech%get('plastic',defaultVal = emptyDict)
active_plastic(ph) = pl%get_asString('type',defaultVal='none') == plastic_label active_plastic(ph) = pl%get_asString('type',defaultVal='none') == plastic_label
enddo end do
end function plastic_active end function plastic_active

View File

@ -9,28 +9,28 @@
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
submodule(phase:plastic) dislotwin submodule(phase:plastic) dislotwin
real(pReal), parameter :: gamma_char_tr = sqrt(0.125_pReal) !< Characteristic shear for transformation
type :: tParameters type :: tParameters
real(pReal) :: & real(pReal) :: &
Q_cl = 1.0_pReal, & !< activation energy for dislocation climb Q_cl = 1.0_pReal, & !< activation energy for dislocation climb
omega = 1.0_pReal, & !< frequency factor for dislocation climb omega = 1.0_pReal, & !< frequency factor for dislocation climb
D = 1.0_pReal, & !< grain size D = 1.0_pReal, & !< grain size
p_sb = 1.0_pReal, & !< p-exponent in shear band velocity p_sb = 1.0_pReal, & !< p-exponent in shear band velocity
q_sb = 1.0_pReal, & !< q-exponent in shear band velocity q_sb = 1.0_pReal, & !< q-exponent in shear band velocity
i_tw = 1.0_pReal, & !< adjustment parameter to calculate MFP for twinning i_tw = 1.0_pReal, & !< adjustment parameter to calculate MFP for twinning
i_tr = 1.0_pReal, & !< adjustment parameter to calculate MFP for transformation i_tr = 1.0_pReal, & !< adjustment parameter to calculate MFP for transformation
L_tw = 1.0_pReal, & !< length of twin nuclei L_tw = 1.0_pReal, & !< length of twin nuclei
L_tr = 1.0_pReal, & !< length of trans nuclei L_tr = 1.0_pReal, & !< length of trans nuclei
x_c = 1.0_pReal, & !< critical distance for formation of twin/trans nucleus x_c = 1.0_pReal, & !< critical distance for formation of twin/trans nucleus
V_cs = 1.0_pReal, & !< cross slip volume V_cs = 1.0_pReal, & !< cross slip volume
xi_sb = 1.0_pReal, & !< value for shearband resistance tau_sb = 1.0_pReal, & !< value for shearband resistance
v_sb = 1.0_pReal, & !< value for shearband velocity_0 gamma_0_sb = 1.0_pReal, & !< value for shearband velocity_0
E_sb = 1.0_pReal, & !< activation energy for shear bands E_sb = 1.0_pReal, & !< activation energy for shear bands
h = 1.0_pReal, & !< stack height of hex nucleus h = 1.0_pReal, & !< stack height of hex nucleus
gamma_char_tr = sqrt(0.125_pReal), & !< Characteristic shear for transformation a_cF = 1.0_pReal, &
a_cF = 1.0_pReal, & cOverA_hP = 1.0_pReal, &
cOverA_hP = 1.0_pReal, & V_mol = 1.0_pReal, &
V_mol = 1.0_pReal, & rho = 1.0_pReal
rho = 1.0_pReal
type(tPolynomial) :: & type(tPolynomial) :: &
Gamma_sf, & !< stacking fault energy Gamma_sf, & !< stacking fault energy
Delta_G !< free energy difference between austensite and martensite Delta_G !< free energy difference between austensite and martensite
@ -331,18 +331,18 @@ module function plastic_dislotwin_init() result(myPlasticity)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! shearband related parameters ! shearband related parameters
prm%v_sb = pl%get_asFloat('v_sb',defaultVal=0.0_pReal) prm%gamma_0_sb = pl%get_asFloat('gamma_0_sb',defaultVal=0.0_pReal)
if (prm%v_sb > 0.0_pReal) then if (prm%gamma_0_sb > 0.0_pReal) then
prm%xi_sb = pl%get_asFloat('xi_sb') prm%tau_sb = pl%get_asFloat('tau_sb')
prm%E_sb = pl%get_asFloat('Q_sb') prm%E_sb = pl%get_asFloat('Q_sb')
prm%p_sb = pl%get_asFloat('p_sb') prm%p_sb = pl%get_asFloat('p_sb')
prm%q_sb = pl%get_asFloat('q_sb') prm%q_sb = pl%get_asFloat('q_sb')
! sanity checks ! sanity checks
if (prm%xi_sb < 0.0_pReal) extmsg = trim(extmsg)//' xi_sb' if (prm%tau_sb < 0.0_pReal) extmsg = trim(extmsg)//' tau_sb'
if (prm%E_sb < 0.0_pReal) extmsg = trim(extmsg)//' Q_sb' if (prm%E_sb < 0.0_pReal) extmsg = trim(extmsg)//' Q_sb'
if (prm%p_sb <= 0.0_pReal) extmsg = trim(extmsg)//' p_sb' if (prm%p_sb <= 0.0_pReal) extmsg = trim(extmsg)//' p_sb'
if (prm%q_sb <= 0.0_pReal) extmsg = trim(extmsg)//' q_sb' if (prm%q_sb <= 0.0_pReal) extmsg = trim(extmsg)//' q_sb'
end if end if
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -364,13 +364,13 @@ module function plastic_dislotwin_init() result(myPlasticity)
prm%h_sl_tw = lattice_interaction_SlipByTwin(N_sl,prm%N_tw,pl%get_as1dFloat('h_sl-tw'), & prm%h_sl_tw = lattice_interaction_SlipByTwin(N_sl,prm%N_tw,pl%get_as1dFloat('h_sl-tw'), &
phase_lattice(ph)) phase_lattice(ph))
if (prm%fccTwinTransNucleation .and. size(prm%N_tw) /= 1) extmsg = trim(extmsg)//' N_tw: nucleation' if (prm%fccTwinTransNucleation .and. size(prm%N_tw) /= 1) extmsg = trim(extmsg)//' N_tw: nucleation'
endif slipAndTwinActive end if slipAndTwinActive
slipAndTransActive: if (prm%sum_N_sl * prm%sum_N_tr > 0) then slipAndTransActive: if (prm%sum_N_sl * prm%sum_N_tr > 0) then
prm%h_sl_tr = lattice_interaction_SlipByTrans(N_sl,prm%N_tr,pl%get_as1dFloat('h_sl-tr'), & prm%h_sl_tr = lattice_interaction_SlipByTrans(N_sl,prm%N_tr,pl%get_as1dFloat('h_sl-tr'), &
phase_lattice(ph)) phase_lattice(ph))
if (prm%fccTwinTransNucleation .and. size(prm%N_tr) /= 1) extmsg = trim(extmsg)//' N_tr: nucleation' if (prm%fccTwinTransNucleation .and. size(prm%N_tr) /= 1) extmsg = trim(extmsg)//' N_tr: nucleation'
endif slipAndTransActive end if slipAndTransActive
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! allocate state arrays ! allocate state arrays
@ -569,7 +569,7 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en)
Lp = Lp * f_matrix Lp = Lp * f_matrix
dLp_dMp = dLp_dMp * f_matrix dLp_dMp = dLp_dMp * f_matrix
shearBandingContribution: if (dNeq0(prm%v_sb)) then shearBandingContribution: if (dNeq0(prm%gamma_0_sb)) then
E_kB_T = prm%E_sb/(K_B*T) E_kB_T = prm%E_sb/(K_B*T)
call math_eigh33(eigValues,eigVectors,Mp) ! is Mp symmetric by design? call math_eigh33(eigValues,eigVectors,Mp) ! is Mp symmetric by design?
@ -580,10 +580,10 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en)
tau = math_tensordot(Mp,P_sb) tau = math_tensordot(Mp,P_sb)
significantShearBandStress: if (abs(tau) > tol_math_check) then significantShearBandStress: if (abs(tau) > tol_math_check) then
StressRatio_p = (abs(tau)/prm%xi_sb)**prm%p_sb StressRatio_p = (abs(tau)/prm%tau_sb)**prm%p_sb
dot_gamma_sb = sign(prm%v_sb*exp(-E_kB_T*(1-StressRatio_p)**prm%q_sb), tau) dot_gamma_sb = sign(prm%gamma_0_sb*exp(-E_kB_T*(1-StressRatio_p)**prm%q_sb), tau)
ddot_gamma_dtau = abs(dot_gamma_sb)*E_kB_T*prm%p_sb*prm%q_sb/prm%xi_sb & ddot_gamma_dtau = abs(dot_gamma_sb)*E_kB_T*prm%p_sb*prm%q_sb/prm%tau_sb &
* (abs(tau)/prm%xi_sb)**(prm%p_sb-1.0_pReal) & * (abs(tau)/prm%tau_sb)**(prm%p_sb-1.0_pReal) &
* (1.0_pReal-StressRatio_p)**(prm%q_sb-1.0_pReal) * (1.0_pReal-StressRatio_p)**(prm%q_sb-1.0_pReal)
Lp = Lp + dot_gamma_sb * P_sb Lp = Lp + dot_gamma_sb * P_sb
@ -697,7 +697,7 @@ module function dislotwin_dotState(Mp,ph,en) result(dotState)
dot_f_tw = f_matrix*dot_gamma_tw/prm%gamma_char_tw dot_f_tw = f_matrix*dot_gamma_tw/prm%gamma_char_tw
if (prm%sum_N_tr > 0) call kinetics_tr(Mp,T,abs_dot_gamma_sl,ph,en,dot_gamma_tr) if (prm%sum_N_tr > 0) call kinetics_tr(Mp,T,abs_dot_gamma_sl,ph,en,dot_gamma_tr)
dot_f_tr = f_matrix*dot_gamma_tr/prm%gamma_char_tr dot_f_tr = f_matrix*dot_gamma_tr/gamma_char_tr
end associate end associate
@ -1026,9 +1026,9 @@ pure subroutine kinetics_tr(Mp,T,abs_dot_gamma_sl,ph,en,&
dP_ncs_dtau = prm%V_cs / (K_B * T) * (P_ncs - 1.0_pReal) dP_ncs_dtau = prm%V_cs / (K_B * T) * (P_ncs - 1.0_pReal)
V = PI/4.0_pReal*dst%Lambda_tr(i,en)**2*prm%t_tr(i) V = PI/4.0_pReal*dst%Lambda_tr(i,en)**2*prm%t_tr(i)
dot_gamma_tr(i) = V*dot_N_0*P_ncs*P*prm%gamma_char_tr dot_gamma_tr(i) = V*dot_N_0*P_ncs*P*gamma_char_tr
if (present(ddot_gamma_dtau_tr)) & if (present(ddot_gamma_dtau_tr)) &
ddot_gamma_dtau_tr(i) = V*dot_N_0*(P*dP_ncs_dtau + P_ncs*dP_dtau)*prm%gamma_char_tr ddot_gamma_dtau_tr(i) = V*dot_N_0*(P*dP_ncs_dtau + P_ncs*dP_dtau)*gamma_char_tr
else else
dot_gamma_tr(i) = 0.0_pReal dot_gamma_tr(i) = 0.0_pReal
if (present(ddot_gamma_dtau_tr)) ddot_gamma_dtau_tr(i) = 0.0_pReal if (present(ddot_gamma_dtau_tr)) ddot_gamma_dtau_tr(i) = 0.0_pReal

View File

@ -122,25 +122,25 @@ module subroutine thermal_init(phases)
allocate(thermalstate(ph)%p(thermal_Nsources(ph))) allocate(thermalstate(ph)%p(thermal_Nsources(ph)))
enddo end do
allocate(thermal_source(maxval(thermal_Nsources),phases%length), source = THERMAL_UNDEFINED_ID) allocate(thermal_source(maxval(thermal_Nsources),phases%length), source = THERMAL_UNDEFINED_ID)
if (maxval(thermal_Nsources) /= 0) then if (maxval(thermal_Nsources) /= 0) then
where(dissipation_init (maxval(thermal_Nsources))) thermal_source = THERMAL_DISSIPATION_ID where(dissipation_init (maxval(thermal_Nsources))) thermal_source = THERMAL_DISSIPATION_ID
where(externalheat_init(maxval(thermal_Nsources))) thermal_source = THERMAL_EXTERNALHEAT_ID where(externalheat_init(maxval(thermal_Nsources))) thermal_source = THERMAL_EXTERNALHEAT_ID
endif end if
thermal_source_maxSizeDotState = 0 thermal_source_maxSizeDotState = 0
do ph = 1,phases%length do ph = 1,phases%length
do so = 1,thermal_Nsources(ph) do so = 1,thermal_Nsources(ph)
thermalState(ph)%p(so)%state = thermalState(ph)%p(so)%state0 thermalState(ph)%p(so)%state = thermalState(ph)%p(so)%state0
enddo end do
thermal_source_maxSizeDotState = max(thermal_source_maxSizeDotState, & thermal_source_maxSizeDotState = max(thermal_source_maxSizeDotState, &
maxval(thermalState(ph)%p%sizeDotState)) maxval(thermalState(ph)%p%sizeDotState))
enddo end do
end subroutine thermal_init end subroutine thermal_init
@ -170,7 +170,7 @@ module function phase_f_T(ph,en) result(f)
end select end select
enddo end do
end function phase_f_T end function phase_f_T
@ -195,7 +195,7 @@ function phase_thermal_collectDotState(ph,en) result(broken)
broken = broken .or. any(IEEE_is_NaN(thermalState(ph)%p(i)%dotState(:,en))) broken = broken .or. any(IEEE_is_NaN(thermalState(ph)%p(i)%dotState(:,en)))
enddo SourceLoop end do SourceLoop
end function phase_thermal_collectDotState end function phase_thermal_collectDotState
@ -262,7 +262,7 @@ function integrateThermalState(Delta_t, ph,en) result(broken)
sizeDotState = thermalState(ph)%p(so)%sizeDotState sizeDotState = thermalState(ph)%p(so)%sizeDotState
thermalState(ph)%p(so)%state(1:sizeDotState,en) = thermalState(ph)%p(so)%state0(1:sizeDotState,en) & thermalState(ph)%p(so)%state(1:sizeDotState,en) = thermalState(ph)%p(so)%state0(1:sizeDotState,en) &
+ thermalState(ph)%p(so)%dotState(1:sizeDotState,en) * Delta_t + thermalState(ph)%p(so)%dotState(1:sizeDotState,en) * Delta_t
enddo end do
end function integrateThermalState end function integrateThermalState
@ -277,7 +277,7 @@ module subroutine thermal_restartWrite(groupHandle,ph)
do so = 1,thermal_Nsources(ph) do so = 1,thermal_Nsources(ph)
call HDF5_write(thermalState(ph)%p(so)%state,groupHandle,'omega_thermal') call HDF5_write(thermalState(ph)%p(so)%state,groupHandle,'omega_thermal')
enddo end do
end subroutine thermal_restartWrite end subroutine thermal_restartWrite
@ -292,7 +292,7 @@ module subroutine thermal_restartRead(groupHandle,ph)
do so = 1,thermal_Nsources(ph) do so = 1,thermal_Nsources(ph)
call HDF5_read(thermalState(ph)%p(so)%state0,groupHandle,'omega_thermal') call HDF5_read(thermalState(ph)%p(so)%state0,groupHandle,'omega_thermal')
enddo end do
end subroutine thermal_restartRead end subroutine thermal_restartRead
@ -305,8 +305,8 @@ module subroutine thermal_forward()
do ph = 1, size(thermalState) do ph = 1, size(thermalState)
do so = 1, size(thermalState(ph)%p) do so = 1, size(thermalState(ph)%p)
thermalState(ph)%p(so)%state0 = thermalState(ph)%p(so)%state thermalState(ph)%p(so)%state0 = thermalState(ph)%p(so)%state
enddo end do
enddo end do
end subroutine thermal_forward end subroutine thermal_forward
@ -380,8 +380,8 @@ function thermal_active(source_label,src_length) result(active_source)
do s = 1, sources%length do s = 1, sources%length
src => sources%get(s) src => sources%get(s)
active_source(s,p) = src%get_asString('type') == source_label active_source(s,p) = src%get_asString('type') == source_label
enddo end do
enddo end do
end function thermal_active end function thermal_active

View File

@ -112,7 +112,7 @@ pure function eval(self,x) result(y)
#else #else
y = IEEE_FMA(y,x-self%x_ref,self%coef(o)) y = IEEE_FMA(y,x-self%x_ref,self%coef(o))
#endif #endif
enddo end do
end function eval end function eval

View File

@ -95,11 +95,11 @@ subroutine results_init(restart)
call results_openJobFile call results_openJobFile
call get_command(commandLine) call get_command(commandLine)
call results_addAttribute('call (restart at '//date//')',trim(commandLine)) call results_addAttribute('call (restart at '//date//')',trim(commandLine))
call h5gmove_f(resultsFile,'setup','tmp',hdferr) call H5Gmove_f(resultsFile,'setup','tmp',hdferr)
call results_addAttribute('description','input data used to run the simulation up to restart at '//date,'tmp') call results_addAttribute('description','input data used to run the simulation up to restart at '//date,'tmp')
call results_closeGroup(results_addGroup('setup')) call results_closeGroup(results_addGroup('setup'))
call results_addAttribute('description','input data used to run the simulation','setup') call results_addAttribute('description','input data used to run the simulation','setup')
call h5gmove_f(resultsFile,'tmp','setup/previous',hdferr) call H5Gmove_f(resultsFile,'tmp','setup/previous',hdferr)
end if end if
call results_closeJobFile call results_closeJobFile
@ -333,8 +333,8 @@ subroutine results_removeLink(link)
integer :: hdferr integer :: hdferr
call h5ldelete_f(resultsFile,link, hdferr) call H5Ldelete_f(resultsFile,link, hdferr)
if (hdferr < 0) call IO_error(1,ext_msg = 'results_removeLink: h5ldelete_soft_f ('//trim(link)//')') if (hdferr < 0) call IO_error(1,ext_msg = 'results_removeLink: H5Ldelete_soft_f ('//trim(link)//')')
end subroutine results_removeLink end subroutine results_removeLink
@ -522,7 +522,7 @@ subroutine results_mapping_phase(ID,entry,label)
writeSize = 0 writeSize = 0
writeSize(worldrank) = size(entry(1,:)) ! total number of entries of this process writeSize(worldrank) = size(entry(1,:)) ! total number of entries of this process
call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, hdferr) call H5Pcreate_f(H5P_DATASET_XFER_F, plist_id, hdferr)
if(hdferr < 0) error stop 'HDF5 error' if(hdferr < 0) error stop 'HDF5 error'
#ifndef PETSC #ifndef PETSC
@ -530,7 +530,7 @@ subroutine results_mapping_phase(ID,entry,label)
#else #else
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! MPI settings and communication ! MPI settings and communication
call h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, hdferr) call H5Pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, hdferr)
if(hdferr < 0) error stop 'HDF5 error' if(hdferr < 0) error stop 'HDF5 error'
call MPI_Allreduce(MPI_IN_PLACE,writeSize,worldsize,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,err_MPI) ! get output at each process call MPI_Allreduce(MPI_IN_PLACE,writeSize,worldsize,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,err_MPI) ! get output at each process
@ -558,82 +558,82 @@ subroutine results_mapping_phase(ID,entry,label)
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
! compound type: label(ID) + entry ! compound type: label(ID) + entry
call h5tcopy_f(H5T_NATIVE_CHARACTER, dt_id, hdferr) call H5Tcopy_f(H5T_NATIVE_CHARACTER, dt_id, hdferr)
if(hdferr < 0) error stop 'HDF5 error' if(hdferr < 0) error stop 'HDF5 error'
call h5tset_size_f(dt_id, int(len(label(1)),SIZE_T), hdferr) call H5Tset_size_f(dt_id, int(len(label(1)),SIZE_T), hdferr)
if(hdferr < 0) error stop 'HDF5 error' if(hdferr < 0) error stop 'HDF5 error'
call h5tget_size_f(dt_id, type_size_string, hdferr) call H5Tget_size_f(dt_id, type_size_string, hdferr)
if(hdferr < 0) error stop 'HDF5 error' if(hdferr < 0) error stop 'HDF5 error'
pI64_t = h5kind_to_type(kind(entryGlobal),H5_INTEGER_KIND) pI64_t = h5kind_to_type(kind(entryGlobal),H5_INTEGER_KIND)
call h5tget_size_f(pI64_t, type_size_int, hdferr) call H5Tget_size_f(pI64_t, type_size_int, hdferr)
if(hdferr < 0) error stop 'HDF5 error' if(hdferr < 0) error stop 'HDF5 error'
call h5tcreate_f(H5T_COMPOUND_F, type_size_string + type_size_int, dtype_id, hdferr) call H5Tcreate_f(H5T_COMPOUND_F, type_size_string + type_size_int, dtype_id, hdferr)
if(hdferr < 0) error stop 'HDF5 error' if(hdferr < 0) error stop 'HDF5 error'
call h5tinsert_f(dtype_id, 'label', 0_SIZE_T, dt_id,hdferr) call H5Tinsert_f(dtype_id, 'label', 0_SIZE_T, dt_id,hdferr)
if(hdferr < 0) error stop 'HDF5 error' if(hdferr < 0) error stop 'HDF5 error'
call h5tinsert_f(dtype_id, 'entry', type_size_string, pI64_t, hdferr) call H5Tinsert_f(dtype_id, 'entry', type_size_string, pI64_t, hdferr)
if(hdferr < 0) error stop 'HDF5 error' if(hdferr < 0) error stop 'HDF5 error'
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! create memory types for each component of the compound type ! create memory types for each component of the compound type
call h5tcreate_f(H5T_COMPOUND_F, type_size_string, label_id, hdferr) call H5Tcreate_f(H5T_COMPOUND_F, type_size_string, label_id, hdferr)
if(hdferr < 0) error stop 'HDF5 error' if(hdferr < 0) error stop 'HDF5 error'
call h5tinsert_f(label_id, 'label', 0_SIZE_T, dt_id, hdferr) call H5Tinsert_f(label_id, 'label', 0_SIZE_T, dt_id, hdferr)
if(hdferr < 0) error stop 'HDF5 error' if(hdferr < 0) error stop 'HDF5 error'
call h5tcreate_f(H5T_COMPOUND_F, type_size_int, entry_id, hdferr) call H5Tcreate_f(H5T_COMPOUND_F, type_size_int, entry_id, hdferr)
if(hdferr < 0) error stop 'HDF5 error' if(hdferr < 0) error stop 'HDF5 error'
call h5tinsert_f(entry_id, 'entry', 0_SIZE_T, pI64_t, hdferr) call H5Tinsert_f(entry_id, 'entry', 0_SIZE_T, pI64_t, hdferr)
if(hdferr < 0) error stop 'HDF5 error' if(hdferr < 0) error stop 'HDF5 error'
call h5tclose_f(dt_id, hdferr) call H5Tclose_f(dt_id, hdferr)
if(hdferr < 0) error stop 'HDF5 error' if(hdferr < 0) error stop 'HDF5 error'
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! create dataspace in memory (local shape = hyperslab) and in file (global shape) ! create dataspace in memory (local shape = hyperslab) and in file (global shape)
call h5screate_simple_f(2,myShape,memspace_id,hdferr,myShape) call H5Screate_simple_f(2,myShape,memspace_id,hdferr,myShape)
if(hdferr < 0) error stop 'HDF5 error' if(hdferr < 0) error stop 'HDF5 error'
call h5screate_simple_f(2,totalShape,filespace_id,hdferr,totalShape) call H5Screate_simple_f(2,totalShape,filespace_id,hdferr,totalShape)
if(hdferr < 0) error stop 'HDF5 error' if(hdferr < 0) error stop 'HDF5 error'
call h5sselect_hyperslab_f(filespace_id, H5S_SELECT_SET_F, myOffset, myShape, hdferr) call H5Sselect_hyperslab_f(filespace_id, H5S_SELECT_SET_F, myOffset, myShape, hdferr)
if(hdferr < 0) error stop 'HDF5 error' if(hdferr < 0) error stop 'HDF5 error'
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! write the components of the compound type individually ! write the components of the compound type individually
call h5pset_preserve_f(plist_id, .true., hdferr) call H5Pset_preserve_f(plist_id, .true., hdferr)
if(hdferr < 0) error stop 'HDF5 error' if(hdferr < 0) error stop 'HDF5 error'
loc_id = results_openGroup('/cell_to') loc_id = results_openGroup('/cell_to')
call h5dcreate_f(loc_id, 'phase', dtype_id, filespace_id, dset_id, hdferr) call H5Dcreate_f(loc_id, 'phase', dtype_id, filespace_id, dset_id, hdferr)
if(hdferr < 0) error stop 'HDF5 error' if(hdferr < 0) error stop 'HDF5 error'
call h5dwrite_f(dset_id, label_id, reshape(label(pack(ID,.true.)),myShape), & call H5Dwrite_f(dset_id, label_id, reshape(label(pack(ID,.true.)),myShape), &
myShape, hdferr, file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id) myShape, hdferr, file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id)
if(hdferr < 0) error stop 'HDF5 error' if(hdferr < 0) error stop 'HDF5 error'
call h5dwrite_f(dset_id, entry_id, reshape(pack(entryGlobal,.true.),myShape), & call H5Dwrite_f(dset_id, entry_id, reshape(pack(entryGlobal,.true.),myShape), &
myShape, hdferr, file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id) myShape, hdferr, file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id)
if(hdferr < 0) error stop 'HDF5 error' if(hdferr < 0) error stop 'HDF5 error'
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! close all ! close all
call HDF5_closeGroup(loc_id) call HDF5_closeGroup(loc_id)
call h5pclose_f(plist_id, hdferr) call H5Pclose_f(plist_id, hdferr)
if(hdferr < 0) error stop 'HDF5 error' if(hdferr < 0) error stop 'HDF5 error'
call h5sclose_f(filespace_id, hdferr) call H5Sclose_f(filespace_id, hdferr)
if(hdferr < 0) error stop 'HDF5 error' if(hdferr < 0) error stop 'HDF5 error'
call h5sclose_f(memspace_id, hdferr) call H5Sclose_f(memspace_id, hdferr)
if(hdferr < 0) error stop 'HDF5 error' if(hdferr < 0) error stop 'HDF5 error'
call h5dclose_f(dset_id, hdferr) call H5Dclose_f(dset_id, hdferr)
if(hdferr < 0) error stop 'HDF5 error' if(hdferr < 0) error stop 'HDF5 error'
call h5tclose_f(dtype_id, hdferr) call H5Tclose_f(dtype_id, hdferr)
if(hdferr < 0) error stop 'HDF5 error' if(hdferr < 0) error stop 'HDF5 error'
call h5tclose_f(label_id, hdferr) call H5Tclose_f(label_id, hdferr)
if(hdferr < 0) error stop 'HDF5 error' if(hdferr < 0) error stop 'HDF5 error'
call h5tclose_f(entry_id, hdferr) call H5Tclose_f(entry_id, hdferr)
call executionStamp('cell_to/phase','cell ID and constituent ID to phase results') call executionStamp('cell_to/phase','cell ID and constituent ID to phase results')
@ -678,7 +678,7 @@ subroutine results_mapping_homogenization(ID,entry,label)
writeSize = 0 writeSize = 0
writeSize(worldrank) = size(entry) ! total number of entries of this process writeSize(worldrank) = size(entry) ! total number of entries of this process
call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, hdferr) call H5Pcreate_f(H5P_DATASET_XFER_F, plist_id, hdferr)
if(hdferr < 0) error stop 'HDF5 error' if(hdferr < 0) error stop 'HDF5 error'
#ifndef PETSC #ifndef PETSC
@ -686,7 +686,7 @@ subroutine results_mapping_homogenization(ID,entry,label)
#else #else
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! MPI settings and communication ! MPI settings and communication
call h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, hdferr) call H5Pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, hdferr)
if(hdferr < 0) error stop 'HDF5 error' if(hdferr < 0) error stop 'HDF5 error'
call MPI_Allreduce(MPI_IN_PLACE,writeSize,worldsize,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,err_MPI) ! get output at each process call MPI_Allreduce(MPI_IN_PLACE,writeSize,worldsize,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,err_MPI) ! get output at each process
@ -710,82 +710,82 @@ subroutine results_mapping_homogenization(ID,entry,label)
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
! compound type: label(ID) + entry ! compound type: label(ID) + entry
call h5tcopy_f(H5T_NATIVE_CHARACTER, dt_id, hdferr) call H5Tcopy_f(H5T_NATIVE_CHARACTER, dt_id, hdferr)
if(hdferr < 0) error stop 'HDF5 error' if(hdferr < 0) error stop 'HDF5 error'
call h5tset_size_f(dt_id, int(len(label(1)),SIZE_T), hdferr) call H5Tset_size_f(dt_id, int(len(label(1)),SIZE_T), hdferr)
if(hdferr < 0) error stop 'HDF5 error' if(hdferr < 0) error stop 'HDF5 error'
call h5tget_size_f(dt_id, type_size_string, hdferr) call H5Tget_size_f(dt_id, type_size_string, hdferr)
if(hdferr < 0) error stop 'HDF5 error' if(hdferr < 0) error stop 'HDF5 error'
pI64_t = h5kind_to_type(kind(entryGlobal),H5_INTEGER_KIND) pI64_t = h5kind_to_type(kind(entryGlobal),H5_INTEGER_KIND)
call h5tget_size_f(pI64_t, type_size_int, hdferr) call H5Tget_size_f(pI64_t, type_size_int, hdferr)
if(hdferr < 0) error stop 'HDF5 error' if(hdferr < 0) error stop 'HDF5 error'
call h5tcreate_f(H5T_COMPOUND_F, type_size_string + type_size_int, dtype_id, hdferr) call H5Tcreate_f(H5T_COMPOUND_F, type_size_string + type_size_int, dtype_id, hdferr)
if(hdferr < 0) error stop 'HDF5 error' if(hdferr < 0) error stop 'HDF5 error'
call h5tinsert_f(dtype_id, 'label', 0_SIZE_T, dt_id,hdferr) call H5Tinsert_f(dtype_id, 'label', 0_SIZE_T, dt_id,hdferr)
if(hdferr < 0) error stop 'HDF5 error' if(hdferr < 0) error stop 'HDF5 error'
call h5tinsert_f(dtype_id, 'entry', type_size_string, pI64_t, hdferr) call H5Tinsert_f(dtype_id, 'entry', type_size_string, pI64_t, hdferr)
if(hdferr < 0) error stop 'HDF5 error' if(hdferr < 0) error stop 'HDF5 error'
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! create memory types for each component of the compound type ! create memory types for each component of the compound type
call h5tcreate_f(H5T_COMPOUND_F, type_size_string, label_id, hdferr) call H5Tcreate_f(H5T_COMPOUND_F, type_size_string, label_id, hdferr)
if(hdferr < 0) error stop 'HDF5 error' if(hdferr < 0) error stop 'HDF5 error'
call h5tinsert_f(label_id, 'label', 0_SIZE_T, dt_id, hdferr) call H5Tinsert_f(label_id, 'label', 0_SIZE_T, dt_id, hdferr)
if(hdferr < 0) error stop 'HDF5 error' if(hdferr < 0) error stop 'HDF5 error'
call h5tcreate_f(H5T_COMPOUND_F, type_size_int, entry_id, hdferr) call H5Tcreate_f(H5T_COMPOUND_F, type_size_int, entry_id, hdferr)
if(hdferr < 0) error stop 'HDF5 error' if(hdferr < 0) error stop 'HDF5 error'
call h5tinsert_f(entry_id, 'entry', 0_SIZE_T, pI64_t, hdferr) call H5Tinsert_f(entry_id, 'entry', 0_SIZE_T, pI64_t, hdferr)
if(hdferr < 0) error stop 'HDF5 error' if(hdferr < 0) error stop 'HDF5 error'
call h5tclose_f(dt_id, hdferr) call H5Tclose_f(dt_id, hdferr)
if(hdferr < 0) error stop 'HDF5 error' if(hdferr < 0) error stop 'HDF5 error'
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! create dataspace in memory (local shape = hyperslab) and in file (global shape) ! create dataspace in memory (local shape = hyperslab) and in file (global shape)
call h5screate_simple_f(1,myShape,memspace_id,hdferr,myShape) call H5Screate_simple_f(1,myShape,memspace_id,hdferr,myShape)
if(hdferr < 0) error stop 'HDF5 error' if(hdferr < 0) error stop 'HDF5 error'
call h5screate_simple_f(1,totalShape,filespace_id,hdferr,totalShape) call H5Screate_simple_f(1,totalShape,filespace_id,hdferr,totalShape)
if(hdferr < 0) error stop 'HDF5 error' if(hdferr < 0) error stop 'HDF5 error'
call h5sselect_hyperslab_f(filespace_id, H5S_SELECT_SET_F, myOffset, myShape, hdferr) call H5Sselect_hyperslab_f(filespace_id, H5S_SELECT_SET_F, myOffset, myShape, hdferr)
if(hdferr < 0) error stop 'HDF5 error' if(hdferr < 0) error stop 'HDF5 error'
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! write the components of the compound type individually ! write the components of the compound type individually
call h5pset_preserve_f(plist_id, .true., hdferr) call H5Pset_preserve_f(plist_id, .true., hdferr)
if(hdferr < 0) error stop 'HDF5 error' if(hdferr < 0) error stop 'HDF5 error'
loc_id = results_openGroup('/cell_to') loc_id = results_openGroup('/cell_to')
call h5dcreate_f(loc_id, 'homogenization', dtype_id, filespace_id, dset_id, hdferr) call H5Dcreate_f(loc_id, 'homogenization', dtype_id, filespace_id, dset_id, hdferr)
if(hdferr < 0) error stop 'HDF5 error' if(hdferr < 0) error stop 'HDF5 error'
call h5dwrite_f(dset_id, label_id, reshape(label(pack(ID,.true.)),myShape), & call H5Dwrite_f(dset_id, label_id, reshape(label(pack(ID,.true.)),myShape), &
myShape, hdferr, file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id) myShape, hdferr, file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id)
if(hdferr < 0) error stop 'HDF5 error' if(hdferr < 0) error stop 'HDF5 error'
call h5dwrite_f(dset_id, entry_id, reshape(pack(entryGlobal,.true.),myShape), & call H5Dwrite_f(dset_id, entry_id, reshape(pack(entryGlobal,.true.),myShape), &
myShape, hdferr, file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id) myShape, hdferr, file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id)
if(hdferr < 0) error stop 'HDF5 error' if(hdferr < 0) error stop 'HDF5 error'
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! close all ! close all
call HDF5_closeGroup(loc_id) call HDF5_closeGroup(loc_id)
call h5pclose_f(plist_id, hdferr) call H5Pclose_f(plist_id, hdferr)
if(hdferr < 0) error stop 'HDF5 error' if(hdferr < 0) error stop 'HDF5 error'
call h5sclose_f(filespace_id, hdferr) call H5Sclose_f(filespace_id, hdferr)
if(hdferr < 0) error stop 'HDF5 error' if(hdferr < 0) error stop 'HDF5 error'
call h5sclose_f(memspace_id, hdferr) call H5Sclose_f(memspace_id, hdferr)
if(hdferr < 0) error stop 'HDF5 error' if(hdferr < 0) error stop 'HDF5 error'
call h5dclose_f(dset_id, hdferr) call H5Dclose_f(dset_id, hdferr)
if(hdferr < 0) error stop 'HDF5 error' if(hdferr < 0) error stop 'HDF5 error'
call h5tclose_f(dtype_id, hdferr) call H5Tclose_f(dtype_id, hdferr)
if(hdferr < 0) error stop 'HDF5 error' if(hdferr < 0) error stop 'HDF5 error'
call h5tclose_f(label_id, hdferr) call H5Tclose_f(label_id, hdferr)
if(hdferr < 0) error stop 'HDF5 error' if(hdferr < 0) error stop 'HDF5 error'
call h5tclose_f(entry_id, hdferr) call H5Tclose_f(entry_id, hdferr)
if(hdferr < 0) error stop 'HDF5 error' if(hdferr < 0) error stop 'HDF5 error'
call executionStamp('cell_to/homogenization','cell ID to homogenization results') call executionStamp('cell_to/homogenization','cell ID to homogenization results')

View File

@ -182,7 +182,7 @@ subroutine fromEulers(self,eu,degrees)
Eulers = eu Eulers = eu
else else
Eulers = merge(eu*INRAD,eu,degrees) Eulers = merge(eu*INRAD,eu,degrees)
endif end if
if (any(Eulers<0.0_pReal) .or. any(Eulers>TAU) .or. Eulers(2) > PI) & if (any(Eulers<0.0_pReal) .or. any(Eulers>TAU) .or. Eulers(2) > PI) &
call IO_error(402,ext_msg='fromEulers') call IO_error(402,ext_msg='fromEulers')
@ -206,14 +206,14 @@ subroutine fromAxisAngle(self,ax,degrees,P)
angle = ax(4) angle = ax(4)
else else
angle = merge(ax(4)*INRAD,ax(4),degrees) angle = merge(ax(4)*INRAD,ax(4),degrees)
endif end if
if (.not. present(P)) then if (.not. present(P)) then
axis = ax(1:3) axis = ax(1:3)
else else
axis = ax(1:3) * merge(-1.0_pReal,1.0_pReal,P == 1) axis = ax(1:3) * merge(-1.0_pReal,1.0_pReal,P == 1)
if(abs(P) /= 1) call IO_error(402,ext_msg='fromAxisAngle (P)') if(abs(P) /= 1) call IO_error(402,ext_msg='fromAxisAngle (P)')
endif end if
if(dNeq(norm2(axis),1.0_pReal) .or. angle < 0.0_pReal .or. angle > PI) & if(dNeq(norm2(axis),1.0_pReal) .or. angle < 0.0_pReal .or. angle > PI) &
call IO_error(402,ext_msg='fromAxisAngle') call IO_error(402,ext_msg='fromAxisAngle')
@ -284,7 +284,7 @@ pure function rotVector(self,v,active) result(vRot)
passive = .not. active passive = .not. active
else else
passive = .true. passive = .true.
endif end if
if (dEq0(norm2(v))) then if (dEq0(norm2(v))) then
vRot = v vRot = v
@ -294,7 +294,7 @@ pure function rotVector(self,v,active) result(vRot)
multiplyQuaternion(conjugateQuaternion(self%q), multiplyQuaternion(v_normed, self%q)), & multiplyQuaternion(conjugateQuaternion(self%q), multiplyQuaternion(v_normed, self%q)), &
passive) passive)
vRot = q(2:4)*norm2(v) vRot = q(2:4)*norm2(v)
endif end if
end function rotVector end function rotVector
@ -318,7 +318,7 @@ pure function rotTensor2(self,T,active) result(tRot)
passive = .not. active passive = .not. active
else else
passive = .true. passive = .true.
endif end if
tRot = merge(matmul(matmul(self%asMatrix(),T),transpose(self%asMatrix())), & tRot = merge(matmul(matmul(self%asMatrix(),T),transpose(self%asMatrix())), &
matmul(matmul(transpose(self%asMatrix()),T),self%asMatrix()), & matmul(matmul(transpose(self%asMatrix()),T),self%asMatrix()), &
@ -347,14 +347,14 @@ pure function rotTensor4(self,T,active) result(tRot)
R = merge(transpose(self%asMatrix()),self%asMatrix(),active) R = merge(transpose(self%asMatrix()),self%asMatrix(),active)
else else
R = self%asMatrix() R = self%asMatrix()
endif end if
tRot = 0.0_pReal tRot = 0.0_pReal
do i = 1,3;do j = 1,3;do k = 1,3;do l = 1,3 do i = 1,3;do j = 1,3;do k = 1,3;do l = 1,3
do m = 1,3;do n = 1,3;do o = 1,3;do p = 1,3 do m = 1,3;do n = 1,3;do o = 1,3;do p = 1,3
tRot(i,j,k,l) = tRot(i,j,k,l) & tRot(i,j,k,l) = tRot(i,j,k,l) &
+ R(i,m) * R(j,n) * R(k,o) * R(l,p) * T(m,n,o,p) + R(i,m) * R(j,n) * R(k,o) * R(l,p) * T(m,n,o,p)
enddo; enddo; enddo; enddo; enddo; enddo; enddo; enddo end do; end do; end do; end do; end do; end do; end do; end do
end function rotTensor4 end function rotTensor4
@ -379,7 +379,7 @@ pure function rotStiffness(self,C,active) result(cRot)
R = merge(transpose(self%asMatrix()),self%asMatrix(),active) R = merge(transpose(self%asMatrix()),self%asMatrix(),active)
else else
R = self%asMatrix() R = self%asMatrix()
endif end if
M = reshape([R(1,1)**2, R(2,1)**2, R(3,1)**2, & M = reshape([R(1,1)**2, R(2,1)**2, R(3,1)**2, &
R(2,1)*R(3,1), R(1,1)*R(3,1), R(1,1)*R(2,1), & R(2,1)*R(3,1), R(1,1)*R(3,1), R(1,1)*R(2,1), &
@ -468,7 +468,7 @@ pure function qu2eu(qu) result(eu)
eu = [atan2((-P*qu(1)*qu(3)+qu(2)*qu(4))*chi, (-P*qu(1)*qu(2)-qu(3)*qu(4))*chi ), & eu = [atan2((-P*qu(1)*qu(3)+qu(2)*qu(4))*chi, (-P*qu(1)*qu(2)-qu(3)*qu(4))*chi ), &
atan2( 2.0_pReal*chi, q03-q12 ), & atan2( 2.0_pReal*chi, q03-q12 ), &
atan2(( P*qu(1)*qu(3)+qu(2)*qu(4))*chi, (-P*qu(1)*qu(2)+qu(3)*qu(4))*chi )] atan2(( P*qu(1)*qu(3)+qu(2)*qu(4))*chi, (-P*qu(1)*qu(2)+qu(3)*qu(4))*chi )]
endif degenerated end if degenerated
where(sign(1.0_pReal,eu)<0.0_pReal) eu = mod(eu+TAU,[TAU,PI,TAU]) where(sign(1.0_pReal,eu)<0.0_pReal) eu = mod(eu+TAU,[TAU,PI,TAU])
end function qu2eu end function qu2eu
@ -526,8 +526,8 @@ pure function om2qu(om) result(qu)
else else
s = 2.0_pReal * sqrt( 1.0_pReal + om(3,3) - om(1,1) - om(2,2) ) s = 2.0_pReal * sqrt( 1.0_pReal + om(3,3) - om(1,1) - om(2,2) )
qu = [ (om(2,1) - om(1,2)) /s,(om(1,3) + om(3,1)) / s,(om(2,3) + om(3,2)) / s,0.25_pReal * s] qu = [ (om(2,1) - om(1,2)) /s,(om(1,3) + om(3,1)) / s,(om(2,3) + om(3,2)) / s,0.25_pReal * s]
endif end if
endif end if
if(sign(1.0_pReal,qu(1))<0.0_pReal) qu =-1.0_pReal * qu if(sign(1.0_pReal,qu(1))<0.0_pReal) qu =-1.0_pReal * qu
qu(2:4) = merge(qu(2:4),qu(2:4)*P,dEq0(qu(2:4))) qu(2:4) = merge(qu(2:4),qu(2:4)*P,dEq0(qu(2:4)))
qu = qu/norm2(qu) qu = qu/norm2(qu)
@ -593,7 +593,7 @@ function om2ax(om) result(ax)
ax(1:3) = VR(1:3,i) ax(1:3) = VR(1:3,i)
where ( dNeq0([om(2,3)-om(3,2), om(3,1)-om(1,3), om(1,2)-om(2,1)])) & where ( dNeq0([om(2,3)-om(3,2), om(3,1)-om(1,3), om(1,2)-om(2,1)])) &
ax(1:3) = sign(ax(1:3),-P *[om(2,3)-om(3,2), om(3,1)-om(1,3), om(1,2)-om(2,1)]) ax(1:3) = sign(ax(1:3),-P *[om(2,3)-om(3,2), om(3,1)-om(1,3), om(1,2)-om(2,1)])
endif end if
end function om2ax end function om2ax
@ -826,7 +826,7 @@ subroutine selfTest()
cos(TAU*x(2))*B,& cos(TAU*x(2))*B,&
sin(TAU*x(1))*A] sin(TAU*x(1))*A]
if(qu(1)<0.0_pReal) qu = qu * (-1.0_pReal) if(qu(1)<0.0_pReal) qu = qu * (-1.0_pReal)
endif end if
if(.not. quaternion_equal(om2qu(qu2om(qu)),qu)) error stop 'om2qu2om' if(.not. quaternion_equal(om2qu(qu2om(qu)),qu)) error stop 'om2qu2om'

View File

@ -123,7 +123,7 @@ function getCWD()
getCWD = c_f_string(getCWD_Cstring) getCWD = c_f_string(getCWD_Cstring)
else else
error stop 'invalid working directory' error stop 'invalid working directory'
endif end if
end function getCWD end function getCWD
@ -145,7 +145,7 @@ function getHostName()
getHostName = c_f_string(getHostName_Cstring) getHostName = c_f_string(getHostName_Cstring)
else else
getHostName = 'n/a (Error!)' getHostName = 'n/a (Error!)'
endif end if
end function getHostName end function getHostName