223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

improved reading in of values, now only warnings in case of problematic entries in material.config 

divergence calculation sqrt scaling optionally introduced for basic scheme spectral solver
This commit is contained in:
Martin Diehl 2013-02-08 15:55:53 +00:00
parent 59e59c90c5
commit e644c6dbc5
5 changed files with 159 additions and 130 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

124
code/IO.f90
View File

@ -844,17 +844,25 @@ end function IO_stringPos
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief read string value at myPos from line !> @brief read string value at myPos from line
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function IO_stringValue(line,positions,myPos) function IO_stringValue(line,positions,myPos,silent)
implicit none implicit none
integer(pInt), dimension(:), intent(in) :: positions integer(pInt), dimension(:), intent(in) :: positions
integer(pInt), intent(in) :: myPos integer(pInt), intent(in) :: myPos
character(len=1+positions(myPos*2+1)-positions(myPos*2)) :: IO_stringValue character(len=1+positions(myPos*2+1)-positions(myPos*2)) :: IO_stringValue
character(len=*), intent(in) :: line character(len=*), intent(in) :: line
logical, optional,intent(in) :: silent
logical :: warn
if (myPos > positions(1)) then if (.not. present(silent)) then
IO_stringValue = '' warn = .false.
call IO_warning(201, e=myPos, ext_msg = trim(line)//' (IO_stringValue)') else
warn = silent
endif
IO_stringValue = ''
if (myPos > positions(1) .or. myPos < 1_pInt) then ! trying to access non-present value
if (warn) call IO_warning(201, e=myPos, ext_msg = 'IO_stringValue: '//trim(line))
else else
IO_stringValue = line(positions(myPos*2):positions(myPos*2+1)) IO_stringValue = line(positions(myPos*2):positions(myPos*2+1))
endif endif
@ -868,7 +876,6 @@ end function IO_stringValue
pure function IO_fixedStringValue (line,ends,myPos) pure function IO_fixedStringValue (line,ends,myPos)
implicit none implicit none
integer(pInt), intent(in) :: myPos integer(pInt), intent(in) :: myPos
integer(pInt), dimension(:), intent(in) :: ends integer(pInt), dimension(:), intent(in) :: ends
character(len=ends(myPos+1)-ends(myPos)) :: IO_fixedStringValue character(len=ends(myPos+1)-ends(myPos)) :: IO_fixedStringValue
@ -888,25 +895,25 @@ real(pReal) function IO_floatValue (line,positions,myPos)
character(len=*), intent(in) :: line character(len=*), intent(in) :: line
integer(pInt), dimension(:), intent(in) :: positions integer(pInt), dimension(:), intent(in) :: positions
integer(pInt), intent(in) :: myPos integer(pInt), intent(in) :: myPos
character(len=64), parameter :: myName = 'IO_floatValue' character(len=15), parameter :: myName = 'IO_floatValue: '
character(len=17), parameter :: validCharacters = '0123456789eEdD.+-' character(len=17), parameter :: validCharacters = '0123456789eEdD.+-'
integer(pInt) :: readStatus, invalidWhere integer(pInt) :: readStatus, invalidWhere
IO_floatValue = 0.0_pReal IO_floatValue = 0.0_pReal
if (myPos > positions(1) .or. myPos < 1_pInt) then ! trying to access non-present value
if (myPos > positions(1) .or. myPos < 1_pInt) then ! trying to access non-present value call IO_warning(201,ext_msg=myName//trim(line))
call IO_warning(201,ext_msg=trim(line)//' ('//trim(myName)//')')
else else
invalidWhere = verify(line(positions(myPos*2):positions(myPos*2+1)),validCharacters) invalidWhere = verify(line(positions(myPos*2):positions(myPos*2+1)),validCharacters) ! search for invalid characters
if (invalidWhere /= 0_pInt) then if (invalidWhere /= 0_pInt) then ! found invaldid character, only read in substring
invalidWhere = invalidWhere-1 invalidWhere = invalidWhere - 1_pInt
call IO_warning(202,ext_msg=line(positions(myPos*2):positions(myPos*2+1))//' ('//trim(myName)//')') call IO_warning(202,ext_msg=myName//line(positions(myPos*2):positions(myPos*2+1)))
else else
invalidWhere = positions(myPos*2+1)-positions(myPos*2)+1 invalidWhere = positions(myPos*2+1)-positions(myPos*2) + 1_pInt ! read until position(myPos*2+1)
endif endif
read(UNIT=line(positions(myPos*2):positions(myPos*2)+invalidWhere-1),iostat=readStatus,FMT=*) IO_floatValue read(UNIT=line(positions(myPos*2):positions(myPos*2)+invalidWhere-1),iostat=readStatus,FMT=*) &
if (readStatus /= 0_pInt) & IO_floatValue
call IO_warning(203,ext_msg=line(positions(myPos*2):positions(myPos*2)+invalidWhere-1)//' ('//trim(myName)//')') if (readStatus /= 0_pInt) & ! error during string to float conversion
call IO_warning(203,ext_msg=myName//line(positions(myPos*2):positions(myPos*2)+invalidWhere-1))
endif endif
end function IO_floatValue end function IO_floatValue
@ -921,24 +928,23 @@ real(pReal) function IO_fixedFloatValue (line,ends,myPos)
character(len=*), intent(in) :: line character(len=*), intent(in) :: line
integer(pInt), intent(in) :: myPos integer(pInt), intent(in) :: myPos
integer(pInt), dimension(:), intent(in) :: ends integer(pInt), dimension(:), intent(in) :: ends
character(len=64), parameter :: myName = 'IO_fixedFloatValue' character(len=20), parameter :: myName = 'IO_fixedFloatValue: '
character(len=17), parameter :: validCharacters = '0123456789eEdD.+-' character(len=17), parameter :: validCharacters = '0123456789eEdD.+-'
integer(pInt) :: readStatus, myStart, invalidWhere integer(pInt) :: readStatus, myStart, invalidWhere
IO_fixedFloatValue = 0.0_pReal IO_fixedFloatValue = 0.0_pReal
myStart = ends(myPos-1) + 1_pInt
myStart = ends(myPos-1)+1 invalidWhere = verify(line(myStart:ends(myPos)),validCharacters) ! search for invalid character
if (invalidWhere /= 0_pInt) then ! found invaldid character, only read in substring
invalidWhere = verify(line(myStart:ends(myPos)),validCharacters) invalidWhere = invalidWhere - 1_pInt
if (invalidWhere /= 0_pInt) then call IO_warning(202,ext_msg=myName//line(myStart:ends(myPos)))
invalidWhere = invalidWhere-1
call IO_warning(202,ext_msg=line(myStart:ends(myPos))//' ('//trim(myName)//')')
else else
invalidWhere = ends(myPos)-myStart+1 invalidWhere = ends(myPos)-myStart + 1_pInt ! read until ends(myPos)
endif endif
read(UNIT=line(myStart:myStart+invalidWhere-1),iostat=readStatus,FMT=*) IO_fixedFloatValue read(UNIT=line(myStart:myStart+invalidWhere-1),iostat=readStatus,FMT=*) IO_fixedFloatValue
if (readStatus /= 0_pInt) & if (readStatus /= 0_pInt) & ! error during string to float conversion
call IO_warning(203,ext_msg=line(myStart:myStart+invalidWhere-1)//' ('//trim(myName)//')') call IO_warning(203,ext_msg=myName//line(myStart:myStart+invalidWhere-1))
end function IO_fixedFloatValue end function IO_fixedFloatValue
@ -952,42 +958,45 @@ real(pReal) function IO_fixedNoEFloatValue (line,ends,myPos)
character(len=*), intent(in) :: line character(len=*), intent(in) :: line
integer(pInt), intent(in) :: myPos integer(pInt), intent(in) :: myPos
integer(pInt), dimension(:), intent(in) :: ends integer(pInt), dimension(:), intent(in) :: ends
character(len=64), parameter :: myName = 'IO_fixedNoEFloatValue' character(len=22), parameter :: myName = 'IO_fixedNoEFloatValue '
character(len=13), parameter :: validBase = '0123456789.+-' character(len=13), parameter :: validBase = '0123456789.+-'
character(len=12), parameter :: validExp = '0123456789+-' character(len=12), parameter :: validExp = '0123456789+-'
integer(pInt) :: expon = 0, myStart, readStatus integer(pInt) :: expon, myStart, readStatus
integer :: pos_exp, end_base, end_exp integer :: pos_exp, end_base, end_exp
real(pReal) :: base = 0.0_pReal real(pReal) :: base
myStart = ends(myPos-1)+1 base = 0.0_pReal
expon = 0_pInt
myStart = ends(myPos-1) + 1_pInt
pos_exp = scan(line(myStart:ends(myPos)),'+-',back=.true.) pos_exp = scan(line(myStart:ends(myPos)),'+-',back=.true.)
if (pos_exp <= 1_pInt) & ! no exponent but only base if (pos_exp <= 1_pInt) & ! no exponent but only base
pos_exp = ends(myPos)-myStart+1 pos_exp = ends(myPos)-myStart + 1_pInt
! --- figure out base --- ! --- figure out base ---
end_base = verify(line(myStart:myStart+pos_exp-1),validBase) end_base = verify(line(myStart:myStart+pos_exp-1),validBase) ! search for invalid character in base
if (end_base /= 0_pInt) then ! invalid base if (end_base /= 0_pInt) then ! found invaldid character, only read in substring
end_base = end_base-1 end_base = end_base-1
call IO_warning(202, ext_msg = line(myStart:myStart+pos_exp-1)//' ('//trim(myName)//':base)') call IO_warning(202, ext_msg = myName//'(base): '//line(myStart:myStart+pos_exp-1))
else else
end_base = pos_exp end_base = pos_exp ! read until begin of exponent
endif endif
read(UNIT=line(myStart:myStart+end_base-1),iostat=readStatus,FMT=*) base read(UNIT=line(myStart:myStart+end_base-1),iostat=readStatus,FMT=*) base
if (readStatus /= 0_pInt) & if (readStatus /= 0_pInt) & ! error during string to float conversion
call IO_warning(203, ext_msg = line(myStart:myStart+end_base-1)//' ('//trim(myName)//':base)') call IO_warning(203, ext_msg = myName//'(base): '//line(myStart:myStart+end_base-1))
! --- figure out exponent --- ! --- figure out exponent ---
end_exp = verify(line(myStart+pos_exp:ends(myPos)),validExp) end_exp = verify(line(myStart+pos_exp:ends(myPos)),validExp) ! search for invalid character in exponent
if (end_exp /= 0_pInt) then ! invalid exponent if (end_exp /= 0_pInt) then ! found invaldid character, only read in substring
end_exp = end_exp-1 end_exp = end_exp - 1_pInt
call IO_warning(202, ext_msg = line(myStart+pos_exp:ends(myPos))//' ('//trim(myName)//':exp)') call IO_warning(202, ext_msg = myName//'(exp): '//line(myStart+pos_exp:ends(myPos)))
else else
end_exp = mystart-ends(myPos)+1 end_exp = mystart-ends(myPos) + 1_pInt ! read until end of string
endif endif
read(UNIT=line(myStart+pos_exp:myStart+end_exp-1),iostat=readStatus,FMT=*) expon read(UNIT=line(myStart+pos_exp:myStart+end_exp-1),iostat=readStatus,FMT=*) expon
if (readStatus /= 0_pInt) & if (readStatus /= 0_pInt) & ! error during string to float conversion
call IO_warning(203, ext_msg = line(myStart+pos_exp:myStart+end_exp-1)//' ('//trim(myName)//':exp)') call IO_warning(203, ext_msg = myName//'(base): '//line(myStart+pos_exp:myStart+end_exp-1))
IO_fixedNoEFloatValue = base*10.0_pReal**expon IO_fixedNoEFloatValue = base*10.0_pReal**expon
@ -1003,25 +1012,26 @@ integer(pInt) function IO_intValue(line,positions,myPos)
character(len=*), intent(in) :: line character(len=*), intent(in) :: line
integer(pInt), dimension(:), intent(in) :: positions integer(pInt), dimension(:), intent(in) :: positions
integer(pInt), intent(in) :: myPos integer(pInt), intent(in) :: myPos
character(len=64), parameter :: myName = 'IO_intValue' character(len=13), parameter :: myName = 'IO_intValue: '
character(len=12), parameter :: validCharacters = '0123456789+-' character(len=12), parameter :: validCharacters = '0123456789+-'
integer(pInt) :: readStatus, invalidWhere integer(pInt) :: readStatus, invalidWhere
IO_intValue = 0_pInt IO_intValue = 0_pInt
if (myPos > positions(1) .or. myPos < 1_pInt) then ! trying to access non-present value if (myPos > positions(1) .or. myPos < 1_pInt) then ! trying to access non-present value
call IO_warning(201,ext_msg=trim(line)//' ('//trim(myName)//')') call IO_warning(201,ext_msg=myName//trim(line))
else else
invalidWhere = verify(line(positions(myPos*2):positions(myPos*2+1)),validCharacters) invalidWhere = verify(line(positions(myPos*2):positions(myPos*2+1)),validCharacters)
if (invalidWhere /= 0_pInt) then if (invalidWhere /= 0_pInt) then ! found invaldid character, only read in substring
invalidWhere = invalidWhere-1 invalidWhere = invalidWhere-1
call IO_warning(202,ext_msg=line(positions(myPos*2):positions(myPos*2+1))//' ('//trim(myName)//')') call IO_warning(202,ext_msg=line(positions(myPos*2):positions(myPos*2+1)))
else else
invalidWhere = positions(myPos*2+1)-positions(myPos*2)+1 invalidWhere = positions(myPos*2+1)-positions(myPos*2)+1
endif endif
read(UNIT=line(positions(myPos*2):positions(myPos*2)+invalidWhere-1),iostat=readStatus,FMT=*) IO_intValue read(UNIT=line(positions(myPos*2):positions(myPos*2)+invalidWhere-1),iostat=readStatus,FMT=*) &
if (readStatus /= 0_pInt) & IO_intValue
call IO_warning(203,ext_msg=line(positions(myPos*2):positions(myPos*2)+invalidWhere-1)//' ('//trim(myName)//')') if (readStatus /= 0_pInt) & ! error during string to int conversion
call IO_warning(203,ext_msg=myName//line(positions(myPos*2):positions(myPos*2)+invalidWhere-1))
endif endif
end function IO_intValue end function IO_intValue
@ -1036,7 +1046,7 @@ integer(pInt) function IO_fixedIntValue(line,ends,myPos)
character(len=*), intent(in) :: line character(len=*), intent(in) :: line
integer(pInt), intent(in) :: myPos integer(pInt), intent(in) :: myPos
integer(pInt), dimension(:), intent(in) :: ends integer(pInt), dimension(:), intent(in) :: ends
character(len=64), parameter :: myName = 'IO_fixedIntValue' character(len=18), parameter :: myName = 'IO_fixedIntValue: '
character(len=13), parameter :: validCharacters = '0123456789.+-' character(len=13), parameter :: validCharacters = '0123456789.+-'
integer(pInt) :: readStatus, myStart, invalidWhere integer(pInt) :: readStatus, myStart, invalidWhere
@ -1047,13 +1057,13 @@ integer(pInt) function IO_fixedIntValue(line,ends,myPos)
invalidWhere = verify(line(myStart:ends(myPos)),validCharacters) invalidWhere = verify(line(myStart:ends(myPos)),validCharacters)
if (invalidWhere /= 0_pInt) then if (invalidWhere /= 0_pInt) then
invalidWhere = invalidWhere-1 invalidWhere = invalidWhere-1
call IO_warning(202,ext_msg=line(myStart:ends(myPos))//' ('//trim(myName)//')') call IO_warning(202,ext_msg=myName//line(myStart:ends(myPos)))
else else
invalidWhere = ends(myPos)-myStart+1 invalidWhere = ends(myPos)-myStart+1
endif endif
read(UNIT=line(myStart:myStart+invalidWhere-1),iostat=readStatus,FMT=*) IO_fixedIntValue read(UNIT=line(myStart:myStart+invalidWhere-1),iostat=readStatus,FMT=*) IO_fixedIntValue
if (readStatus /= 0_pInt) & if (readStatus /= 0_pInt) &
call IO_warning(203,ext_msg=line(myStart:myStart+invalidWhere-1)//' ('//trim(myName)//')') call IO_warning(203,ext_msg=myName//line(myStart:myStart+invalidWhere-1))
end function IO_fixedIntValue end function IO_fixedIntValue

59
code/constitutive_phenopowerlaw.f90
View File

@ -125,20 +125,15 @@ subroutine constitutive_phenopowerlaw_init(myFile)
use debug, only: debug_level,& use debug, only: debug_level,&
debug_constitutive,& debug_constitutive,&
debug_levelBasic debug_levelBasic
use lattice, only: lattice_initializeStructure, lattice_symmetrizeC66, & use lattice
lattice_maxNslipFamily, lattice_maxNtwinFamily, &
lattice_maxNinteraction, lattice_NslipSystem, lattice_NtwinSystem, &
lattice_maxNonSchmid, &
lattice_interactionSlipSlip, &
lattice_interactionSlipTwin, &
lattice_interactionTwinSlip, &
lattice_interactionTwinTwin
implicit none implicit none
integer(pInt), intent(in) :: myFile integer(pInt), intent(in) :: myFile
integer(pInt), parameter :: maxNchunks = lattice_maxNinteraction + 1_pInt integer(pInt), parameter :: maxNchunks = lattice_maxNinteraction + 1_pInt
integer(pInt), dimension(1+2*maxNchunks) :: positions integer(pInt), dimension(1+2*maxNchunks) :: positions
integer(pInt) section, maxNinstance, i,j,k, f,o, & integer(pInt) section, maxNinstance, i,j,k, f,o, &
Nchunks_SlipSlip, Nchunks_SlipTwin, Nchunks_TwinSlip, Nchunks_TwinTwin, &
Nchunks_SlipFamilies, Nchunks_TwinFamilies, &
mySize, myStructure, index_myFamily, index_otherFamily mySize, myStructure, index_myFamily, index_otherFamily
character(len=64) :: tag character(len=64) :: tag
character(len=1024) :: line = '' ! to start initialized character(len=1024) :: line = '' ! to start initialized
@ -151,6 +146,13 @@ subroutine constitutive_phenopowerlaw_init(myFile)
maxNinstance = int(count(phase_plasticity == constitutive_phenopowerlaw_label),pInt) maxNinstance = int(count(phase_plasticity == constitutive_phenopowerlaw_label),pInt)
if (maxNinstance == 0) return if (maxNinstance == 0) return
Nchunks_SlipSlip = lattice_maxNinteraction
Nchunks_SlipTwin = lattice_maxNinteraction
Nchunks_TwinSlip = lattice_maxNinteraction
Nchunks_TwinTwin = lattice_maxNinteraction
Nchunks_SlipFamilies = lattice_maxNslipFamily
Nchunks_TwinFamilies = lattice_maxNtwinFamily
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) then if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) then
write(6,'(a16,1x,i5)') '# instances:',maxNinstance write(6,'(a16,1x,i5)') '# instances:',maxNinstance
write(6,*) write(6,*)
@ -262,6 +264,29 @@ subroutine constitutive_phenopowerlaw_init(myFile)
IO_lc(IO_stringValue(line,positions,2_pInt)) IO_lc(IO_stringValue(line,positions,2_pInt))
case ('lattice_structure') case ('lattice_structure')
constitutive_phenopowerlaw_structureName(i) = IO_lc(IO_stringValue(line,positions,2_pInt)) constitutive_phenopowerlaw_structureName(i) = IO_lc(IO_stringValue(line,positions,2_pInt))
select case (constitutive_phenopowerlaw_structureName(i))
case ('fcc')
Nchunks_SlipSlip = maxval(lattice_fcc_interactionSlipSlip)
Nchunks_SlipTwin = maxval(lattice_fcc_interactionSlipTwin)
Nchunks_TwinSlip = maxval(lattice_fcc_interactionTwinSlip)
Nchunks_TwinTwin = maxval(lattice_fcc_interactionTwinTwin)
Nchunks_SlipFamilies = count(lattice_fcc_NslipSystem > 0_pInt)
Nchunks_TwinFamilies = count(lattice_fcc_NtwinSystem > 0_pInt)
case ('bcc')
Nchunks_SlipSlip = maxval(lattice_bcc_interactionSlipSlip)
Nchunks_SlipTwin = maxval(lattice_bcc_interactionSlipTwin)
Nchunks_TwinSlip = maxval(lattice_bcc_interactionTwinSlip)
Nchunks_TwinTwin = maxval(lattice_bcc_interactionTwinTwin)
Nchunks_SlipFamilies = count(lattice_bcc_NslipSystem > 0_pInt)
Nchunks_TwinFamilies = count(lattice_bcc_NtwinSystem > 0_pInt)
case ('hex')
Nchunks_SlipSlip = maxval(lattice_hex_interactionSlipSlip)
Nchunks_SlipTwin = maxval(lattice_hex_interactionSlipTwin)
Nchunks_TwinSlip = maxval(lattice_hex_interactionTwinSlip)
Nchunks_TwinTwin = maxval(lattice_hex_interactionTwinTwin)
Nchunks_SlipFamilies = count(lattice_hex_NslipSystem > 0_pInt)
Nchunks_TwinFamilies = count(lattice_hex_NtwinSystem > 0_pInt)
end select
case ('covera_ratio') case ('covera_ratio')
constitutive_phenopowerlaw_CoverA(i) = IO_floatValue(line,positions,2_pInt) constitutive_phenopowerlaw_CoverA(i) = IO_floatValue(line,positions,2_pInt)
case ('c11') case ('c11')
@ -283,7 +308,7 @@ subroutine constitutive_phenopowerlaw_init(myFile)
case ('c66') case ('c66')
constitutive_phenopowerlaw_Cslip_66(6,6,i) = IO_floatValue(line,positions,2_pInt) constitutive_phenopowerlaw_Cslip_66(6,6,i) = IO_floatValue(line,positions,2_pInt)
case ('nslip') case ('nslip')
do j = 1_pInt, lattice_maxNslipFamily do j = 1_pInt, Nchunks_SlipFamilies
constitutive_phenopowerlaw_Nslip(j,i) = IO_intValue(line,positions,1_pInt+j) constitutive_phenopowerlaw_Nslip(j,i) = IO_intValue(line,positions,1_pInt+j)
enddo enddo
case ('gdot0_slip') case ('gdot0_slip')
@ -291,17 +316,17 @@ subroutine constitutive_phenopowerlaw_init(myFile)
case ('n_slip') case ('n_slip')
constitutive_phenopowerlaw_n_slip(i) = IO_floatValue(line,positions,2_pInt) constitutive_phenopowerlaw_n_slip(i) = IO_floatValue(line,positions,2_pInt)
case ('tau0_slip') case ('tau0_slip')
do j = 1_pInt, lattice_maxNslipFamily do j = 1_pInt, Nchunks_SlipFamilies
constitutive_phenopowerlaw_tau0_slip(j,i) = IO_floatValue(line,positions,1_pInt+j) constitutive_phenopowerlaw_tau0_slip(j,i) = IO_floatValue(line,positions,1_pInt+j)
enddo enddo
case ('tausat_slip') case ('tausat_slip')
do j = 1_pInt, lattice_maxNslipFamily do j = 1_pInt, Nchunks_SlipFamilies
constitutive_phenopowerlaw_tausat_slip(j,i) = IO_floatValue(line,positions,1_pInt+j) constitutive_phenopowerlaw_tausat_slip(j,i) = IO_floatValue(line,positions,1_pInt+j)
enddo enddo
case ('a_slip', 'w0_slip') case ('a_slip', 'w0_slip')
constitutive_phenopowerlaw_a_slip(i) = IO_floatValue(line,positions,2_pInt) constitutive_phenopowerlaw_a_slip(i) = IO_floatValue(line,positions,2_pInt)
case ('ntwin') case ('ntwin')
do j = 1_pInt, lattice_maxNtwinFamily do j = 1_pInt, Nchunks_TwinFamilies
constitutive_phenopowerlaw_Ntwin(j,i) = IO_intValue(line,positions,1_pInt+j) constitutive_phenopowerlaw_Ntwin(j,i) = IO_intValue(line,positions,1_pInt+j)
enddo enddo
case ('gdot0_twin') case ('gdot0_twin')
@ -309,7 +334,7 @@ subroutine constitutive_phenopowerlaw_init(myFile)
case ('n_twin') case ('n_twin')
constitutive_phenopowerlaw_n_twin(i) = IO_floatValue(line,positions,2_pInt) constitutive_phenopowerlaw_n_twin(i) = IO_floatValue(line,positions,2_pInt)
case ('tau0_twin') case ('tau0_twin')
do j = 1_pInt, lattice_maxNtwinFamily do j = 1_pInt, Nchunks_TwinFamilies
constitutive_phenopowerlaw_tau0_twin(j,i) = IO_floatValue(line,positions,1_pInt+j) constitutive_phenopowerlaw_tau0_twin(j,i) = IO_floatValue(line,positions,1_pInt+j)
enddo enddo
case ('s_pr') case ('s_pr')
@ -338,19 +363,19 @@ subroutine constitutive_phenopowerlaw_init(myFile)
case ('atol_twinfrac') case ('atol_twinfrac')
constitutive_phenopowerlaw_aTolTwinfrac(i) = IO_floatValue(line,positions,2_pInt) constitutive_phenopowerlaw_aTolTwinfrac(i) = IO_floatValue(line,positions,2_pInt)
case ('interaction_slipslip') case ('interaction_slipslip')
do j = 1_pInt, lattice_maxNinteraction do j = 1_pInt, Nchunks_SlipSlip
constitutive_phenopowerlaw_interaction_SlipSlip(j,i) = IO_floatValue(line,positions,1_pInt+j) constitutive_phenopowerlaw_interaction_SlipSlip(j,i) = IO_floatValue(line,positions,1_pInt+j)
enddo enddo
case ('interaction_sliptwin') case ('interaction_sliptwin')
do j = 1_pInt, lattice_maxNinteraction do j = 1_pInt, Nchunks_SlipTwin
constitutive_phenopowerlaw_interaction_SlipTwin(j,i) = IO_floatValue(line,positions,1_pInt+j) constitutive_phenopowerlaw_interaction_SlipTwin(j,i) = IO_floatValue(line,positions,1_pInt+j)
enddo enddo
case ('interaction_twinslip') case ('interaction_twinslip')
do j = 1_pInt, lattice_maxNinteraction do j = 1_pInt, Nchunks_TwinSlip
constitutive_phenopowerlaw_interaction_TwinSlip(j,i) = IO_floatValue(line,positions,1_pInt+j) constitutive_phenopowerlaw_interaction_TwinSlip(j,i) = IO_floatValue(line,positions,1_pInt+j)
enddo enddo
case ('interaction_twintwin') case ('interaction_twintwin')
do j = 1_pInt, lattice_maxNinteraction do j = 1_pInt, Nchunks_TwinTwin
constitutive_phenopowerlaw_interaction_TwinTwin(j,i) = IO_floatValue(line,positions,1_pInt+j) constitutive_phenopowerlaw_interaction_TwinTwin(j,i) = IO_floatValue(line,positions,1_pInt+j)
enddo enddo
case ('nonschmid_coefficients') case ('nonschmid_coefficients')

28
code/debug.f90
View File

@ -144,12 +144,9 @@ subroutine debug_init
integer(pInt), dimension(1+2*maxNchunks) :: positions integer(pInt), dimension(1+2*maxNchunks) :: positions
character(len=64) :: tag character(len=64) :: tag
character(len=1024) :: line character(len=1024) :: line
!$OMP CRITICAL (write2out) write(6,'(/,a)') ' <<<+- debug init -+>>>'
write(6,*) write(6,'(a)') ' $Id$'
write(6,*) '<<<+- debug init -+>>>'
write(6,*) '$Id$'
#include "compilation_info.f90" #include "compilation_info.f90"
!$OMP END CRITICAL (write2out)
if (allocated(debug_StressLoopDistribution)) & if (allocated(debug_StressLoopDistribution)) &
deallocate(debug_StressLoopDistribution) deallocate(debug_StressLoopDistribution)
@ -223,7 +220,7 @@ subroutine debug_init
what = debug_maxNtype + 2_pInt what = debug_maxNtype + 2_pInt
end select end select
if(what /= 0) then if(what /= 0) then
do i = 2_pInt, maxNchunks do i = 2_pInt, positions(1)
select case(IO_lc(IO_stringValue(line,positions,i))) select case(IO_lc(IO_stringValue(line,positions,i)))
case('basic') case('basic')
debug_level(what) = ior(debug_level(what), debug_levelBasic) debug_level(what) = ior(debug_level(what), debug_levelBasic)
@ -254,21 +251,14 @@ subroutine debug_init
debug_level(i) = ior(debug_level(i), debug_level(debug_maxNtype + 1_pInt)) ! fill all debug types with levels specified by "all" debug_level(i) = ior(debug_level(i), debug_level(debug_maxNtype + 1_pInt)) ! fill all debug types with levels specified by "all"
enddo enddo
if (iand(debug_level(debug_debug),debug_levelBasic) /= 0) then if (iand(debug_level(debug_debug),debug_levelBasic) /= 0) &
!$OMP CRITICAL (write2out) write(6,'(a,/)') ' using values from config file'
write(6,*) 'using values from config file'
write(6,*)
!$OMP END CRITICAL (write2out)
endif
! no config file, so we use standard values ! no config file, so we use standard values
else else
if (iand(debug_level(debug_debug),debug_levelBasic) /= 0) then if (iand(debug_level(debug_debug),debug_levelBasic) /= 0) &
!$OMP CRITICAL (write2out) write(6,'(a,/)') ' using standard values'
write(6,*) 'using standard values'
write(6,*)
!$OMP END CRITICAL (write2out)
endif
endif endif
!output switched on (debug level for debug must be extensive) !output switched on (debug level for debug must be extensive)
@ -302,7 +292,6 @@ subroutine debug_init
end select end select
if(debug_level(i) /= 0) then if(debug_level(i) /= 0) then
!$OMP CRITICAL (write2out)
write(6,'(a,a)') tag,' debugging:' write(6,'(a,a)') tag,' debugging:'
if(iand(debug_level(i),debug_levelBasic) /= 0) write(6,'(a)') ' basic' if(iand(debug_level(i),debug_levelBasic) /= 0) write(6,'(a)') ' basic'
if(iand(debug_level(i),debug_levelExtensive) /= 0) write(6,'(a)') ' extensive' if(iand(debug_level(i),debug_levelExtensive) /= 0) write(6,'(a)') ' extensive'
@ -317,7 +306,6 @@ subroutine debug_init
if(iand(debug_level(i),debug_spectralDivergence)/= 0) write(6,'(a)') ' divergence' if(iand(debug_level(i),debug_spectralDivergence)/= 0) write(6,'(a)') ' divergence'
if(iand(debug_level(i),debug_spectralRotation) /= 0) write(6,'(a)') ' rotation' if(iand(debug_level(i),debug_spectralRotation) /= 0) write(6,'(a)') ' rotation'
if(iand(debug_level(i),debug_spectralPETSc) /= 0) write(6,'(a)') ' PETSc' if(iand(debug_level(i),debug_spectralPETSc) /= 0) write(6,'(a)') ' PETSc'
!$OMP END CRITICAL (write2out)
endif endif
enddo enddo
endif endif

36
code/lattice.f90
View File

@ -96,10 +96,10 @@ module lattice
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! fcc (1) ! fcc (1)
integer(pInt), dimension(lattice_maxNslipFamily), parameter, private :: & integer(pInt), dimension(lattice_maxNslipFamily), parameter, public :: &
lattice_fcc_NslipSystem = int([12, 0, 0, 0, 0],pInt) lattice_fcc_NslipSystem = int([12, 0, 0, 0, 0],pInt)
integer(pInt), dimension(lattice_maxNtwinFamily), parameter, private :: & integer(pInt), dimension(lattice_maxNtwinFamily), parameter, public :: &
lattice_fcc_NtwinSystem = int([12, 0, 0, 0],pInt) lattice_fcc_NtwinSystem = int([12, 0, 0, 0],pInt)
integer(pInt), parameter, private :: & integer(pInt), parameter, private :: &
@ -157,7 +157,7 @@ module lattice
0.7071067812_pReal & 0.7071067812_pReal &
],[lattice_fcc_Ntwin]) !< Twin system <112>{111} Sorted according to Eisenlohr & Hantcherli ],[lattice_fcc_Ntwin]) !< Twin system <112>{111} Sorted according to Eisenlohr & Hantcherli
integer(pInt), dimension(lattice_fcc_Nslip,lattice_fcc_Nslip), target, private :: & integer(pInt), dimension(lattice_fcc_Nslip,lattice_fcc_Nslip), target, public :: &
lattice_fcc_interactionSlipSlip = reshape(int( [& lattice_fcc_interactionSlipSlip = reshape(int( [&
1,2,2,4,6,5,3,5,5,4,5,6, & ! ---> slip 1,2,2,4,6,5,3,5,5,4,5,6, & ! ---> slip
2,1,2,6,4,5,5,4,6,5,3,5, & ! | 2,1,2,6,4,5,5,4,6,5,3,5, & ! |
@ -180,7 +180,7 @@ module lattice
!< 5 --- glissile junctions !< 5 --- glissile junctions
!< 6 --- Lomer locks !< 6 --- Lomer locks
integer(pInt), dimension(lattice_fcc_Nslip,lattice_fcc_Ntwin), target, private :: & integer(pInt), dimension(lattice_fcc_Nslip,lattice_fcc_Ntwin), target, public :: &
lattice_fcc_interactionSlipTwin = reshape(int( [& lattice_fcc_interactionSlipTwin = reshape(int( [&
1,1,1,3,3,3,2,2,2,3,3,3, & ! ---> twin 1,1,1,3,3,3,2,2,2,3,3,3, & ! ---> twin
1,1,1,3,3,3,3,3,3,2,2,2, & ! | 1,1,1,3,3,3,3,3,3,2,2,2, & ! |
@ -200,10 +200,10 @@ module lattice
!< 2 --- screw trace between slip system and twin habit plane (easy cross slip) !< 2 --- screw trace between slip system and twin habit plane (easy cross slip)
!< 3 --- other interaction !< 3 --- other interaction
integer(pInt), dimension(lattice_fcc_Ntwin,lattice_fcc_Nslip), target, private :: & integer(pInt), dimension(lattice_fcc_Ntwin,lattice_fcc_Nslip), target, public :: &
lattice_fcc_interactionTwinSlip = 0_pInt lattice_fcc_interactionTwinSlip = 0_pInt
integer(pInt), dimension(lattice_fcc_Ntwin,lattice_fcc_Ntwin), target, private :: & integer(pInt), dimension(lattice_fcc_Ntwin,lattice_fcc_Ntwin), target, public :: &
lattice_fcc_interactionTwinTwin = reshape(int( [& lattice_fcc_interactionTwinTwin = reshape(int( [&
1,1,1,2,2,2,2,2,2,2,2,2, & ! ---> twin 1,1,1,2,2,2,2,2,2,2,2,2, & ! ---> twin
1,1,1,2,2,2,2,2,2,2,2,2, & ! | 1,1,1,2,2,2,2,2,2,2,2,2, & ! |
@ -230,10 +230,10 @@ module lattice
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! bcc (2) ! bcc (2)
integer(pInt), dimension(lattice_maxNslipFamily), parameter, private :: & integer(pInt), dimension(lattice_maxNslipFamily), parameter, public :: &
lattice_bcc_NslipSystem = int([ 12, 12, 0, 0, 0], pInt) lattice_bcc_NslipSystem = int([ 12, 12, 0, 0, 0], pInt)
integer(pInt), dimension(lattice_maxNtwinFamily), parameter, private :: & integer(pInt), dimension(lattice_maxNtwinFamily), parameter, public :: &
lattice_bcc_NtwinSystem = int([ 12, 0, 0, 0], pInt) lattice_bcc_NtwinSystem = int([ 12, 0, 0, 0], pInt)
integer(pInt), parameter, private :: & integer(pInt), parameter, private :: &
@ -334,7 +334,7 @@ module lattice
],[lattice_bcc_Ntwin]) ],[lattice_bcc_Ntwin])
! slip--slip interactions for BCC structures (2) from Lee et al. Int J Plast 15 (1999) 625-645 ! slip--slip interactions for BCC structures (2) from Lee et al. Int J Plast 15 (1999) 625-645
integer(pInt), dimension(lattice_bcc_Nslip,lattice_bcc_Nslip), target, private :: & integer(pInt), dimension(lattice_bcc_Nslip,lattice_bcc_Nslip), target, public :: &
lattice_bcc_interactionSlipSlip = reshape(int( [& lattice_bcc_interactionSlipSlip = reshape(int( [&
1,3,6,6,5,4,4,2,4,2,5,4, 6,6,4,2,2,4,6,6,4,2,6,6, & ! ---> slip 1,3,6,6,5,4,4,2,4,2,5,4, 6,6,4,2,2,4,6,6,4,2,6,6, & ! ---> slip
3,1,6,6,4,2,5,4,5,4,4,2, 6,6,2,4,4,2,6,6,2,4,6,6, & ! | 3,1,6,6,4,2,5,4,5,4,4,2, 6,6,2,4,4,2,6,6,2,4,6,6, & ! |
@ -370,7 +370,7 @@ module lattice
!< 5 --- weak sessile interaction !< 5 --- weak sessile interaction
!< 6 --- strong sessile interaction !< 6 --- strong sessile interaction
integer(pInt), dimension(lattice_bcc_Nslip,lattice_bcc_Ntwin), target, private :: & integer(pInt), dimension(lattice_bcc_Nslip,lattice_bcc_Ntwin), target, public :: &
lattice_bcc_interactionSlipTwin = reshape(int( [& lattice_bcc_interactionSlipTwin = reshape(int( [&
3,3,3,2,2,3,3,3,3,2,3,3, & ! ---> twin 3,3,3,2,2,3,3,3,3,2,3,3, & ! ---> twin
3,3,2,3,3,2,3,3,2,3,3,3, & ! | 3,3,2,3,3,2,3,3,2,3,3,3, & ! |
@ -404,11 +404,11 @@ module lattice
!< 3 --- other interaction !< 3 --- other interaction
! twin--slip interactions for BCC structures (2) MISSING: not implemented yet ! twin--slip interactions for BCC structures (2) MISSING: not implemented yet
integer(pInt), dimension(lattice_bcc_Ntwin,lattice_bcc_Nslip), target, private :: & integer(pInt), dimension(lattice_bcc_Ntwin,lattice_bcc_Nslip), target, public :: &
lattice_bcc_interactionTwinSlip = 0_pInt lattice_bcc_interactionTwinSlip = 0_pInt
! twin--twin interactions for BCC structures (2) ! twin--twin interactions for BCC structures (2)
integer(pInt), dimension(lattice_bcc_Ntwin,lattice_bcc_Ntwin), target, private :: & integer(pInt), dimension(lattice_bcc_Ntwin,lattice_bcc_Ntwin), target, public :: &
lattice_bcc_interactionTwinTwin = reshape(int( [& lattice_bcc_interactionTwinTwin = reshape(int( [&
1,3,3,3,3,3,3,2,3,3,2,3, & ! ---> twin 1,3,3,3,3,3,3,2,3,3,2,3, & ! ---> twin
3,1,3,3,3,3,2,3,3,3,3,2, & ! | 3,1,3,3,3,3,2,3,3,3,3,2, & ! |
@ -438,10 +438,10 @@ module lattice
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! hex (3+) ! hex (3+)
integer(pInt), dimension(lattice_maxNslipFamily), parameter, private :: & integer(pInt), dimension(lattice_maxNslipFamily), parameter, public :: &
lattice_hex_NslipSystem = int([ 3, 3, 6, 12, 6],pInt) lattice_hex_NslipSystem = int([ 3, 3, 6, 12, 6],pInt)
integer(pInt), dimension(lattice_maxNtwinFamily), parameter, private :: & integer(pInt), dimension(lattice_maxNtwinFamily), parameter, public :: &
lattice_hex_NtwinSystem = int([ 6, 6, 6, 6],pInt) lattice_hex_NtwinSystem = int([ 6, 6, 6, 6],pInt)
integer(pInt), parameter , private :: & integer(pInt), parameter , private :: &
@ -555,7 +555,7 @@ module lattice
!* 3. twin-twin interaction - 20 types !* 3. twin-twin interaction - 20 types
!* 4. twin-slip interaction - 16 types !* 4. twin-slip interaction - 16 types
integer(pInt), dimension(lattice_hex_Nslip,lattice_hex_Nslip), target, private :: & integer(pInt), dimension(lattice_hex_Nslip,lattice_hex_Nslip), target, public :: &
lattice_hex_interactionSlipSlip = reshape(int( [& lattice_hex_interactionSlipSlip = reshape(int( [&
1, 6, 6, 11,11,11, 15,15,15,15,15,15, 18,18,18,18,18,18,18,18,18,18,18,18, 20,20,20,20,20,20, & ! ---> slip 1, 6, 6, 11,11,11, 15,15,15,15,15,15, 18,18,18,18,18,18,18,18,18,18,18,18, 20,20,20,20,20,20, & ! ---> slip
6, 1, 6, 11,11,11, 15,15,15,15,15,15, 18,18,18,18,18,18,18,18,18,18,18,18, 20,20,20,20,20,20, & ! | 6, 1, 6, 11,11,11, 15,15,15,15,15,15, 18,18,18,18,18,18,18,18,18,18,18,18, 20,20,20,20,20,20, & ! |
@ -594,7 +594,7 @@ module lattice
],pInt),[lattice_hex_Nslip,lattice_hex_Nslip],order=[2,1]) ],pInt),[lattice_hex_Nslip,lattice_hex_Nslip],order=[2,1])
!* isotropic interaction at the moment !* isotropic interaction at the moment
integer(pInt), dimension(lattice_hex_Nslip,lattice_hex_Ntwin), target, private :: & integer(pInt), dimension(lattice_hex_Nslip,lattice_hex_Ntwin), target, public :: &
lattice_hex_interactionSlipTwin = reshape(int( [& lattice_hex_interactionSlipTwin = reshape(int( [&
1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, & ! --> twin 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, & ! --> twin
1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, & ! | 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, & ! |
@ -633,7 +633,7 @@ module lattice
],pInt),[lattice_hex_Nslip,lattice_hex_Ntwin],order=[2,1]) ],pInt),[lattice_hex_Nslip,lattice_hex_Ntwin],order=[2,1])
!* isotropic interaction at the moment !* isotropic interaction at the moment
integer(pInt), dimension(lattice_hex_Ntwin,lattice_hex_Nslip), target, private :: & integer(pInt), dimension(lattice_hex_Ntwin,lattice_hex_Nslip), target, public :: &
lattice_hex_interactionTwinSlip = reshape(int( [& lattice_hex_interactionTwinSlip = reshape(int( [&
1, 1, 1, 5, 5, 5, 9, 9, 9, 9, 9, 9, 13,13,13,13,13,13,13,13,13,13,13,13, 17,17,17,17,17,17, & ! --> slip 1, 1, 1, 5, 5, 5, 9, 9, 9, 9, 9, 9, 13,13,13,13,13,13,13,13,13,13,13,13, 17,17,17,17,17,17, & ! --> slip
1, 1, 1, 5, 5, 5, 9, 9, 9, 9, 9, 9, 13,13,13,13,13,13,13,13,13,13,13,13, 17,17,17,17,17,17, & ! | 1, 1, 1, 5, 5, 5, 9, 9, 9, 9, 9, 9, 13,13,13,13,13,13,13,13,13,13,13,13, 17,17,17,17,17,17, & ! |
@ -665,7 +665,7 @@ module lattice
],pInt),[lattice_hex_Ntwin,lattice_hex_Nslip],order=[2,1]) ],pInt),[lattice_hex_Ntwin,lattice_hex_Nslip],order=[2,1])
integer(pInt), dimension(lattice_hex_Ntwin,lattice_hex_Ntwin), target, private :: & integer(pInt), dimension(lattice_hex_Ntwin,lattice_hex_Ntwin), target, public :: &
lattice_hex_interactionTwinTwin = reshape(int( [& lattice_hex_interactionTwinTwin = reshape(int( [&
1, 5, 5, 5, 5, 5, 9, 9, 9, 9, 9, 9, 12,12,12,12,12,12, 14,14,14,14,14,14, & ! ---> twin 1, 5, 5, 5, 5, 5, 9, 9, 9, 9, 9, 9, 12,12,12,12,12,12, 14,14,14,14,14,14, & ! ---> twin
5, 1, 5, 5, 5, 5, 9, 9, 9, 9, 9, 9, 12,12,12,12,12,12, 14,14,14,14,14,14, & ! | 5, 1, 5, 5, 5, 5, 9, 9, 9, 9, 9, 9, 12,12,12,12,12,12, 14,14,14,14,14,14, & ! |

40
code/mesh.f90
View File

@ -426,11 +426,10 @@ subroutine mesh_init(ip,el)
implicit none implicit none
integer(pInt), parameter :: fileUnit = 222_pInt integer(pInt), parameter :: fileUnit = 222_pInt
integer(pInt), intent(in) :: el, ip integer(pInt), intent(in) :: el, ip
integer(pInt) :: e integer(pInt) :: j
write(6,*) write(6,'(/,a)') ' <<<+- mesh init -+>>>'
write(6,*) '<<<+- mesh init -+>>>' write(6,'(a)') ' $Id$'
write(6,*) '$Id$'
#include "compilation_info.f90" #include "compilation_info.f90"
if (allocated(mesh_mapFEtoCPelem)) deallocate(mesh_mapFEtoCPelem) if (allocated(mesh_mapFEtoCPelem)) deallocate(mesh_mapFEtoCPelem)
@ -463,12 +462,17 @@ subroutine mesh_init(ip,el)
! scale dimension to calculate either uncorrected, dimension-independent, or dimension- and reso- ! scale dimension to calculate either uncorrected, dimension-independent, or dimension- and reso-
! lution-independent divergence ! lution-independent divergence
if (divergence_correction == 1_pInt) then if (divergence_correction == 1_pInt) then
do e = 1_pInt, 3_pInt do j = 1_pInt, 3_pInt
if (e /= minloc(geomdim,1) .and. e /= maxloc(geomdim,1)) scaledDim = geomdim/geomdim(e) if (j /= minloc(geomdim,1) .and. j /= maxloc(geomdim,1)) scaledDim = geomdim/geomdim(j)
enddo enddo
elseif (divergence_correction == 2_pInt) then elseif (divergence_correction == 2_pInt) then
do e = 1_pInt, 3_pInt do j = 1_pInt, 3_pInt
if (e /= minloc(geomdim/res,1) .and. e /= maxloc(geomdim/res,1)) scaledDim = geomdim/geomdim(e)*res(e) if (j /= minloc(geomdim/res,1) .and. j /= maxloc(geomdim/res,1)) scaledDim = geomdim/geomdim(j)*res(j)
enddo
elseif (divergence_correction == 3_pInt) then
do j = 1_pInt, 3_pInt
if (j/=minloc(geomdim/sqrt(real(res,pReal)),1) .and. j/=maxloc(geomdim/sqrt(real(res,pReal)),1))&
scaledDim=geomdim/geomdim(j)*sqrt(real(res(j),pReal))
enddo enddo
else else
scaledDim = geomdim scaledDim = geomdim
@ -530,7 +534,7 @@ subroutine mesh_init(ip,el)
FEsolving_execElem = [ 1_pInt,mesh_NcpElems ] FEsolving_execElem = [ 1_pInt,mesh_NcpElems ]
if (allocated(FEsolving_execIP)) deallocate(FEsolving_execIP) if (allocated(FEsolving_execIP)) deallocate(FEsolving_execIP)
allocate(FEsolving_execIP(2_pInt,mesh_NcpElems)); FEsolving_execIP = 1_pInt allocate(FEsolving_execIP(2_pInt,mesh_NcpElems)); FEsolving_execIP = 1_pInt
forall (e = 1_pInt:mesh_NcpElems) FEsolving_execIP(2,e) = FE_Nips(FE_geomtype(mesh_element(2,e))) forall (j = 1_pInt:mesh_NcpElems) FEsolving_execIP(2,j) = FE_Nips(FE_geomtype(mesh_element(2,j)))
if (allocated(calcMode)) deallocate(calcMode) if (allocated(calcMode)) deallocate(calcMode)
allocate(calcMode(mesh_maxNips,mesh_NcpElems)) allocate(calcMode(mesh_maxNips,mesh_NcpElems))
@ -809,7 +813,7 @@ function mesh_spectral_getResolution(fileUnit)
read(myUnit,'(a1024)') line read(myUnit,'(a1024)') line
positions = IO_stringPos(line,7_pInt) positions = IO_stringPos(line,7_pInt)
keyword = IO_lc(IO_StringValue(line,positions,2_pInt)) keyword = IO_lc(IO_StringValue(line,positions,2_pInt,.true.))
if (keyword(1:4) == 'head') then if (keyword(1:4) == 'head') then
headerLength = IO_intValue(line,positions,1_pInt) + 1_pInt headerLength = IO_intValue(line,positions,1_pInt) + 1_pInt
else else
@ -819,7 +823,7 @@ function mesh_spectral_getResolution(fileUnit)
do i = 1_pInt, headerLength do i = 1_pInt, headerLength
read(myUnit,'(a1024)') line read(myUnit,'(a1024)') line
positions = IO_stringPos(line,7_pInt) positions = IO_stringPos(line,7_pInt)
select case ( IO_lc(IO_StringValue(line,positions,1_pInt)) ) select case ( IO_lc(IO_StringValue(line,positions,1_pInt,.true.)) )
case ('resolution') case ('resolution')
gotResolution = .true. gotResolution = .true.
do j = 2_pInt,6_pInt,2_pInt do j = 2_pInt,6_pInt,2_pInt
@ -890,7 +894,7 @@ function mesh_spectral_getDimension(fileUnit)
read(myUnit,'(a1024)') line read(myUnit,'(a1024)') line
positions = IO_stringPos(line,7_pInt) positions = IO_stringPos(line,7_pInt)
keyword = IO_lc(IO_StringValue(line,positions,2_pInt)) keyword = IO_lc(IO_StringValue(line,positions,2_pInt,.true.))
if (keyword(1:4) == 'head') then if (keyword(1:4) == 'head') then
headerLength = IO_intValue(line,positions,1_pInt) + 1_pInt headerLength = IO_intValue(line,positions,1_pInt) + 1_pInt
else else
@ -900,7 +904,7 @@ function mesh_spectral_getDimension(fileUnit)
do i = 1_pInt, headerLength do i = 1_pInt, headerLength
read(myUnit,'(a1024)') line read(myUnit,'(a1024)') line
positions = IO_stringPos(line,7_pInt) positions = IO_stringPos(line,7_pInt)
select case ( IO_lc(IO_StringValue(line,positions,1)) ) select case ( IO_lc(IO_StringValue(line,positions,1,.true.)) )
case ('dimension') case ('dimension')
gotDimension = .true. gotDimension = .true.
do j = 2_pInt,6_pInt,2_pInt do j = 2_pInt,6_pInt,2_pInt
@ -964,7 +968,7 @@ function mesh_spectral_getHomogenization(fileUnit)
read(myUnit,'(a1024)') line read(myUnit,'(a1024)') line
positions = IO_stringPos(line,7_pInt) positions = IO_stringPos(line,7_pInt)
keyword = IO_lc(IO_StringValue(line,positions,2_pInt)) keyword = IO_lc(IO_StringValue(line,positions,2_pInt,.true.))
if (keyword(1:4) == 'head') then if (keyword(1:4) == 'head') then
headerLength = IO_intValue(line,positions,1_pInt) + 1_pInt headerLength = IO_intValue(line,positions,1_pInt) + 1_pInt
else else
@ -974,7 +978,7 @@ function mesh_spectral_getHomogenization(fileUnit)
do i = 1_pInt, headerLength do i = 1_pInt, headerLength
read(myUnit,'(a1024)') line read(myUnit,'(a1024)') line
positions = IO_stringPos(line,7_pInt) positions = IO_stringPos(line,7_pInt)
select case ( IO_lc(IO_StringValue(line,positions,1)) ) select case ( IO_lc(IO_StringValue(line,positions,1,.true.)) )
case ('homogenization') case ('homogenization')
gotHomogenization = .true. gotHomogenization = .true.
mesh_spectral_getHomogenization = IO_intValue(line,positions,2_pInt) mesh_spectral_getHomogenization = IO_intValue(line,positions,2_pInt)
@ -989,6 +993,8 @@ function mesh_spectral_getHomogenization(fileUnit)
call IO_error(error_ID = 842_pInt, ext_msg='mesh_spectral_getHomogenization') call IO_error(error_ID = 842_pInt, ext_msg='mesh_spectral_getHomogenization')
end function mesh_spectral_getHomogenization end function mesh_spectral_getHomogenization
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Count overall number of nodes and elements in mesh and stores them in !> @brief Count overall number of nodes and elements in mesh and stores them in
!! 'mesh_Nelems' and 'mesh_Nnodes' !! 'mesh_Nelems' and 'mesh_Nnodes'
@ -1129,7 +1135,7 @@ subroutine mesh_spectral_build_elements(myUnit)
read(myUnit,'(a65536)') line read(myUnit,'(a65536)') line
myPos = IO_stringPos(line,7_pInt) myPos = IO_stringPos(line,7_pInt)
keyword = IO_lc(IO_StringValue(line,myPos,2_pInt)) keyword = IO_lc(IO_StringValue(line,myPos,2_pInt,.true.))
if (keyword(1:4) == 'head') then if (keyword(1:4) == 'head') then
headerLength = IO_intValue(line,myPos,1_pInt) + 1_pInt headerLength = IO_intValue(line,myPos,1_pInt) + 1_pInt
else else
@ -1735,7 +1741,7 @@ function mesh_deformedCoordsLinear(gDim,F,FavgIn) result(coords)
real(pReal), dimension(3,0:size(F,3)-1,0:size(F,4)-1,0:size(F,5)-1,0:7) :: & real(pReal), dimension(3,0:size(F,3)-1,0:size(F,4)-1,0:size(F,5)-1,0:7) :: &
coordsAvgOrder coordsAvgOrder
integer(pInt), parameter, dimension(3) :: & integer(pInt), parameter, dimension(3) :: &
iOnes = 1.0_pReal iOnes = 1_pInt
real(pReal), parameter, dimension(3) :: & real(pReal), parameter, dimension(3) :: &
fOnes = 1.0_pReal fOnes = 1.0_pReal
real(pReal), dimension(3) :: & real(pReal), dimension(3) :: &