From e644c6dbc5cf9a4d85322b1ca775cd6058d1d41d Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 8 Feb 2013 15:55:53 +0000 Subject: [PATCH] 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 --- code/IO.f90 | 126 +++++++++++++++------------- code/constitutive_phenopowerlaw.f90 | 59 +++++++++---- code/debug.f90 | 28 ++----- code/lattice.f90 | 36 ++++---- code/mesh.f90 | 40 +++++---- 5 files changed, 159 insertions(+), 130 deletions(-) diff --git a/code/IO.f90 b/code/IO.f90 index 8f5a0eaff..0a4f6fc2a 100644 --- a/code/IO.f90 +++ b/code/IO.f90 @@ -844,17 +844,25 @@ end function IO_stringPos !-------------------------------------------------------------------------------------------------- !> @brief read string value at myPos from line !-------------------------------------------------------------------------------------------------- -function IO_stringValue(line,positions,myPos) +function IO_stringValue(line,positions,myPos,silent) implicit none integer(pInt), dimension(:), intent(in) :: positions integer(pInt), intent(in) :: myPos character(len=1+positions(myPos*2+1)-positions(myPos*2)) :: IO_stringValue character(len=*), intent(in) :: line - - if (myPos > positions(1)) then - IO_stringValue = '' - call IO_warning(201, e=myPos, ext_msg = trim(line)//' (IO_stringValue)') + logical, optional,intent(in) :: silent + logical :: warn + + if (.not. present(silent)) then + warn = .false. + 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 IO_stringValue = line(positions(myPos*2):positions(myPos*2+1)) endif @@ -868,7 +876,6 @@ end function IO_stringValue pure function IO_fixedStringValue (line,ends,myPos) implicit none - integer(pInt), intent(in) :: myPos integer(pInt), dimension(:), intent(in) :: ends 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 integer(pInt), dimension(:), intent(in) :: positions 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.+-' integer(pInt) :: readStatus, invalidWhere IO_floatValue = 0.0_pReal - - 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)//')') + if (myPos > positions(1) .or. myPos < 1_pInt) then ! trying to access non-present value + call IO_warning(201,ext_msg=myName//trim(line)) else - invalidWhere = verify(line(positions(myPos*2):positions(myPos*2+1)),validCharacters) - if (invalidWhere /= 0_pInt) then - invalidWhere = invalidWhere-1 - call IO_warning(202,ext_msg=line(positions(myPos*2):positions(myPos*2+1))//' ('//trim(myName)//')') + invalidWhere = verify(line(positions(myPos*2):positions(myPos*2+1)),validCharacters) ! search for invalid characters + if (invalidWhere /= 0_pInt) then ! found invaldid character, only read in substring + invalidWhere = invalidWhere - 1_pInt + call IO_warning(202,ext_msg=myName//line(positions(myPos*2):positions(myPos*2+1))) 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 - read(UNIT=line(positions(myPos*2):positions(myPos*2)+invalidWhere-1),iostat=readStatus,FMT=*) IO_floatValue - if (readStatus /= 0_pInt) & - call IO_warning(203,ext_msg=line(positions(myPos*2):positions(myPos*2)+invalidWhere-1)//' ('//trim(myName)//')') + read(UNIT=line(positions(myPos*2):positions(myPos*2)+invalidWhere-1),iostat=readStatus,FMT=*) & + IO_floatValue + 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 end function IO_floatValue @@ -921,24 +928,23 @@ real(pReal) function IO_fixedFloatValue (line,ends,myPos) character(len=*), intent(in) :: line integer(pInt), intent(in) :: myPos 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.+-' integer(pInt) :: readStatus, myStart, invalidWhere IO_fixedFloatValue = 0.0_pReal - - myStart = ends(myPos-1)+1 + myStart = ends(myPos-1) + 1_pInt - invalidWhere = verify(line(myStart:ends(myPos)),validCharacters) - if (invalidWhere /= 0_pInt) then - invalidWhere = invalidWhere-1 - call IO_warning(202,ext_msg=line(myStart:ends(myPos))//' ('//trim(myName)//')') + invalidWhere = verify(line(myStart:ends(myPos)),validCharacters) ! search for invalid character + if (invalidWhere /= 0_pInt) then ! found invaldid character, only read in substring + invalidWhere = invalidWhere - 1_pInt + call IO_warning(202,ext_msg=myName//line(myStart:ends(myPos))) else - invalidWhere = ends(myPos)-myStart+1 + invalidWhere = ends(myPos)-myStart + 1_pInt ! read until ends(myPos) endif read(UNIT=line(myStart:myStart+invalidWhere-1),iostat=readStatus,FMT=*) IO_fixedFloatValue - if (readStatus /= 0_pInt) & - call IO_warning(203,ext_msg=line(myStart:myStart+invalidWhere-1)//' ('//trim(myName)//')') + if (readStatus /= 0_pInt) & ! error during string to float conversion + call IO_warning(203,ext_msg=myName//line(myStart:myStart+invalidWhere-1)) end function IO_fixedFloatValue @@ -952,42 +958,45 @@ real(pReal) function IO_fixedNoEFloatValue (line,ends,myPos) character(len=*), intent(in) :: line integer(pInt), intent(in) :: myPos 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=12), parameter :: validExp = '0123456789+-' - integer(pInt) :: expon = 0, myStart, readStatus + integer(pInt) :: expon, myStart, readStatus integer :: pos_exp, end_base, end_exp - real(pReal) :: base = 0.0_pReal + real(pReal) :: base + + base = 0.0_pReal + expon = 0_pInt - myStart = ends(myPos-1)+1 + myStart = ends(myPos-1) + 1_pInt pos_exp = scan(line(myStart:ends(myPos)),'+-',back=.true.) - if (pos_exp <= 1_pInt) & ! no exponent but only base - pos_exp = ends(myPos)-myStart+1 + if (pos_exp <= 1_pInt) & ! no exponent but only base + pos_exp = ends(myPos)-myStart + 1_pInt ! --- figure out base --- - end_base = verify(line(myStart:myStart+pos_exp-1),validBase) - if (end_base /= 0_pInt) then ! invalid base + end_base = verify(line(myStart:myStart+pos_exp-1),validBase) ! search for invalid character in base + if (end_base /= 0_pInt) then ! found invaldid character, only read in substring 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 - end_base = pos_exp + end_base = pos_exp ! read until begin of exponent endif read(UNIT=line(myStart:myStart+end_base-1),iostat=readStatus,FMT=*) base - if (readStatus /= 0_pInt) & - call IO_warning(203, ext_msg = line(myStart:myStart+end_base-1)//' ('//trim(myName)//':base)') + if (readStatus /= 0_pInt) & ! error during string to float conversion + call IO_warning(203, ext_msg = myName//'(base): '//line(myStart:myStart+end_base-1)) ! --- figure out exponent --- - end_exp = verify(line(myStart+pos_exp:ends(myPos)),validExp) - if (end_exp /= 0_pInt) then ! invalid exponent - end_exp = end_exp-1 - call IO_warning(202, ext_msg = line(myStart+pos_exp:ends(myPos))//' ('//trim(myName)//':exp)') + end_exp = verify(line(myStart+pos_exp:ends(myPos)),validExp) ! search for invalid character in exponent + if (end_exp /= 0_pInt) then ! found invaldid character, only read in substring + end_exp = end_exp - 1_pInt + call IO_warning(202, ext_msg = myName//'(exp): '//line(myStart+pos_exp:ends(myPos))) else - end_exp = mystart-ends(myPos)+1 + end_exp = mystart-ends(myPos) + 1_pInt ! read until end of string endif read(UNIT=line(myStart+pos_exp:myStart+end_exp-1),iostat=readStatus,FMT=*) expon - if (readStatus /= 0_pInt) & - call IO_warning(203, ext_msg = line(myStart+pos_exp:myStart+end_exp-1)//' ('//trim(myName)//':exp)') + if (readStatus /= 0_pInt) & ! error during string to float conversion + call IO_warning(203, ext_msg = myName//'(base): '//line(myStart+pos_exp:myStart+end_exp-1)) IO_fixedNoEFloatValue = base*10.0_pReal**expon @@ -1003,25 +1012,26 @@ integer(pInt) function IO_intValue(line,positions,myPos) character(len=*), intent(in) :: line integer(pInt), dimension(:), intent(in) :: positions 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+-' integer(pInt) :: readStatus, invalidWhere IO_intValue = 0_pInt - 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)//')') + if (myPos > positions(1) .or. myPos < 1_pInt) then ! trying to access non-present value + call IO_warning(201,ext_msg=myName//trim(line)) else 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 - 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 invalidWhere = positions(myPos*2+1)-positions(myPos*2)+1 endif - read(UNIT=line(positions(myPos*2):positions(myPos*2)+invalidWhere-1),iostat=readStatus,FMT=*) IO_intValue - if (readStatus /= 0_pInt) & - call IO_warning(203,ext_msg=line(positions(myPos*2):positions(myPos*2)+invalidWhere-1)//' ('//trim(myName)//')') + read(UNIT=line(positions(myPos*2):positions(myPos*2)+invalidWhere-1),iostat=readStatus,FMT=*) & + IO_intValue + 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 end function IO_intValue @@ -1036,7 +1046,7 @@ integer(pInt) function IO_fixedIntValue(line,ends,myPos) character(len=*), intent(in) :: line integer(pInt), intent(in) :: myPos 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.+-' integer(pInt) :: readStatus, myStart, invalidWhere @@ -1047,13 +1057,13 @@ integer(pInt) function IO_fixedIntValue(line,ends,myPos) invalidWhere = verify(line(myStart:ends(myPos)),validCharacters) if (invalidWhere /= 0_pInt) then 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 invalidWhere = ends(myPos)-myStart+1 endif read(UNIT=line(myStart:myStart+invalidWhere-1),iostat=readStatus,FMT=*) IO_fixedIntValue 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 diff --git a/code/constitutive_phenopowerlaw.f90 b/code/constitutive_phenopowerlaw.f90 index 137e49f0a..f878b9a13 100644 --- a/code/constitutive_phenopowerlaw.f90 +++ b/code/constitutive_phenopowerlaw.f90 @@ -125,20 +125,15 @@ subroutine constitutive_phenopowerlaw_init(myFile) use debug, only: debug_level,& debug_constitutive,& debug_levelBasic - use lattice, only: lattice_initializeStructure, lattice_symmetrizeC66, & - lattice_maxNslipFamily, lattice_maxNtwinFamily, & - lattice_maxNinteraction, lattice_NslipSystem, lattice_NtwinSystem, & - lattice_maxNonSchmid, & - lattice_interactionSlipSlip, & - lattice_interactionSlipTwin, & - lattice_interactionTwinSlip, & - lattice_interactionTwinTwin + use lattice implicit none integer(pInt), intent(in) :: myFile integer(pInt), parameter :: maxNchunks = lattice_maxNinteraction + 1_pInt integer(pInt), dimension(1+2*maxNchunks) :: positions 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 character(len=64) :: tag 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) 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 write(6,'(a16,1x,i5)') '# instances:',maxNinstance write(6,*) @@ -262,6 +264,29 @@ subroutine constitutive_phenopowerlaw_init(myFile) IO_lc(IO_stringValue(line,positions,2_pInt)) case ('lattice_structure') 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') constitutive_phenopowerlaw_CoverA(i) = IO_floatValue(line,positions,2_pInt) case ('c11') @@ -283,7 +308,7 @@ subroutine constitutive_phenopowerlaw_init(myFile) case ('c66') constitutive_phenopowerlaw_Cslip_66(6,6,i) = IO_floatValue(line,positions,2_pInt) 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) enddo case ('gdot0_slip') @@ -291,17 +316,17 @@ subroutine constitutive_phenopowerlaw_init(myFile) case ('n_slip') constitutive_phenopowerlaw_n_slip(i) = IO_floatValue(line,positions,2_pInt) 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) enddo 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) enddo case ('a_slip', 'w0_slip') constitutive_phenopowerlaw_a_slip(i) = IO_floatValue(line,positions,2_pInt) 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) enddo case ('gdot0_twin') @@ -309,7 +334,7 @@ subroutine constitutive_phenopowerlaw_init(myFile) case ('n_twin') constitutive_phenopowerlaw_n_twin(i) = IO_floatValue(line,positions,2_pInt) 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) enddo case ('s_pr') @@ -338,19 +363,19 @@ subroutine constitutive_phenopowerlaw_init(myFile) case ('atol_twinfrac') constitutive_phenopowerlaw_aTolTwinfrac(i) = IO_floatValue(line,positions,2_pInt) 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) enddo 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) enddo 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) enddo 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) enddo case ('nonschmid_coefficients') diff --git a/code/debug.f90 b/code/debug.f90 index cf81df5b8..7e5ce98e0 100644 --- a/code/debug.f90 +++ b/code/debug.f90 @@ -144,12 +144,9 @@ subroutine debug_init integer(pInt), dimension(1+2*maxNchunks) :: positions character(len=64) :: tag character(len=1024) :: line - !$OMP CRITICAL (write2out) - write(6,*) - write(6,*) '<<<+- debug init -+>>>' - write(6,*) '$Id$' + write(6,'(/,a)') ' <<<+- debug init -+>>>' + write(6,'(a)') ' $Id$' #include "compilation_info.f90" - !$OMP END CRITICAL (write2out) if (allocated(debug_StressLoopDistribution)) & deallocate(debug_StressLoopDistribution) @@ -223,7 +220,7 @@ subroutine debug_init what = debug_maxNtype + 2_pInt end select if(what /= 0) then - do i = 2_pInt, maxNchunks + do i = 2_pInt, positions(1) select case(IO_lc(IO_stringValue(line,positions,i))) case('basic') 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" enddo - if (iand(debug_level(debug_debug),debug_levelBasic) /= 0) then - !$OMP CRITICAL (write2out) - write(6,*) 'using values from config file' - write(6,*) - !$OMP END CRITICAL (write2out) - endif + if (iand(debug_level(debug_debug),debug_levelBasic) /= 0) & + write(6,'(a,/)') ' using values from config file' + ! no config file, so we use standard values else - if (iand(debug_level(debug_debug),debug_levelBasic) /= 0) then - !$OMP CRITICAL (write2out) - write(6,*) 'using standard values' - write(6,*) - !$OMP END CRITICAL (write2out) - endif + if (iand(debug_level(debug_debug),debug_levelBasic) /= 0) & + write(6,'(a,/)') ' using standard values' endif !output switched on (debug level for debug must be extensive) @@ -302,7 +292,6 @@ subroutine debug_init end select if(debug_level(i) /= 0) then - !$OMP CRITICAL (write2out) write(6,'(a,a)') tag,' debugging:' if(iand(debug_level(i),debug_levelBasic) /= 0) write(6,'(a)') ' basic' 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_spectralRotation) /= 0) write(6,'(a)') ' rotation' if(iand(debug_level(i),debug_spectralPETSc) /= 0) write(6,'(a)') ' PETSc' - !$OMP END CRITICAL (write2out) endif enddo endif diff --git a/code/lattice.f90 b/code/lattice.f90 index 4a9a63d59..82b931153 100644 --- a/code/lattice.f90 +++ b/code/lattice.f90 @@ -96,10 +96,10 @@ module lattice !-------------------------------------------------------------------------------------------------- ! 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) - integer(pInt), dimension(lattice_maxNtwinFamily), parameter, private :: & + integer(pInt), dimension(lattice_maxNtwinFamily), parameter, public :: & lattice_fcc_NtwinSystem = int([12, 0, 0, 0],pInt) integer(pInt), parameter, private :: & @@ -157,7 +157,7 @@ module lattice 0.7071067812_pReal & ],[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( [& 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, & ! | @@ -180,7 +180,7 @@ module lattice !< 5 --- glissile junctions !< 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( [& 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, & ! | @@ -200,10 +200,10 @@ module lattice !< 2 --- screw trace between slip system and twin habit plane (easy cross slip) !< 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 - 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( [& 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, & ! | @@ -230,10 +230,10 @@ module lattice !-------------------------------------------------------------------------------------------------- ! 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) - integer(pInt), dimension(lattice_maxNtwinFamily), parameter, private :: & + integer(pInt), dimension(lattice_maxNtwinFamily), parameter, public :: & lattice_bcc_NtwinSystem = int([ 12, 0, 0, 0], pInt) integer(pInt), parameter, private :: & @@ -334,7 +334,7 @@ module lattice ],[lattice_bcc_Ntwin]) ! 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( [& 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, & ! | @@ -370,7 +370,7 @@ module lattice !< 5 --- weak 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( [& 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, & ! | @@ -404,11 +404,11 @@ module lattice !< 3 --- other interaction ! 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 ! 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( [& 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, & ! | @@ -438,10 +438,10 @@ module lattice !-------------------------------------------------------------------------------------------------- ! 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) - integer(pInt), dimension(lattice_maxNtwinFamily), parameter, private :: & + integer(pInt), dimension(lattice_maxNtwinFamily), parameter, public :: & lattice_hex_NtwinSystem = int([ 6, 6, 6, 6],pInt) integer(pInt), parameter , private :: & @@ -555,7 +555,7 @@ module lattice !* 3. twin-twin interaction - 20 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( [& 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, & ! | @@ -594,7 +594,7 @@ module lattice ],pInt),[lattice_hex_Nslip,lattice_hex_Nslip],order=[2,1]) !* 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( [& 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, & ! | @@ -633,7 +633,7 @@ module lattice ],pInt),[lattice_hex_Nslip,lattice_hex_Ntwin],order=[2,1]) !* 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( [& 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, & ! | @@ -665,7 +665,7 @@ module lattice ],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( [& 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, & ! | diff --git a/code/mesh.f90 b/code/mesh.f90 index 43a9e3c0b..3a579498e 100644 --- a/code/mesh.f90 +++ b/code/mesh.f90 @@ -426,11 +426,10 @@ subroutine mesh_init(ip,el) implicit none integer(pInt), parameter :: fileUnit = 222_pInt integer(pInt), intent(in) :: el, ip - integer(pInt) :: e + integer(pInt) :: j - write(6,*) - write(6,*) '<<<+- mesh init -+>>>' - write(6,*) '$Id$' + write(6,'(/,a)') ' <<<+- mesh init -+>>>' + write(6,'(a)') ' $Id$' #include "compilation_info.f90" 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- ! lution-independent divergence if (divergence_correction == 1_pInt) then - do e = 1_pInt, 3_pInt - if (e /= minloc(geomdim,1) .and. e /= maxloc(geomdim,1)) scaledDim = geomdim/geomdim(e) + do j = 1_pInt, 3_pInt + if (j /= minloc(geomdim,1) .and. j /= maxloc(geomdim,1)) scaledDim = geomdim/geomdim(j) enddo elseif (divergence_correction == 2_pInt) then - do e = 1_pInt, 3_pInt - if (e /= minloc(geomdim/res,1) .and. e /= maxloc(geomdim/res,1)) scaledDim = geomdim/geomdim(e)*res(e) + do j = 1_pInt, 3_pInt + 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 else scaledDim = geomdim @@ -530,7 +534,7 @@ subroutine mesh_init(ip,el) FEsolving_execElem = [ 1_pInt,mesh_NcpElems ] if (allocated(FEsolving_execIP)) deallocate(FEsolving_execIP) 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) allocate(calcMode(mesh_maxNips,mesh_NcpElems)) @@ -809,7 +813,7 @@ function mesh_spectral_getResolution(fileUnit) read(myUnit,'(a1024)') line 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 headerLength = IO_intValue(line,positions,1_pInt) + 1_pInt else @@ -819,7 +823,7 @@ function mesh_spectral_getResolution(fileUnit) do i = 1_pInt, headerLength read(myUnit,'(a1024)') line 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') gotResolution = .true. do j = 2_pInt,6_pInt,2_pInt @@ -890,7 +894,7 @@ function mesh_spectral_getDimension(fileUnit) read(myUnit,'(a1024)') line 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 headerLength = IO_intValue(line,positions,1_pInt) + 1_pInt else @@ -900,7 +904,7 @@ function mesh_spectral_getDimension(fileUnit) do i = 1_pInt, headerLength read(myUnit,'(a1024)') line 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') gotDimension = .true. do j = 2_pInt,6_pInt,2_pInt @@ -964,7 +968,7 @@ function mesh_spectral_getHomogenization(fileUnit) read(myUnit,'(a1024)') line 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 headerLength = IO_intValue(line,positions,1_pInt) + 1_pInt else @@ -974,7 +978,7 @@ function mesh_spectral_getHomogenization(fileUnit) do i = 1_pInt, headerLength read(myUnit,'(a1024)') line 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') gotHomogenization = .true. 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') end function mesh_spectral_getHomogenization + + !-------------------------------------------------------------------------------------------------- !> @brief Count overall number of nodes and elements in mesh and stores them in !! 'mesh_Nelems' and 'mesh_Nnodes' @@ -1129,7 +1135,7 @@ subroutine mesh_spectral_build_elements(myUnit) read(myUnit,'(a65536)') line 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 headerLength = IO_intValue(line,myPos,1_pInt) + 1_pInt 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) :: & coordsAvgOrder integer(pInt), parameter, dimension(3) :: & - iOnes = 1.0_pReal + iOnes = 1_pInt real(pReal), parameter, dimension(3) :: & fOnes = 1.0_pReal real(pReal), dimension(3) :: &