diff --git a/code/DAMASK_spectral.f90 b/code/DAMASK_spectral.f90 index c36cccf88..b855aaff0 100644 --- a/code/DAMASK_spectral.f90 +++ b/code/DAMASK_spectral.f90 @@ -53,7 +53,6 @@ program DAMASK_spectral use numerics, only: err_div_tol, err_stress_tolrel, rotation_tol, itmax, & memory_efficient, update_gamma, & simplified_algorithm, divergence_correction, & - cut_off_value, & DAMASK_NumThreadsInt, & fftw_planner_flag, fftw_timelimit use homogenization, only: materialpoint_sizeResults, materialpoint_results @@ -121,18 +120,9 @@ program DAMASK_spectral s0_reference real(pReal), dimension(6) :: cstress ! cauchy stress real(pReal), dimension(6,6) :: dsde, c0_66, s0_66 ! small strain stiffness - real(pReal), dimension(9,9) :: s_prev99, c_prev99, c0_99, s0_99 ! compliance and stiffness in matrix notation + real(pReal), dimension(9,9) :: s_prev99, c_prev99 ! compliance and stiffness in matrix notation real(pReal), dimension(:,:), allocatable :: s_reduced, c_reduced ! reduced compliance and stiffness (only for stress BC) - real(pReal), dimension(6,6) :: mask_inversion = reshape([& - 1.0_pReal, 1.0_pReal, 1.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal,& - 1.0_pReal, 1.0_pReal, 1.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal,& - 1.0_pReal, 1.0_pReal, 1.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal,& - 0.0_pReal, 0.0_pReal, 0.0_pReal, 1.0_pReal, 0.0_pReal, 0.0_pReal,& - 0.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal, 1.0_pReal, 0.0_pReal,& - 0.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal, 1.0_pReal],& - [ 6_pInt, 6_pInt]) - real(pReal), dimension(3,3,3,3) :: temp_3333 = 0.0_pReal - integer(pInt) :: size_reduced = 0_pInt ! number of stress BCs + integer(pInt) :: size_reduced = 0_pInt ! number of stress BCs !-------------------------------------------------------------------------------------------------- ! pointwise data @@ -151,7 +141,7 @@ program DAMASK_spectral real(pReal), dimension(3,3) :: xiDyad ! product of wave vectors real(pReal), dimension(:,:,:,:,:,:,:), allocatable :: gamma_hat ! gamma operator (field) for spectral method real(pReal), dimension(:,:,:,:), allocatable :: xi ! wave vector field for divergence and for gamma operator - integer(pInt), dimension(3) :: k_s, cutting_freq + integer(pInt), dimension(3) :: k_s !-------------------------------------------------------------------------------------------------- ! loop variables, convergence etc. @@ -174,7 +164,7 @@ program DAMASK_spectral !-------------------------------------------------------------------------------------------------- !variables for additional output due to general debugging - real(pReal) :: defgradDetMax, defgradDetMin, maxCorrectionSym, maxCorrectionSkew, max_diag, max_offdiag + real(pReal) :: defgradDetMax, defgradDetMin, maxCorrectionSym, maxCorrectionSkew !-------------------------------------------------------------------------------------------------- ! variables for additional output of divergence calculations @@ -378,8 +368,6 @@ program DAMASK_spectral res1_red = res(1)/2_pInt + 1_pInt ! size of complex array in first dimension (c2r, r2c) Npoints = res(1)*res(2)*res(3) wgt = 1.0_pReal/real(Npoints, pReal) - if (cut_off_value <0.0_pReal .or. cut_off_value >0.9_pReal) stop - cutting_freq = nint(real(res,pReal)*cut_off_value,pInt) ! for cut_off_value=0.0 just the highest freq. is removed !-------------------------------------------------------------------------------------------------- ! output of geometry @@ -394,7 +382,6 @@ program DAMASK_spectral print '(a,3(i12 ))','resolution a b c:', res print '(a,3(f12.5))','dimension x y z:', geomdim print '(a,i5)','homogenization: ',homog - if(cut_off_value/=0.0_pReal) print '(a,3(i12),a)', 'cutting away ', cutting_freq, ' frequencies' print '(a)', '#############################################################' print '(a,a)', 'loadcase file: ',trim(getLoadcaseName()) @@ -421,8 +408,8 @@ program DAMASK_spectral write (*,'(3(3(f12.7,1x)/))',advance='no') merge(math_transpose33(bc(loadcase)%deformation),& reshape(spread(DAMASK_NaN,1,9),[ 3,3]),transpose(bc(loadcase)%maskDeformation)) write (*,'(a,/,3(3(f12.7,1x)/))',advance='no') ' stress / GPa:',& - 1e-9*merge(math_transpose33(bc(loadcase)%stress),reshape(spread(DAMASK_NaN,1,9),[ 3,3])& - ,transpose(bc(loadcase)%maskStress)) + 1e-9_pReal*merge(math_transpose33(bc(loadcase)%stress),& + reshape(spread(DAMASK_NaN,1,9),[ 3,3]),transpose(bc(loadcase)%maskStress)) if (any(bc(loadcase)%rotation /= math_I3)) & write (*,'(a,/,3(3(f12.7,1x)/))',advance='no') ' rotation of loadframe:',& math_transpose33(bc(loadcase)%rotation) @@ -473,7 +460,7 @@ program DAMASK_spectral ielem = ielem + 1_pInt defgrad(i,j,k,1:3,1:3) = math_I3 defgradold(i,j,k,1:3,1:3) = math_I3 - coordinates(i,j,k,1:3) = geomdim/real(res, pReal)*[i,j,k] - geomdim/real(2_pInt*res,pReal) + coordinates(i,j,k,1:3) = geomdim/real(res * [i,j,k], pReal) - geomdim/real(2_pInt*res,pReal) call CPFEM_general(2_pInt,coordinates(i,j,k,1:3),math_I3,math_I3,temperature(i,j,k),& 0.0_pReal,ielem,1_pInt,cstress,dsde,pstress,dPdF) c_current = c_current + dPdF @@ -511,16 +498,16 @@ program DAMASK_spectral do j = 1_pInt, res(2) k_s(2) = j - 1_pInt if(j > res(2)/2_pInt + 1_pInt) k_s(2) = k_s(2) - res(2) - do i = 1, res1_red + do i = 1_pInt, res1_red k_s(1) = i - 1_pInt xi(1:3,i,j,k) = real(k_s, pReal)/geomdim enddo; enddo; enddo !-------------------------------------------------------------------------------------------------- ! calculate the gamma operator - if(memory_efficient) then ! allocate just single fourth order tensor + if(memory_efficient) then ! allocate just single fourth order tensor allocate (gamma_hat(1,1,1,3,3,3,3)); gamma_hat = 0.0_pReal - else ! precalculation of gamma_hat field + else ! precalculation of gamma_hat field allocate (gamma_hat(res1_red ,res(2),res(3),3,3,3,3)); gamma_hat = 0.0_pReal do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res1_red ! if(k==res(3)/2 .or. k==res(3)/2+2 .or.& @@ -528,21 +515,22 @@ program DAMASK_spectral ! i==res(1)/2 .or. i==res(1)/2+2) then ! gamma_hat(i,j,k,1:3,1:3,1:3,1:3) = s0_reference ! else - forall(l = 1_pInt:3_pInt, m = 1_pInt:3_pInt) & - xiDyad(l,m) = xi(l, i,j,k)*xi(m, i,j,k) - forall(l = 1_pInt:3_pInt, m = 1_pInt:3_pInt) & - temp33_Real(l,m) = sum(c0_reference(l,m,1:3,1:3)*xiDyad) - temp33_Real = math_inv33(temp33_Real) - forall(l=1_pInt:3_pInt, m=1_pInt:3_pInt, n=1_pInt:3_pInt, p=1_pInt:3_pInt)& - gamma_hat(i,j,k, l,m,n,p) = temp33_Real(l,n)*xiDyad(m,p) - ! endif + if(any([i,j,k] /= 1_pInt)) then ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1 + forall(l = 1_pInt:3_pInt, m = 1_pInt:3_pInt) & + xiDyad(l,m) = xi(l, i,j,k)*xi(m, i,j,k) + forall(l = 1_pInt:3_pInt, m = 1_pInt:3_pInt) & + temp33_Real(l,m) = sum(c0_reference(l,m,1:3,1:3)*xiDyad) + temp33_Real = math_inv33(temp33_Real) + forall(l=1_pInt:3_pInt, m=1_pInt:3_pInt, n=1_pInt:3_pInt, p=1_pInt:3_pInt)& + gamma_hat(i,j,k, l,m,n,p) = temp33_Real(l,n)*xiDyad(m,p) + endif enddo; enddo; enddo gamma_hat(1,1,1, 1:3,1:3,1:3,1:3) = 0.0_pReal ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1 endif !-------------------------------------------------------------------------------------------------- ! general initialization of fftw (see manual on fftw.org for more details) - if (pReal /= C_DOUBLE .or. pInt /= C_INT) call IO_error(error_ID=808_pInt) ! check for correct precision in C + if (pReal /= C_DOUBLE .or. pInt /= C_INT) call IO_error(error_ID=808_pInt) ! check for correct precision in C #ifdef _OPENMP if(DAMASK_NumThreadsInt > 0_pInt) then ierr = fftw_init_threads() @@ -550,7 +538,7 @@ program DAMASK_spectral call fftw_plan_with_nthreads(DAMASK_NumThreadsInt) endif #endif - call fftw_set_timelimit(fftw_timelimit) ! set timelimit for plan creation + call fftw_set_timelimit(fftw_timelimit) ! set timelimit for plan creation !-------------------------------------------------------------------------------------------------- ! creating plans @@ -739,7 +727,7 @@ program DAMASK_spectral !-------------------------------------------------------------------------------------------------- ! report begin of new increment print '(a)', '##################################################################' - print '(A,I5.5,A,es12.6)', 'Increment ', totalIncsCounter, ' Time ',time + print '(A,I5.5,A,es12.5)', 'Increment ', totalIncsCounter, ' Time ',time guessmode = 1.0_pReal ! keep guessing along former trajectory during same loadcase CPFEM_mode = 1_pInt ! winding forward @@ -801,7 +789,7 @@ program DAMASK_spectral row = (mod(totalIncsCounter+iter-2_pInt,9_pInt))/3_pInt + 1_pInt ! go through the elements of the tensors, controlled by totalIncsCounter and iter, starting at 1 column = (mod(totalIncsCounter+iter-2_pInt,3_pInt)) + 1_pInt scalarField_real(1:res(1),1:res(2),1:res(3)) =& ! store the selected component - tensorField_real(1:res(1),1:res(2),1:res(3),row,column) + cmplx(tensorField_real(1:res(1),1:res(2),1:res(3),row,column),0.0_pReal,pReal) endif !-------------------------------------------------------------------------------------------------- @@ -829,7 +817,7 @@ program DAMASK_spectral pstress_av_lab = real(tensorField_fourier(1,1,1,1:3,1:3),pReal)*wgt pstress_av = math_rotate_forward33(pstress_av_lab,bc(loadcase)%rotation) write (*,'(a,/,3(3(f12.7,1x)/))',advance='no') 'Piola-Kirchhoff stress / MPa:',& - math_transpose33(pstress_av)/1.e6 + math_transpose33(pstress_av)/1.e6_pReal !-------------------------------------------------------------------------------------------------- ! comparing 1 and 3x3 FT results @@ -940,20 +928,17 @@ program DAMASK_spectral print '(a,es11.4)', 'error divergence FT max = ',err_div_max print '(a,es11.4)', 'error divergence Real RMS = ',err_real_div_RMS print '(a,es11.4)', 'error divergence Real max = ',err_real_div_max - print '(a,es11.4)', 'divergence RMS FT/real = ',err_div_RMS/err_real_div_RMS - print '(a,es11.4)', 'divergence max FT/real = ',err_div_max/err_real_div_max print '(a,es11.4)', 'max deviat. from postProc = ',max_div_error endif - print '(a,f6.2,a,es11.4,a)', 'error divergence = ', err_div/err_div_tol, ' (',err_div,' 1/m)' + print '(a,f6.2,a,es11.4,3a)','error divergence = ', err_div/err_div_tol, & + ' (',err_div_RMS,' N/m',char(179),')' !-------------------------------------------------------------------------------------------------- ! divergence is calculated from FT(stress), depending on algorithm use field for spectral method if (.not. simplified_algorithm) tensorField_fourier = tau_fourier - max_diag = tiny(1.0_pReal) - max_offdiag = tiny(1.0_pReal) !-------------------------------------------------------------------------------------------------- ! to the actual spectral method calculation (mechanical equilibrium) - if(memory_efficient) then ! memory saving version, on-the-fly calculation of gamma_hat + if(memory_efficient) then ! memory saving version, on-the-fly calculation of gamma_hat do k = 1_pInt, res(3); do j = 1_pInt, res(2) ;do i = 1_pInt, res1_red ! if(k==res(3)/2 .or. k==res(3)/2+2 .or.& @@ -962,34 +947,38 @@ program DAMASK_spectral ! forall( m = 1_pInt:3_pInt, n = 1_pInt:3_pInt)& ! temp33_Complex(m,n) = sum(s0_reference(m,n, 1:3,1:3)* tensorField_fourier(i,j,k,1:3,1:3)) ! else - forall(l = 1_pInt:3_pInt, m = 1_pInt:3_pInt) & - xiDyad(l,m) = xi(l, i,j,k)*xi(m, i,j,k) - forall(l = 1_pInt:3_pInt, m = 1_pInt:3_pInt) & - temp33_Real(l,m) = sum(c0_reference(l,m,1:3,1:3)*xiDyad) - temp33_Real = math_inv33(temp33_Real) - forall(l=1_pInt:3_pInt, m=1_pInt:3_pInt, n=1_pInt:3_pInt, p=1_pInt:3_pInt)& - gamma_hat(1,1,1, l,m,n,p) = temp33_Real(l,n)*xiDyad(m,p) - forall(l = 1_pInt:3_pInt, m = 1_pInt:3_pInt) & - temp33_Complex(l,m) = sum(gamma_hat(1,1,1, l,m, 1:3,1:3) * tensorField_fourier(i,j,k,1:3,1:3)) - tensorField_fourier(i,j,k,1:3,1:3) = temp33_Complex - ! endif + if(any([i,j,k] /= 1_pInt)) then ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1 + forall(l = 1_pInt:3_pInt, m = 1_pInt:3_pInt) & + xiDyad(l,m) = xi(l, i,j,k)*xi(m, i,j,k) + forall(l = 1_pInt:3_pInt, m = 1_pInt:3_pInt) & + temp33_Real(l,m) = sum(c0_reference(l,m,1:3,1:3)*xiDyad) + temp33_Real = math_inv33(temp33_Real) + forall(l=1_pInt:3_pInt, m=1_pInt:3_pInt, n=1_pInt:3_pInt, p=1_pInt:3_pInt)& + gamma_hat(1,1,1, l,m,n,p) = temp33_Real(l,n)*xiDyad(m,p) + forall(l = 1_pInt:3_pInt, m = 1_pInt:3_pInt) & + temp33_Complex(l,m) = sum(gamma_hat(1,1,1, l,m, 1:3,1:3) *& + tensorField_fourier(i,j,k,1:3,1:3)) + tensorField_fourier(i,j,k,1:3,1:3) = temp33_Complex + endif enddo; enddo; enddo - else ! use precalculated gamma-operator + else ! use precalculated gamma-operator do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt,res1_red forall( m = 1_pInt:3_pInt, n = 1_pInt:3_pInt) & - temp33_Complex(m,n) = sum(gamma_hat(i,j,k, m,n, 1:3,1:3) * tensorField_fourier(i,j,k,1:3,1:3)) + temp33_Complex(m,n) = sum(gamma_hat(i,j,k, m,n, 1:3,1:3) *& + tensorField_fourier(i,j,k,1:3,1:3)) tensorField_fourier(i,j,k, 1:3,1:3) = temp33_Complex enddo; enddo; enddo endif - if (simplified_algorithm) then ! do not use the polarization field based algorithm - tensorField_fourier(1,1,1,1:3,1:3) = (defgrad_av_lab - defgradAim_lab) & ! assign (negative) average deformation gradient change to zero frequency (real part) - * real(Npoints,pReal) ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1 + if (simplified_algorithm) then ! do not use the polarization field based algorithm + tensorField_fourier(1,1,1,1:3,1:3) = cmplx((defgrad_av_lab - defgradAim_lab) & ! assign (negative) average deformation gradient change to zero frequency (real part) + * real(Npoints,pReal),0.0_pReal,pReal) ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1 else - tensorField_fourier(1,1,1,1:3,1:3) = defgradAim_lab * real(Npoints,pReal) ! assign deformation aim to zero frequency (real part) + tensorField_fourier(1,1,1,1:3,1:3) = cmplx(defgradAim_lab*real(Npoints,pReal),& ! assign deformation aim to zero frequency (real part) + 0.0_pReal,pReal) endif !-------------------------------------------------------------------------------------------------- @@ -998,7 +987,7 @@ program DAMASK_spectral do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res1_red scalarField_fourier(i,j,k) = tensorField_fourier(i,j,k,row,column) enddo; enddo; enddo - do i = 0_pInt, res(1)/2_pInt-2_pInt !unpack fft data for conj complex symmetric part. can be directly used in calculation of cstress_field + do i = 0_pInt, res(1)/2_pInt-2_pInt !unpack fft data for conj complex symmetric part m = 1_pInt do k = 1_pInt, res(3) n = 1_pInt @@ -1136,11 +1125,13 @@ end program DAMASK_spectral ! quit subroutine to satisfy IO_error ! !******************************************************************** -subroutine quit(id) +subroutine quit(stop_id) use prec implicit none - integer(pInt) id + integer(pInt), intent(in) :: stop_id + print*, stop_id + + stop 'abnormal termination of DAMASK_spectral' - stop end subroutine diff --git a/code/DAMASK_spectral_interface.f90 b/code/DAMASK_spectral_interface.f90 index deda85309..877b32767 100644 --- a/code/DAMASK_spectral_interface.f90 +++ b/code/DAMASK_spectral_interface.f90 @@ -39,7 +39,7 @@ subroutine DAMASK_interface_init() implicit none character(len=1024) commandLine, hostName, userName - integer(pInt):: i, start = 0_pInt, length=0_pInt + integer :: i, start = 0, length=0 integer, dimension(8) :: date_and_time_values ! type default integer call get_command(commandLine) call DATE_AND_TIME(VALUES=date_and_time_values) @@ -47,7 +47,7 @@ subroutine DAMASK_interface_init() if(640_pInt .or. index(commandLine,' --help ',.true.)>0_pInt) then ! search for ' -h ' or '--help' + if(index(commandLine,' -h ',.true.)>0 .or. index(commandLine,' --help ',.true.)>0) then ! search for ' -h ' or '--help' write(6,*) '$Id$' #include "compilation_info.f90" print '(a)', '#############################################################' @@ -87,12 +87,12 @@ subroutine DAMASK_interface_init() endif if (.not.(command_argument_count()==4 .or. command_argument_count()==6)) & ! check for correct number of given arguments (no --help) stop 'Wrong Nr. of Arguments. Run DAMASK_spectral.exe --help' ! Could not find valid keyword (position 0 +3). Functions from IO.f90 are not available - start = index(commandLine,'-g',.true.) + 3_pInt ! search for '-g' and jump to first char of geometry + start = index(commandLine,'-g',.true.) + 3 ! search for '-g' and jump to first char of geometry if (index(commandLine,'--geom',.true.)>0) then ! if '--geom' is found, use that (contains '-g') - start = index(commandLine,'--geom',.true.) + 7_pInt + start = index(commandLine,'--geom',.true.) + 7 endif if (index(commandLine,'--geometry',.true.)>0) then ! again, now searching for --geometry' - start = index(commandLine,'--geometry',.true.) + 11_pInt + start = index(commandLine,'--geometry',.true.) + 11 endif if(start==3_pInt) stop 'No Geometry specified, terminating DAMASK'! Could not find valid keyword (position 0 +3). Functions from IO.f90 are not available length = index(commandLine(start:len(commandLine)),' ',.false.) @@ -105,12 +105,12 @@ subroutine DAMASK_interface_init() if(640) then ! if '--load' is found, use that (contains '-l') - start = index(commandLine,'--load',.true.) + 7_pInt + start = index(commandLine,'--load',.true.) + 7 endif if (index(commandLine,'--loadcase',.true.)>0) then ! again, now searching for --loadcase' - start = index(commandLine,'--loadcase',.true.) + 11_pInt + start = index(commandLine,'--loadcase',.true.) + 11 endif if(start==3_pInt) stop 'No Loadcase specified, terminating DAMASK'! Could not find valid keyword (position 0 +3). Functions from IO.f90 are not available length = index(commandLine(start:len(commandLine)),' ',.false.) @@ -123,9 +123,9 @@ subroutine DAMASK_interface_init() if(640) then ! if '--restart' is found, use that (contains '-l') - start = index(commandLine,'--restart',.true.) + 7_pInt + start = index(commandLine,'--restart',.true.) + 7 endif length = index(commandLine(start:len(commandLine)),' ',.false.) @@ -201,12 +201,12 @@ function getModelName() character(1024) getModelName, cwd character(len=*), parameter :: pathSep = achar(47)//achar(92) ! forwardslash, backwardslash - integer(pInt) :: posExt,posSep + integer :: posExt,posSep posExt = scan(geometryParameter,'.',back=.true.) posSep = scan(geometryParameter,pathSep,back=.true.) - if (posExt <= posSep) posExt = len_trim(geometryParameter)+1_pInt ! no extension present + if (posExt <= posSep) posExt = len_trim(geometryParameter)+1 ! no extension present getModelName = geometryParameter(1:posExt-1_pInt) ! path to geometry file (excl. extension) if (scan(getModelName,pathSep) /= 1) then ! relative path given as command line argument @@ -227,19 +227,17 @@ endfunction getModelName !******************************************************************** function getLoadCase() - use prec, only: pInt - implicit none - character(1024) getLoadCase + character(1024) :: getLoadCase character(len=*), parameter :: pathSep = achar(47)//achar(92) ! forwardslash, backwardslash - integer(pInt) posExt,posSep + integer :: posExt,posSep posExt = scan(loadcaseParameter,'.',back=.true.) posSep = scan(loadcaseParameter,pathSep,back=.true.) - if (posExt <= posSep) posExt = len_trim(loadcaseParameter)+1_pInt ! no extension present - getLoadCase = loadcaseParameter(posSep+1_pInt:posExt-1_pInt) ! name of load case file exluding extension + if (posExt <= posSep) posExt = len_trim(loadcaseParameter)+1 ! no extension present + getLoadCase = loadcaseParameter(posSep+1:posExt-1) ! name of load case file exluding extension endfunction getLoadCase @@ -254,11 +252,9 @@ function getLoadcaseName() implicit none - character(len=1024) getLoadcaseName,cwd + character(len=1024) :: getLoadcaseName,cwd character(len=*), parameter :: pathSep = achar(47)//achar(92) ! forwardslash, backwardslash - integer(pInt) posExt,posSep - posExt = 0_pInt - + integer :: posExt = 0, posSep getLoadcaseName = loadcaseParameter posExt = scan(getLoadcaseName,'.',back=.true.) posSep = scan(getLoadcaseName,pathSep,back=.true.) @@ -286,31 +282,31 @@ function rectifyPath(path) implicit none - character(len=*) path - character(len=len_trim(path)) rectifyPath - integer(pInt) i,j,k,l + character(len=*) :: path + character(len=len_trim(path)) :: rectifyPath + integer :: i,j,k,l !no pInt !remove ./ from path l = len_trim(path) rectifyPath = path - do i = l,3_pInt,-1_pInt - if ( rectifyPath(i-1_pInt:i) == './' .and. rectifyPath(i-2_pInt:i-2_pInt) /= '.' ) & - rectifyPath(i-1_pInt:l) = rectifyPath(i+1_pInt:l)//' ' + do i = l,3,-1 + if ( rectifyPath(i-1:i) == './' .and. rectifyPath(i-2:i-2) /= '.' ) & + rectifyPath(i-1:l) = rectifyPath(i+1:l)//' ' enddo !remove ../ and corresponding directory from rectifyPath l = len_trim(rectifyPath) i = index(rectifyPath(i:l),'../') - j = 0_pInt + j = 0 do while (i > j) - j = scan(rectifyPath(1:i-2_pInt),'/',back=.true.) - rectifyPath(j+1_pInt:l) = rectifyPath(i+3_pInt:l)//repeat(' ',2_pInt+i-j) - if (rectifyPath(j+1_pInt:j+1_pInt) == '/') then !search for '//' that appear in case of XXX/../../XXX + j = scan(rectifyPath(1:i-2),'/',back=.true.) + rectifyPath(j+1:l) = rectifyPath(i+3:l)//repeat(' ',2+i-j) + if (rectifyPath(j+1:j+1) == '/') then !search for '//' that appear in case of XXX/../../XXX k = len_trim(rectifyPath) - rectifyPath(j+1_pInt:k-1_pInt) = rectifyPath(j+2_pInt:k) + rectifyPath(j+1:k-1) = rectifyPath(j+2:k) rectifyPath(k:k) = ' ' endif - i = j+index(rectifyPath(j+1_pInt:l),'../') + i = j+index(rectifyPath(j+1:l),'../') enddo if(len_trim(rectifyPath) == 0) rectifyPath = '/' @@ -330,18 +326,18 @@ function makeRelativePath(a,b) character (len=*) :: a,b character (len=1024) :: makeRelativePath - integer(pInt) i,posLastCommonSlash,remainingSlashes + integer :: i,posLastCommonSlash,remainingSlashes !no pInt - posLastCommonSlash = 0_pInt - remainingSlashes = 0_pInt - do i = 1_pInt,min(1024,len_trim(a),len_trim(b)) + posLastCommonSlash = 0 + remainingSlashes = 0 + do i = 1, min(1024,len_trim(a),len_trim(b)) if (a(i:i) /= b(i:i)) exit if (a(i:i) == '/') posLastCommonSlash = i enddo - do i = posLastCommonSlash+1_pInt,len_trim(a) - if (a(i:i) == '/') remainingSlashes = remainingSlashes + 1_pInt + do i = posLastCommonSlash+1,len_trim(a) + if (a(i:i) == '/') remainingSlashes = remainingSlashes + 1 enddo - makeRelativePath = repeat('../',remainingSlashes)//b(posLastCommonSlash+1_pInt:len_trim(b)) + makeRelativePath = repeat('../',remainingSlashes)//b(posLastCommonSlash+1:len_trim(b)) endfunction makeRelativePath diff --git a/code/FEsolving.f90 b/code/FEsolving.f90 index 3d8e03bc8..4cfd09c64 100644 --- a/code/FEsolving.f90 +++ b/code/FEsolving.f90 @@ -70,9 +70,9 @@ commandLine(i:i) = achar(iachar(commandLine(i:i))+32) ! make lowercase enddo if (index(commandLine,'-r ',.true.)>0) & ! look for -r - start = index(commandLine,'-r ',.true.) + 3_pInt ! set to position after trailing space + start = index(commandLine,'-r ',.true.) + 3 ! set to position after trailing space if (index(commandLine,'--restart ',.true.)>0) & ! look for --restart - start = index(commandLine,'--restart ',.true.) + 10_pInt ! set to position after trailing space + start = index(commandLine,'--restart ',.true.) + 10 ! set to position after trailing space if(start /= 0_pInt) then ! found something length = verify(commandLine(start:len(commandLine)),'0123456789',.false.) ! where is first non number after argument? read(commandLine(start:start+length),'(I12)') restartInc ! read argument @@ -99,12 +99,12 @@ restartWrite = iand(IO_intValue(line,positions,1_pInt),1_pInt) > 0_pInt restartRead = iand(IO_intValue(line,positions,1_pInt),2_pInt) > 0_pInt case ('*restart') - do i=2,positions(1) - restartWrite = (IO_lc(IO_StringValue(line,positions,i)) == 'write') .or. restartWrite - restartRead = (IO_lc(IO_StringValue(line,positions,i)) == 'read') .or. restartRead + do j=2_pInt,positions(1) + restartWrite = (IO_lc(IO_StringValue(line,positions,j)) == 'write') .or. restartWrite + restartRead = (IO_lc(IO_StringValue(line,positions,j)) == 'read') .or. restartRead enddo if(restartWrite) then - do j=2,positions(1) + do j=2_pInt,positions(1) restartWrite = (IO_lc(IO_StringValue(line,positions,j)) /= 'frequency=0') .and. restartWrite enddo endif diff --git a/code/IO.f90 b/code/IO.f90 index da523ca6c..908e01fa8 100644 --- a/code/IO.f90 +++ b/code/IO.f90 @@ -25,20 +25,20 @@ CONTAINS !--------------------------- ! function IO_abaqus_assembleInputFile -! subroutine IO_open_file(unit,relPath) -! subroutine IO_open_inputFile(unit, model) -! subroutine IO_open_logFile(unit) +! subroutine IO_open_file(myUnit,relPath) +! subroutine IO_open_inputFile(myUnit, model) +! subroutine IO_open_logFile(myUnit) ! function IO_hybridIA(Nast,ODFfileName) ! private function hybridIA_reps(dV_V,steps,C) ! function IO_stringPos(line,maxN) -! function IO_stringValue(line,positions,pos) -! function IO_floatValue(line,positions,pos) -! function IO_intValue(line,positions,pos) -! function IO_fixedStringValue(line,ends,pos) -! function IO_fixedFloatValue(line,ends,pos) -! function IO_fixedFloatNoEValue(line,ends,pos) -! function IO_fixedIntValue(line,ends,pos) -! function IO_continousIntValues(unit,maxN) +! function IO_stringValue(line,positions,myPos) +! function IO_floatValue(line,positions,myPos) +! function IO_intValue(line,positions,myPos) +! function IO_fixedStringValue(line,ends,myPos) +! function IO_fixedFloatValue(line,ends,myPos) +! function IO_fixedFloatNoEValue(line,ends,myPos) +! function IO_fixedIntValue(line,ends,myPos) +! function IO_continousIntValues(myUnit,maxN) ! function IO_lc(line) ! subroutine IO_lcInplace(line) ! subroutine IO_error(ID) @@ -76,7 +76,7 @@ recursive function IO_abaqus_assembleInputFile(unit1,unit2) result(createSuccess character(len=300) line,fname integer(pInt), intent(in) :: unit1, unit2 logical createSuccess,fexist - integer(pInt), parameter :: maxNchunks = 6 + integer(pInt), parameter :: maxNchunks = 6_pInt integer(pInt), dimension(1+2*maxNchunks) :: positions @@ -88,7 +88,7 @@ recursive function IO_abaqus_assembleInputFile(unit1,unit2) result(createSuccess ! call IO_lcInPlace(line) if (IO_lc(IO_StringValue(line,positions,1_pInt))=='*include') then - fname = trim(getSolverWorkingDirectoryName())//trim(line(9_pInt+scan(line(9_pInt:),'='):)) + fname = trim(getSolverWorkingDirectoryName())//trim(line(9+scan(line(9:),'='):)) inquire(file=fname, exist=fexist) if (.not.(fexist)) then !$OMP CRITICAL (write2out) @@ -121,25 +121,25 @@ end function !*********************************************************** ! check if the input file for Abaqus contains part info !*********************************************************** - function IO_abaqus_hasNoPart(unit) + function IO_abaqus_hasNoPart(myUnit) use prec, only: pInt implicit none - integer(pInt) unit - integer(pInt), parameter :: maxNchunks = 1 - integer(pInt), dimension(1+2*maxNchunks) :: pos + integer(pInt) myUnit + integer(pInt), parameter :: maxNchunks = 1_pInt + integer(pInt), dimension(1+2*maxNchunks) :: myPos logical IO_abaqus_hasNoPart character(len=300) line IO_abaqus_hasNoPart = .true. 610 FORMAT(A300) - rewind(unit) + rewind(myUnit) do - read(unit,610,END=620) line - pos = IO_stringPos(line,maxNchunks) - if (IO_lc(IO_stringValue(line,pos,1_pInt)) == '*part' ) then + read(myUnit,610,END=620) line + myPos = IO_stringPos(line,maxNchunks) + if (IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*part' ) then IO_abaqus_hasNoPart = .false. exit endif @@ -150,22 +150,22 @@ end function !******************************************************************** -! open existing file to given unit +! open existing file to given myUnit ! path to file is relative to working directory !******************************************************************** - logical function IO_open_file_stat(unit,relPath) + logical function IO_open_file_stat(myUnit,relPath) use prec, only: pInt use DAMASK_interface implicit none - integer(pInt), intent(in) :: unit + integer(pInt), intent(in) :: myUnit character(len=*), intent(in) :: relPath character(len=1024) path integer(pInt) stat path = trim(getSolverWorkingDirectoryName())//relPath - open(unit,status='old',iostat=stat,file=path) + open(myUnit,status='old',iostat=stat,file=path) IO_open_file_stat = (stat == 0_pInt) endfunction @@ -175,37 +175,37 @@ end function ! open existing file to given unit ! path to file is relative to working directory !******************************************************************** - subroutine IO_open_file(unit,relPath) + subroutine IO_open_file(myUnit,relPath) use prec, only: pInt use DAMASK_interface implicit none - integer(pInt), intent(in) :: unit + integer(pInt), intent(in) ::myUnit character(len=*), intent(in) :: relPath character(len=1024) path integer(pInt) stat path = trim(getSolverWorkingDirectoryName())//relPath - open(unit,status='old',iostat=stat,file=path) + open(myUnit,status='old',iostat=stat,file=path) if (stat /= 0_pInt) call IO_error(100_pInt,ext_msg=path) endsubroutine !******************************************************************** -! open FEM inputfile to given unit +! open FEM inputfile to given myUnit ! AP: 12.07.10 ! : changed the function to open *.inp_assembly, which is basically ! the input file without comment lines and possibly assembled includes !******************************************************************** - subroutine IO_open_inputFile(unit,model) + subroutine IO_open_inputFile(myUnit,model) use prec, only: pReal, pInt use DAMASK_interface implicit none - integer(pInt), intent(in) :: unit + integer(pInt), intent(in) :: myUnit character(len=*), intent(in) :: model character(len=1024) path integer(pInt) stat @@ -214,19 +214,18 @@ end function if (FEsolver == 'Abaqus') then path = trim(getSolverWorkingDirectoryName())//trim(model)//InputFileExtension - open(unit+1,status='old',iostat=stat,file=path) + open(myUnit+1,status='old',iostat=stat,file=path) if (stat /= 0_pInt) call IO_error(100_pInt,ext_msg=path) path = trim(getSolverWorkingDirectoryName())//trim(model)//InputFileExtension//'_assembly' - open(unit,iostat=stat,file=path) + open(myUnit,iostat=stat,file=path) if (stat /= 0_pInt) call IO_error(100_pInt,ext_msg=path) - if (IO_abaqus_assembleInputFile(unit,unit+1_pInt)) call IO_error(103_pInt) ! strip comments and concatenate any "include"s - close(unit+1_pInt) + if (IO_abaqus_assembleInputFile(myUnit,myUnit+1_pInt)) call IO_error(103_pInt) ! strip comments and concatenate any "include"s + close(myUnit+1_pInt) else - path = trim(getSolverWorkingDirectoryName())//trim(model)//InputFileExtension - open(unit,status='old',iostat=stat,file=path) + open(myUnit,status='old',iostat=stat,file=path) if (stat /= 0_pInt) call IO_error(100_pInt,ext_msg=path) endif @@ -235,20 +234,19 @@ end function !******************************************************************** -! open FEM logfile to given unit +! open FEM logfile to given myUnit !******************************************************************** - subroutine IO_open_logFile(unit) - + subroutine IO_open_logFile(myUnit) use prec, only: pReal, pInt use DAMASK_interface implicit none - integer(pInt), intent(in) :: unit + integer(pInt), intent(in) :: myUnit character(len=1024) path integer(pInt) stat path = trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//LogFileExtension - open(unit,status='old',iostat=stat,file=path) + open(myUnit,status='old',iostat=stat,file=path) if (stat /= 0) call IO_error(100_pInt,ext_msg=path) endsubroutine @@ -256,23 +254,23 @@ end function !******************************************************************** ! open (write) file related to current job -! but with different extension to given unit +! but with different extension to given myUnit !******************************************************************** - logical function IO_open_jobFile_stat(unit,newExt) + logical function IO_open_jobFile_stat(myUnit,newExt) use prec, only: pReal, pInt use DAMASK_interface implicit none - integer(pInt), intent(in) :: unit + integer(pInt), intent(in) :: myUnit character(len=*), intent(in) :: newExt character(len=1024) path integer(pInt) stat path = trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//'.'//newExt - open(unit,status='old',iostat=stat,file=path) + open(myUnit,status='old',iostat=stat,file=path) IO_open_jobFile_stat = (stat == 0_pInt) - + endfunction @@ -280,19 +278,19 @@ end function ! open (write) file related to current job ! but with different extension to given unit !******************************************************************** - subroutine IO_open_jobFile(unit,newExt) + subroutine IO_open_jobFile(myUnit,newExt) use prec, only: pReal, pInt use DAMASK_interface implicit none - integer(pInt), intent(in) :: unit + integer(pInt), intent(in) :: myUnit character(len=*), intent(in) :: newExt character(len=1024) path integer(pInt) stat path = trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//'.'//newExt - open(unit,status='old',iostat=stat,file=path) + open(myUnit,status='old',iostat=stat,file=path) if (stat /= 0_pInt) call IO_error(100_pInt,ext_msg=path) endsubroutine @@ -300,21 +298,21 @@ end function !******************************************************************** ! open (write) file related to current job -! but with different extension to given unit +! but with different extension to given myUnit !******************************************************************** - subroutine IO_write_jobFile(unit,newExt) + subroutine IO_write_jobFile(myUnit,newExt) use prec, only: pReal, pInt use DAMASK_interface implicit none - integer(pInt), intent(in) :: unit + integer(pInt), intent(in) :: myUnit character(len=*), intent(in) :: newExt character(len=1024) path integer(pInt) stat path = trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//'.'//newExt - open(unit,status='replace',iostat=stat,file=path) + open(myUnit,status='replace',iostat=stat,file=path) if (stat /= 0_pInt) call IO_error(100_pInt,ext_msg=path) endsubroutine @@ -322,15 +320,15 @@ end function !******************************************************************** ! open (write) binary file related to current job -! but with different extension to given unit +! but with different extension to given myUnit !******************************************************************** - subroutine IO_write_jobBinaryFile(unit,newExt,recMultiplier) + subroutine IO_write_jobBinaryFile(myUnit,newExt,recMultiplier) use prec, only: pReal, pInt use DAMASK_interface implicit none - integer(pInt), intent(in) :: unit + integer(pInt), intent(in) :: myUnit integer(pInt), intent(in), optional :: recMultiplier character(len=*), intent(in) :: newExt character(len=1024) path @@ -338,9 +336,9 @@ end function path = trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//'.'//newExt if (present(recMultiplier)) then - open(unit,status='replace',form='unformatted',access='direct',recl=pReal*recMultiplier,iostat=stat,file=path) + open(myUnit,status='replace',form='unformatted',access='direct',recl=pReal*recMultiplier,iostat=stat,file=path) else - open(unit,status='replace',form='unformatted',access='direct',recl=pReal,iostat=stat,file=path) + open(myUnit,status='replace',form='unformatted',access='direct',recl=pReal,iostat=stat,file=path) endif if (stat /= 0_pInt) call IO_error(100_pInt,ext_msg=path) @@ -349,15 +347,15 @@ end function !******************************************************************** ! open (read) binary file related to restored job -! and with different extension to given unit +! and with different extension to given myUnit !******************************************************************** - subroutine IO_read_jobBinaryFile(unit,newExt,jobName,recMultiplier) + subroutine IO_read_jobBinaryFile(myUnit,newExt,jobName,recMultiplier) use prec, only: pReal, pInt use DAMASK_interface implicit none - integer(pInt), intent(in) :: unit + integer(pInt), intent(in) :: myUnit integer(pInt), intent(in), optional :: recMultiplier character(len=*), intent(in) :: newExt, jobName character(len=1024) path @@ -365,9 +363,9 @@ end function path = trim(getSolverWorkingDirectoryName())//trim(jobName)//'.'//newExt if (present(recMultiplier)) then - open(unit,status='old',form='unformatted',access='direct',recl=pReal*recMultiplier,iostat=stat,file=path) + open(myUnit,status='old',form='unformatted',access='direct',recl=pReal*recMultiplier,iostat=stat,file=path) else - open(unit,status='old',form='unformatted',access='direct',recl=pReal,iostat=stat,file=path) + open(myUnit,status='old',form='unformatted',access='direct',recl=pReal,iostat=stat,file=path) endif if (stat /= 0) then call IO_error(100_pInt,ext_msg=path) @@ -390,9 +388,9 @@ end function real(pReal), intent(in) :: C hybridIA_reps = 0_pInt - do phi1=1,steps(1) - do Phi =1,steps(2) - do phi2=1,steps(3) + do phi1=1_pInt,steps(1) + do Phi =1_pInt,steps(2) + do phi2=1_pInt,steps(3) hybridIA_reps = hybridIA_reps+nint(C*dV_V(phi2,Phi,phi1), pInt) enddo enddo @@ -416,7 +414,7 @@ end function character(len=80) line character(len=*), parameter :: fileFormat = '(A80)' integer(pInt) i,j,bin,Nast,NnonZero,Nset,Nreps,reps,phi1,Phi,phi2 - integer(pInt), dimension(7) :: pos + integer(pInt), dimension(7) :: myPos integer(pInt), dimension(3) :: steps integer(pInt), dimension(:), allocatable :: binSet real(pReal) center,sum_dV_V,prob,dg_0,C,lowerC,upperC,rnd @@ -429,18 +427,18 @@ end function !--- parse header of ODF file --- !--- limits in phi1, Phi, phi2 --- read(999,fmt=fileFormat,end=100) line - pos = IO_stringPos(line,3_pInt) - if (pos(1).ne.3) goto 100 - do i=1,3 - limits(i) = IO_floatValue(line,pos,i)*inRad + myPos = IO_stringPos(line,3_pInt) + if (myPos(1).ne.3) goto 100 + do i=1_pInt,3_pInt + limits(i) = IO_floatValue(line,myPos,i)*inRad enddo !--- deltas in phi1, Phi, phi2 --- read(999,fmt=fileFormat,end=100) line - pos = IO_stringPos(line,3_pInt) - if (pos(1).ne.3) goto 100 - do i=1,3 - deltas(i) = IO_floatValue(line,pos,i)*inRad + myPos = IO_stringPos(line,3_pInt) + if (myPos(1).ne.3) goto 100 + do i=1_pInt,3_pInt + deltas(i) = IO_floatValue(line,myPos,i)*inRad enddo steps = nint(limits/deltas,pInt) allocate(dV_V(steps(3),steps(2),steps(1))) @@ -461,12 +459,12 @@ end function dg_0 = deltas(1)*deltas(3)*2.0_pReal*sin(deltas(2)/2.0_pReal) NnonZero = 0_pInt - do phi1=1,steps(1) - do Phi=1,steps(2) - do phi2=1,steps(3) + do phi1=1_pInt,steps(1) + do Phi=1_pInt,steps(2) + do phi2=1_pInt,steps(3) read(999,fmt=*,end=100) prob if (prob > 0.0_pReal) then - NnonZero = NnonZero+1 + NnonZero = NnonZero+1_pInt sum_dV_V = sum_dV_V+prob else prob = 0.0_pReal @@ -506,19 +504,19 @@ end function allocate(binSet(Nreps)) bin = 0_pInt ! bin counter - i = 1 ! set counter - do phi1=1,steps(1) - do Phi=1,steps(2) - do phi2=1,steps(3) + i = 1_pInt ! set counter + do phi1=1_pInt,steps(1) + do Phi=1_pInt,steps(2) + do phi2=1_pInt,steps(3) reps = nint(C*dV_V(phi2,Phi,phi1), pInt) binSet(i:i+reps-1) = bin - bin = bin+1 ! advance bin + bin = bin+1_pInt ! advance bin i = i+reps ! advance set enddo enddo enddo - do i=1,Nast + do i=1_pInt,Nast if (i < Nast) then call random_number(rnd) j = nint(rnd*(Nreps-i)+i+0.5_pReal,pInt) @@ -552,7 +550,7 @@ end function character(len=*), intent(in) :: line character(len=*), parameter :: blank = achar(32)//achar(9)//achar(10)//achar(13) ! whitespaces character(len=*), parameter :: comment = achar(35) ! comment id '#' - integer(pInt) posNonBlank, posComment + integer :: posNonBlank, posComment !no pInt logical IO_isBlank posNonBlank = verify(line,blank) @@ -572,7 +570,7 @@ end function character(len=*), intent(in) :: line,openChar,closeChar character(len=*), parameter :: sep=achar(32)//achar(9)//achar(10)//achar(13) ! whitespaces character(len=len_trim(line)) IO_getTag - integer(pInt) left,right + integer :: left,right !no pInt IO_getTag = '' left = scan(line,openChar) @@ -596,7 +594,7 @@ end function integer(pInt) IO_countSections character(len=1024) line - IO_countSections = 0 + IO_countSections = 0_pInt line = '' rewind(file) @@ -609,7 +607,7 @@ end function if (IO_isBlank(line)) cycle ! skip empty lines if (IO_getTag(line,'<','>') /= '') exit ! stop at next part if (IO_getTag(line,'[',']') /= '') & ! found [section] identifier - IO_countSections = IO_countSections + 1 + IO_countSections = IO_countSections + 1_pInt enddo 100 endfunction @@ -627,7 +625,7 @@ end function integer(pInt), intent(in) :: file, Nsections character(len=*), intent(in) :: part, myTag integer(pInt), dimension(Nsections) :: IO_countTagInPart, counter - integer(pInt), parameter :: maxNchunks = 1 + integer(pInt), parameter :: maxNchunks = 1_pInt integer(pInt), dimension(1+2*maxNchunks) :: positions integer(pInt) section character(len=1024) line,tag @@ -646,7 +644,7 @@ end function if (IO_isBlank(line)) cycle ! skip empty lines if (IO_getTag(line,'<','>') /= '') exit ! stop at next part if (IO_getTag(line,'[',']') /= '') & ! found [section] identifier - section = section + 1 + section = section + 1_pInt if (section > 0) then positions = IO_stringPos(line,maxNchunks) tag = IO_lc(IO_stringValue(line,positions,1_pInt)) ! extract key @@ -672,7 +670,7 @@ endfunction integer(pInt), intent(in) :: file, Nsections character(len=*), intent(in) :: part, myTag logical, dimension(Nsections) :: IO_spotTagInPart - integer(pInt), parameter :: maxNchunks = 1 + integer(pInt), parameter :: maxNchunks = 1_pInt integer(pInt), dimension(1+2*maxNchunks) :: positions integer(pInt) section character(len=1024) line,tag @@ -691,8 +689,8 @@ endfunction if (IO_isBlank(line)) cycle ! skip empty lines if (IO_getTag(line,'<','>') /= '') exit ! stop at next part if (IO_getTag(line,'[',']') /= '') & ! found [section] identifier - section = section + 1 - if (section > 0) then + section = section + 1_pInt + if (section > 0_pInt) then positions = IO_stringPos(line,maxNchunks) tag = IO_lc(IO_stringValue(line,positions,1_pInt)) ! extract key if (tag == myTag) & ! match @@ -716,11 +714,11 @@ endfunction character(len=*), intent(in) :: line character(len=*), parameter :: sep=achar(44)//achar(32)//achar(9)//achar(10)//achar(13) ! comma and whitespaces integer(pInt), intent(in) :: N - integer(pInt) left,right - integer(pInt) IO_stringPos(1+N*2) + integer :: left,right !no pInt (verify and scan return default integer) + integer(pInt) :: IO_stringPos(1_pInt+N*2_pInt) - IO_stringPos = -1 - IO_stringPos(1) = 0 + IO_stringPos = -1_pInt + IO_stringPos(1) = 0_pInt right = 0 do while (verify(line(right+1:),sep)>0) @@ -730,55 +728,55 @@ endfunction exit endif if ( IO_stringPos(1) 1) then - read(UNIT=line(ends(pos)+1:ends(pos)+pos_exp-1),ERR=100,FMT=*) base - read(UNIT=line(ends(pos)+pos_exp:ends(pos+1)),ERR=100,FMT=*) expon + read(UNIT=line(ends(myPos)+1:ends(myPos)+pos_exp-1),ERR=100,FMT=*) base + read(UNIT=line(ends(myPos)+pos_exp:ends(myPos+1)),ERR=100,FMT=*) expon else - read(UNIT=line(ends(pos)+1:ends(pos+1)),ERR=100,FMT=*) base + read(UNIT=line(ends(myPos)+1:ends(myPos+1)),ERR=100,FMT=*) base expon = 0_pInt endif IO_fixedNoEFloatValue = base*10.0_pReal**expon @@ -848,21 +847,21 @@ endfunction !******************************************************************** -! read int value at pos from line +! read int value at myPos from line !******************************************************************** - pure function IO_intValue (line,positions,pos) + pure function IO_intValue (line,positions,myPos) use prec, only: pReal,pInt implicit none character(len=*), intent(in) :: line - integer(pInt), intent(in) :: positions(*),pos - integer(pInt) IO_intValue + integer(pInt), intent(in) :: positions(*),myPos + integer(pInt) :: IO_intValue - if (positions(1) < pos) then + if (positions(1) < myPos) then IO_intValue = 0_pInt else - read(UNIT=line(positions(pos*2):positions(pos*2+1)),ERR=100,FMT=*) IO_intValue + read(UNIT=line(positions(myPos*2):positions(myPos*2+1)),ERR=100,FMT=*) IO_intValue endif return 100 IO_intValue = huge(1_pInt) @@ -871,18 +870,18 @@ endfunction !******************************************************************** -! read int value at pos from fixed format line +! read int value at myPos from fixed format line !******************************************************************** - pure function IO_fixedIntValue (line,ends,pos) + pure function IO_fixedIntValue (line,ends,myPos) use prec, only: pReal,pInt implicit none character(len=*), intent(in) :: line - integer(pInt), intent(in) :: ends(*),pos + integer(pInt), intent(in) :: ends(*),myPos integer(pInt) IO_fixedIntValue - read(UNIT=line(ends(pos)+1:ends(pos+1)),ERR=100,FMT=*) IO_fixedIntValue + read(UNIT=line(ends(myPos)+1:ends(myPos+1)),ERR=100,FMT=*) IO_fixedIntValue return 100 IO_fixedIntValue = huge(1_pInt) @@ -899,7 +898,7 @@ endfunction character (len=*), intent(in) :: line character (len=len(line)) IO_lc - integer(pInt) i + integer :: i !no pInt (len returns default integer) IO_lc = line do i=1,len(line) @@ -919,7 +918,7 @@ endfunction character (len=*) line character (len=len(line)) IO_lc - integer(pInt) i + integer :: i !no pInt (len returns default integer) IO_lc = line do i=1,len(line) @@ -939,15 +938,15 @@ endfunction implicit none integer(pInt) remainingChunks,unit,N - integer(pInt), parameter :: maxNchunks = 64 - integer(pInt), dimension(1+2*maxNchunks) :: pos + integer(pInt), parameter :: maxNchunks = 64_pInt + integer(pInt), dimension(1+2*maxNchunks) :: myPos character(len=300) line remainingChunks = N do while (remainingChunks > 0) read(unit,'(A300)',end=100) line - pos = IO_stringPos(line,maxNchunks) - remainingChunks = remainingChunks - pos(1) + myPos = IO_stringPos(line,maxNchunks) + remainingChunks = remainingChunks - myPos(1) enddo 100 endsubroutine @@ -957,19 +956,18 @@ endfunction !******************************************************************** pure function IO_extractValue (line,key) - use prec, only: pReal,pInt implicit none character(len=*), intent(in) :: line,key character(len=*), parameter :: sep = achar(61) ! '=' - integer(pInt) pos + integer :: myPos ! no pInt (scan returns default integer) character(len=300) IO_extractValue IO_extractValue = '' - pos = scan(line,sep) - if (pos > 0 .and. line(:pos-1) == key(:pos-1)) & ! key matches expected key - IO_extractValue = line(pos+1:) ! extract value + myPos = scan(line,sep) + if (myPos > 0 .and. line(:myPos-1) == key(:myPos-1)) & ! key matches expected key + IO_extractValue = line(myPos+1:) ! extract value endfunction @@ -980,28 +978,28 @@ endfunction ! : is not changed back to the original version since *.inp_assembly does not ! : contain any comment lines (12.07.2010) !******************************************************************** - function IO_countDataLines (unit) + function IO_countDataLines (myUnit) use prec, only: pReal,pInt implicit none - integer(pInt) IO_countDataLines,unit - integer(pInt), parameter :: maxNchunks = 1 - integer(pInt), dimension(1+2*maxNchunks) :: pos - character(len=300) line,tmp + integer(pInt) :: IO_countDataLines,myUnit + integer(pInt), parameter :: maxNchunks = 1_pInt + integer(pInt), dimension(1+2*maxNchunks) :: myPos + character(len=300) :: line,tmp - IO_countDataLines = 0 + IO_countDataLines = 0_pInt do - read(unit,'(A300)',end=100) line - pos = IO_stringPos(line,maxNchunks) - tmp = IO_lc(IO_stringValue(line,pos,1_pInt)) + read(myUnit,'(A300)',end=100) line + myPos = IO_stringPos(line,maxNchunks) + tmp = IO_lc(IO_stringValue(line,myPos,1_pInt)) if (tmp(1:1) == '*' .and. tmp(2:2) /= '*') then ! found keyword exit else if (tmp(2:2) /= '*') IO_countDataLines = IO_countDataLines + 1_pInt endif enddo -100 backspace(unit) +100 backspace(myUnit) endfunction @@ -1011,16 +1009,16 @@ endfunction ! Marc: ints concatenated by "c" as last char or range of values a "to" b ! Abaqus: triplet of start,stop,inc !******************************************************************** - function IO_countContinousIntValues (unit) + function IO_countContinousIntValues (myUnit) use DAMASK_interface use prec, only: pReal,pInt implicit none - integer(pInt) unit,l,count - integer(pInt) IO_countContinousIntValues - integer(pInt), parameter :: maxNchunks = 8192 - integer(pInt), dimension(1+2*maxNchunks) :: pos + integer(pInt) :: myUnit,l,c + integer(pInt) :: IO_countContinousIntValues + integer(pInt), parameter :: maxNchunks = 8192_pInt + integer(pInt), dimension(1+2*maxNchunks) :: myPos character(len=65536) line IO_countContinousIntValues = 0_pInt @@ -1029,15 +1027,15 @@ endfunction case ('Marc','Spectral') do - read(unit,'(A300)',end=100) line - pos = IO_stringPos(line,maxNchunks) - if (IO_lc(IO_stringValue(line,pos,2_pInt)) == 'to' ) then ! found range indicator - IO_countContinousIntValues = 1_pInt + IO_intValue(line,pos,3_pInt) - IO_intValue(line,pos,1_pInt) + read(myUnit,'(A300)',end=100) line + myPos = IO_stringPos(line,maxNchunks) + if (IO_lc(IO_stringValue(line,myPos,2_pInt)) == 'to' ) then ! found range indicator + IO_countContinousIntValues = 1_pInt + IO_intValue(line,myPos,3_pInt) - IO_intValue(line,myPos,1_pInt) exit ! only one single range indicator allowed else - IO_countContinousIntValues = IO_countContinousIntValues+pos(1)-1_pInt ! add line's count when assuming 'c' - if ( IO_lc(IO_stringValue(line,pos,pos(1))) /= 'c' ) then ! line finished, read last value - IO_countContinousIntValues = IO_countContinousIntValues+1 + IO_countContinousIntValues = IO_countContinousIntValues+myPos(1)-1_pInt ! add line's count when assuming 'c' + if ( IO_lc(IO_stringValue(line,myPos,myPos(1))) /= 'c' ) then ! line finished, read last value + IO_countContinousIntValues = IO_countContinousIntValues+1_pInt exit ! data ended endif endif @@ -1045,17 +1043,17 @@ endfunction case('Abaqus') - count = IO_countDataLines(unit) - do l = 1,count - backspace(unit) + c = IO_countDataLines(myUnit) + do l = 1_pInt,c + backspace(myUnit) enddo - do l = 1,count - read(unit,'(A300)',end=100) line - pos = IO_stringPos(line,maxNchunks) - IO_countContinousIntValues = IO_countContinousIntValues + 1 + & ! assuming range generation - (IO_intValue(line,pos,2_pInt)-IO_intValue(line,pos,1_pInt))/& - max(1_pInt,IO_intValue(line,pos,3_pInt)) + do l = 1_pInt,c + read(myUnit,'(A300)',end=100) line + myPos = IO_stringPos(line,maxNchunks) + IO_countContinousIntValues = IO_countContinousIntValues + 1_pInt + & ! assuming range generation + (IO_intValue(line,myPos,2_pInt)-IO_intValue(line,myPos,1_pInt))/& + max(1_pInt,IO_intValue(line,myPos,3_pInt)) enddo endselect @@ -1068,53 +1066,53 @@ endfunction ! Marc: ints concatenated by "c" as last char, range of a "to" b, or named set ! Abaqus: triplet of start,stop,inc or named set !******************************************************************** - function IO_continousIntValues (unit,maxN,lookupName,lookupMap,lookupMaxN) + function IO_continousIntValues (myUnit,maxN,lookupName,lookupMap,lookupMaxN) use DAMASK_interface use prec, only: pReal,pInt implicit none - integer(pInt) unit,maxN,i,j,l,count,first,last + integer(pInt) myUnit,maxN,i,j,l,c,first,last integer(pInt), dimension(1+maxN) :: IO_continousIntValues integer(pInt), parameter :: maxNchunks = 8192_pInt - integer(pInt), dimension(1+2*maxNchunks) :: pos + integer(pInt), dimension(1+2*maxNchunks) :: myPos character(len=64), dimension(:) :: lookupName integer(pInt) :: lookupMaxN integer(pInt), dimension(:,:) :: lookupMap character(len=65536) line logical rangeGeneration - IO_continousIntValues = 0 + IO_continousIntValues = 0_pInt rangeGeneration = .false. select case (FEsolver) case ('Marc','Spectral') do - read(unit,'(A65536)',end=100) line - pos = IO_stringPos(line,maxNchunks) - if (verify(IO_stringValue(line,pos,1_pInt),'0123456789') > 0) then ! a non-int, i.e. set name - do i = 1_pInt,lookupMaxN ! loop over known set names - if (IO_stringValue(line,pos,1_pInt) == lookupName(i)) then ! found matching name + read(myUnit,'(A65536)',end=100) line + myPos = IO_stringPos(line,maxNchunks) + if (verify(IO_stringValue(line,myPos,1_pInt),'0123456789') > 0) then ! a non-int, i.e. set name + do i = 1_pInt, lookupMaxN ! loop over known set names + if (IO_stringValue(line,myPos,1_pInt) == lookupName(i)) then ! found matching name IO_continousIntValues = lookupMap(:,i) ! return resp. entity list exit endif enddo exit - else if (pos(1) > 2_pInt .and. IO_lc(IO_stringValue(line,pos,2_pInt)) == 'to' ) then ! found range indicator - do i = IO_intValue(line,pos,1_pInt),IO_intValue(line,pos,3_pInt) + else if (myPos(1) > 2_pInt .and. IO_lc(IO_stringValue(line,myPos,2_pInt)) == 'to' ) then ! found range indicator + do i = IO_intValue(line,myPos,1_pInt),IO_intValue(line,myPos,3_pInt) IO_continousIntValues(1) = IO_continousIntValues(1) + 1_pInt IO_continousIntValues(1+IO_continousIntValues(1)) = i enddo exit else - do i = 1,pos(1)-1 ! interpret up to second to last value + do i = 1_pInt ,myPos(1)-1_pInt ! interpret up to second to last value IO_continousIntValues(1) = IO_continousIntValues(1) + 1_pInt - IO_continousIntValues(1+IO_continousIntValues(1)) = IO_intValue(line,pos,i) + IO_continousIntValues(1+IO_continousIntValues(1)) = IO_intValue(line,myPos,i) enddo - if ( IO_lc(IO_stringValue(line,pos,pos(1))) /= 'c' ) then ! line finished, read last value + if ( IO_lc(IO_stringValue(line,myPos,myPos(1))) /= 'c' ) then ! line finished, read last value IO_continousIntValues(1) = IO_continousIntValues(1) + 1_pInt - IO_continousIntValues(1+IO_continousIntValues(1)) = IO_intValue(line,pos,pos(1)) + IO_continousIntValues(1+IO_continousIntValues(1)) = IO_intValue(line,myPos,myPos(1)) exit endif endif @@ -1122,43 +1120,43 @@ endfunction case('Abaqus') - count = IO_countDataLines(unit) - do l = 1,count - backspace(unit) + c = IO_countDataLines(myUnit) + do l = 1_pInt,c + backspace(myUnit) enddo ! check if the element values in the elset are auto generated - backspace(unit) - read(unit,'(A65536)',end=100) line - pos = IO_stringPos(line,maxNchunks) - do i = 1,pos(1) - if (IO_lc(IO_stringValue(line,pos,i)) == 'generate') rangeGeneration = .true. + backspace(myUnit) + read(myUnit,'(A65536)',end=100) line + myPos = IO_stringPos(line,maxNchunks) + do i = 1_pInt,myPos(1) + if (IO_lc(IO_stringValue(line,myPos,i)) == 'generate') rangeGeneration = .true. enddo - do l = 1,count - read(unit,'(A65536)',end=100) line - pos = IO_stringPos(line,maxNchunks) - if (verify(IO_stringValue(line,pos,1_pInt),'0123456789') > 0_pInt) then ! a non-int, i.e. set names follow on this line - do i = 1,pos(1) ! loop over set names in line - do j = 1,lookupMaxN ! look thru known set names - if (IO_stringValue(line,pos,i) == lookupName(j)) then ! found matching name - first = 2 + IO_continousIntValues(1) ! where to start appending data - last = first + lookupMap(1,j) - 1 ! up to where to append data + do l = 1_pInt,c + read(myUnit,'(A65536)',end=100) line + myPos = IO_stringPos(line,maxNchunks) + if (verify(IO_stringValue(line,myPos,1_pInt),'0123456789') > 0) then ! a non-int, i.e. set names follow on this line + do i = 1_pInt,myPos(1) ! loop over set names in line + do j = 1_pInt,lookupMaxN ! look thru known set names + if (IO_stringValue(line,myPos,i) == lookupName(j)) then ! found matching name + first = 2_pInt + IO_continousIntValues(1) ! where to start appending data + last = first + lookupMap(1,j) - 1_pInt ! up to where to append data IO_continousIntValues(first:last) = lookupMap(2:1+lookupMap(1,j),j) ! add resp. entity list IO_continousIntValues(1) = IO_continousIntValues(1) + lookupMap(1,j) ! count them endif enddo enddo else if (rangeGeneration) then ! range generation - do i = IO_intValue(line,pos,1_pInt),IO_intValue(line,pos,2_pInt),max(1_pInt,IO_intValue(line,pos,3_pInt)) - IO_continousIntValues(1) = IO_continousIntValues(1) + 1 + do i = IO_intValue(line,myPos,1_pInt),IO_intValue(line,myPos,2_pInt),max(1_pInt,IO_intValue(line,myPos,3_pInt)) + IO_continousIntValues(1) = IO_continousIntValues(1) + 1_pInt IO_continousIntValues(1+IO_continousIntValues(1)) = i enddo else ! read individual elem nums - do i = 1,pos(1) -! write(*,*)'IO_CIV-int',IO_intValue(line,pos,i) - IO_continousIntValues(1) = IO_continousIntValues(1) + 1 - IO_continousIntValues(1+IO_continousIntValues(1)) = IO_intValue(line,pos,i) + do i = 1_pInt,myPos(1) +! write(*,*)'IO_CIV-int',IO_intValue(line,myPos,i) + IO_continousIntValues(1) = IO_continousIntValues(1) + 1_pInt + IO_continousIntValues(1+IO_continousIntValues(1)) = IO_intValue(line,myPos,i) enddo endif enddo diff --git a/code/crystallite.f90 b/code/crystallite.f90 index f5ee810e9..4956f349a 100644 --- a/code/crystallite.f90 +++ b/code/crystallite.f90 @@ -3178,7 +3178,7 @@ function crystallite_postResults(& crystallite_postResults = 0.0_pReal c = 0_pInt - crystallite_postResults(c+1) = crystallite_sizePostResults(crystID) ! size of results from cryst + crystallite_postResults(c+1) = real(crystallite_sizePostResults(crystID),pReal) ! size of results from cryst c = c + 1_pInt do o = 1,crystallite_Noutput(crystID) @@ -3186,10 +3186,10 @@ function crystallite_postResults(& select case(crystallite_output(o,crystID)) case ('phase') mySize = 1_pInt - crystallite_postResults(c+1) = material_phase(g,i,e) ! phaseID of grain + crystallite_postResults(c+1) = real(material_phase(g,i,e),pReal) ! phaseID of grain case ('texture') mySize = 1_pInt - crystallite_postResults(c+1) = material_texture(g,i,e) ! textureID of grain + crystallite_postResults(c+1) = real(material_texture(g,i,e),pReal) ! textureID of grain case ('volume') mySize = 1_pInt detF = math_det33(crystallite_partionedF(1:3,1:3,g,i,e)) ! V_current = det(F) * V_reference @@ -3210,36 +3210,36 @@ function crystallite_postResults(& case ('defgrad','f') mySize = 9_pInt - crystallite_postResults(c+1:c+mySize) = reshape(math_transpose33(crystallite_partionedF(1:3,1:3,g,i,e)),(/mySize/)) + crystallite_postResults(c+1:c+mySize) = reshape(math_transpose33(crystallite_partionedF(1:3,1:3,g,i,e)),[mySize]) case ('e') mySize = 9_pInt crystallite_postResults(c+1:c+mySize) = 0.5_pReal * reshape((math_mul33x33( & math_transpose33(crystallite_partionedF(1:3,1:3,g,i,e)), & - crystallite_partionedF(1:3,1:3,g,i,e)) - math_I3),(/mySize/)) + crystallite_partionedF(1:3,1:3,g,i,e)) - math_I3),[mySize]) case ('fe') mySize = 9_pInt - crystallite_postResults(c+1:c+mySize) = reshape(math_transpose33(crystallite_Fe(1:3,1:3,g,i,e)),(/mySize/)) + crystallite_postResults(c+1:c+mySize) = reshape(math_transpose33(crystallite_Fe(1:3,1:3,g,i,e)),[mySize]) case ('ee') Ee = 0.5_pReal * (math_mul33x33(math_transpose33(crystallite_Fe(1:3,1:3,g,i,e)), crystallite_Fe(1:3,1:3,g,i,e)) - math_I3) mySize = 9_pInt - crystallite_postResults(c+1:c+mySize) = reshape(Ee,(/mySize/)) + crystallite_postResults(c+1:c+mySize) = reshape(Ee,[mySize]) case ('fp') mySize = 9_pInt - crystallite_postResults(c+1:c+mySize) = reshape(math_transpose33(crystallite_Fp(1:3,1:3,g,i,e)),(/mySize/)) + crystallite_postResults(c+1:c+mySize) = reshape(math_transpose33(crystallite_Fp(1:3,1:3,g,i,e)),[mySize]) case ('lp') mySize = 9_pInt - crystallite_postResults(c+1:c+mySize) = reshape(math_transpose33(crystallite_Lp(1:3,1:3,g,i,e)),(/mySize/)) + crystallite_postResults(c+1:c+mySize) = reshape(math_transpose33(crystallite_Lp(1:3,1:3,g,i,e)),[mySize]) case ('p','firstpiola','1stpiola') mySize = 9_pInt - crystallite_postResults(c+1:c+mySize) = reshape(math_transpose33(crystallite_P(1:3,1:3,g,i,e)),(/mySize/)) + crystallite_postResults(c+1:c+mySize) = reshape(math_transpose33(crystallite_P(1:3,1:3,g,i,e)),[mySize]) case ('s','tstar','secondpiola','2ndpiola') mySize = 9_pInt - crystallite_postResults(c+1:c+mySize) = reshape(math_Mandel6to33(crystallite_Tstar_v(1:6,g,i,e)),(/mySize/)) + crystallite_postResults(c+1:c+mySize) = reshape(math_Mandel6to33(crystallite_Tstar_v(1:6,g,i,e)),[mySize]) end select c = c + mySize enddo - crystallite_postResults(c+1) = constitutive_sizePostResults(g,i,e) ! size of constitutive results + crystallite_postResults(c+1) = real(constitutive_sizePostResults(g,i,e),pReal) ! size of constitutive results c = c + 1_pInt crystallite_postResults(c+1:c+constitutive_sizePostResults(g,i,e)) = constitutive_postResults(crystallite_Tstar_v(1:6,g,i,e), & crystallite_Fe, & diff --git a/code/homogenization.f90 b/code/homogenization.f90 index d9455de1d..d4ec99bad 100644 --- a/code/homogenization.f90 +++ b/code/homogenization.f90 @@ -562,7 +562,7 @@ subroutine materialpoint_postResults(dt) thePos = 0_pInt theSize = homogenization_sizePostResults(i,e) - materialpoint_results(thePos+1,i,e) = theSize ! tell size of homogenization results + materialpoint_results(thePos+1,i,e) = real(theSize,pReal) ! tell size of homogenization results thePos = thePos + 1_pInt if (theSize > 0_pInt) then ! any homogenization results to mention? @@ -570,7 +570,7 @@ subroutine materialpoint_postResults(dt) thePos = thePos + theSize endif - materialpoint_results(thePos+1,i,e) = myNgrains ! tell number of grains at materialpoint + materialpoint_results(thePos+1,i,e) = real(myNgrains,pReal) ! tell number of grains at materialpoint thePos = thePos + 1_pInt do g = 1,myNgrains ! loop over all grains diff --git a/code/homogenization_RGC.f90 b/code/homogenization_RGC.f90 index e54fbdf05..087e4b347 100644 --- a/code/homogenization_RGC.f90 +++ b/code/homogenization_RGC.f90 @@ -281,7 +281,7 @@ subroutine homogenization_RGC_partitionDeformation(& write(6,'(1x,a,i3,a,i3,a)')'========== Increment: ',theInc,' Cycle: ',cycleCounter,' ==========' write(6,'(1x,a32)')'Overall deformation gradient: ' do i = 1_pInt,3_pInt - write(6,'(1x,3(e14.8,1x))')(avgF(i,j), j = 1_pInt,3_pInt) + write(6,'(1x,3(e15.8,1x))')(avgF(i,j), j = 1_pInt,3_pInt) enddo write(6,*)' ' call flush(6) @@ -307,7 +307,7 @@ subroutine homogenization_RGC_partitionDeformation(& !$OMP CRITICAL (write2out) write(6,'(1x,a32,1x,i3)')'Deformation gradient of grain: ',iGrain do i = 1_pInt,3_pInt - write(6,'(1x,3(e14.8,1x))')(F(i,j,iGrain), j = 1_pInt,3_pInt) + write(6,'(1x,3(e15.8,1x))')(F(i,j,iGrain), j = 1_pInt,3_pInt) enddo write(6,*)' ' call flush(6) @@ -394,7 +394,7 @@ function homogenization_RGC_updateState(& !$OMP CRITICAL (write2out) write(6,'(1x,a30)')'Obtained state: ' do i = 1,3*nIntFaceTot - write(6,'(1x,2(e14.8,1x))')state%p(i) + write(6,'(1x,2(e15.8,1x))')state%p(i) enddo write(6,*)' ' !$OMP END CRITICAL (write2out) @@ -410,11 +410,11 @@ function homogenization_RGC_updateState(& if (debug_verbosity == 4) then !$OMP CRITICAL (write2out) do iGrain = 1,nGrain - write(6,'(1x,a30,1x,i3,1x,a4,3(1x,e14.8))')'Mismatch magnitude of grain(',iGrain,') :',NN(1,iGrain),NN(2,iGrain),NN(3,iGrain) + write(6,'(1x,a30,1x,i3,1x,a4,3(1x,e15.8))')'Mismatch magnitude of grain(',iGrain,') :',NN(1,iGrain),NN(2,iGrain),NN(3,iGrain) write(6,*)' ' write(6,'(1x,a30,1x,i3)')'Stress and penalties of grain: ',iGrain do i = 1,3 - write(6,'(1x,3(e14.8,1x),1x,3(e14.8,1x),1x,3(e14.8,1x))')(P(i,j,iGrain), j = 1,3), & + write(6,'(1x,3(e15.8,1x),1x,3(e15.8,1x),1x,3(e15.8,1x))')(P(i,j,iGrain), j = 1,3), & (R(i,j,iGrain), j = 1,3), & (D(i,j,iGrain), j = 1,3) enddo @@ -459,7 +459,7 @@ function homogenization_RGC_updateState(& if (debug_verbosity == 4) then !$OMP CRITICAL (write2out) write(6,'(1x,a30,1x,i3)')'Traction at interface: ',iNum - write(6,'(1x,3(e14.8,1x))')(tract(iNum,j), j = 1,3) + write(6,'(1x,3(e15.8,1x))')(tract(iNum,j), j = 1,3) write(6,*)' ' !$OMP END CRITICAL (write2out) endif @@ -478,9 +478,9 @@ function homogenization_RGC_updateState(& !$OMP CRITICAL (write2out) write(6,'(1x,a)')' ' write(6,'(1x,a,1x,i2,1x,i4)')'RGC residual check ...',ip,el - write(6,'(1x,a15,1x,e14.8,1x,a7,i3,1x,a12,i2,i2)')'Max stress: ',stresMax, & + write(6,'(1x,a15,1x,e15.8,1x,a7,i3,1x,a12,i2,i2)')'Max stress: ',stresMax, & '@ grain',stresLoc(3),'in component',stresLoc(1),stresLoc(2) - write(6,'(1x,a15,1x,e14.8,1x,a7,i3,1x,a12,i2)')'Max residual: ',residMax, & + write(6,'(1x,a15,1x,e15.8,1x,a7,i3,1x,a12,i2)')'Max residual: ',residMax, & '@ iface',residLoc(1),'in direction',residLoc(2) call flush(6) !$OMP END CRITICAL (write2out) @@ -523,15 +523,15 @@ function homogenization_RGC_updateState(& if (debug_verbosity == 4 .and. debug_e == el .and. debug_i == ip) then !$OMP CRITICAL (write2out) - write(6,'(1x,a30,1x,e14.8)')'Constitutive work: ',constitutiveWork - write(6,'(1x,a30,3(1x,e14.8))')'Magnitude mismatch: ',sum(NN(1,:))/dble(nGrain), & + write(6,'(1x,a30,1x,e15.8)')'Constitutive work: ',constitutiveWork + write(6,'(1x,a30,3(1x,e15.8))')'Magnitude mismatch: ',sum(NN(1,:))/dble(nGrain), & sum(NN(2,:))/dble(nGrain), & sum(NN(3,:))/dble(nGrain) - write(6,'(1x,a30,1x,e14.8)')'Penalty energy: ',penaltyEnergy - write(6,'(1x,a30,1x,e14.8)')'Volume discrepancy: ',volDiscrep + write(6,'(1x,a30,1x,e15.8)')'Penalty energy: ',penaltyEnergy + write(6,'(1x,a30,1x,e15.8)')'Volume discrepancy: ',volDiscrep write(6,*)'' - write(6,'(1x,a30,1x,e14.8)')'Maximum relaxation rate: ',maxval(abs(drelax))/dt - write(6,'(1x,a30,1x,e14.8)')'Average relaxation rate: ',sum(abs(drelax))/dt/dble(3*nIntFaceTot) + write(6,'(1x,a30,1x,e15.8)')'Maximum relaxation rate: ',maxval(abs(drelax))/dt + write(6,'(1x,a30,1x,e15.8)')'Average relaxation rate: ',sum(abs(drelax))/dt/dble(3*nIntFaceTot) write(6,*)'' call flush(6) !$OMP END CRITICAL (write2out) @@ -619,7 +619,7 @@ function homogenization_RGC_updateState(& !$OMP CRITICAL (write2out) write(6,'(1x,a30)')'Jacobian matrix of stress' do i = 1,3*nIntFaceTot - write(6,'(1x,100(e10.4,1x))')(smatrix(i,j), j = 1,3*nIntFaceTot) + write(6,'(1x,100(e11.4,1x))')(smatrix(i,j), j = 1,3*nIntFaceTot) enddo write(6,*)' ' call flush(6) @@ -675,7 +675,7 @@ function homogenization_RGC_updateState(& !$OMP CRITICAL (write2out) write(6,'(1x,a30)')'Jacobian matrix of penalty' do i = 1,3*nIntFaceTot - write(6,'(1x,100(e10.4,1x))')(pmatrix(i,j), j = 1,3*nIntFaceTot) + write(6,'(1x,100(e11.4,1x))')(pmatrix(i,j), j = 1,3*nIntFaceTot) enddo write(6,*)' ' call flush(6) @@ -695,7 +695,7 @@ function homogenization_RGC_updateState(& !$OMP CRITICAL (write2out) write(6,'(1x,a30)')'Jacobian matrix of penalty' do i = 1,3*nIntFaceTot - write(6,'(1x,100(e10.4,1x))')(rmatrix(i,j), j = 1,3*nIntFaceTot) + write(6,'(1x,100(e11.4,1x))')(rmatrix(i,j), j = 1,3*nIntFaceTot) enddo write(6,*)' ' call flush(6) @@ -709,7 +709,7 @@ function homogenization_RGC_updateState(& !$OMP CRITICAL (write2out) write(6,'(1x,a30)')'Jacobian matrix (total)' do i = 1,3*nIntFaceTot - write(6,'(1x,100(e10.4,1x))')(jmatrix(i,j), j = 1,3*nIntFaceTot) + write(6,'(1x,100(e11.4,1x))')(jmatrix(i,j), j = 1,3*nIntFaceTot) enddo write(6,*)' ' call flush(6) @@ -728,7 +728,7 @@ function homogenization_RGC_updateState(& !$OMP CRITICAL (write2out) write(6,'(1x,a30)')'Jacobian inverse' do i = 1,3*nIntFaceTot - write(6,'(1x,100(e10.4,1x))')(jnverse(i,j), j = 1_pInt,3_pInt*nIntFaceTot) + write(6,'(1x,100(e11.4,1x))')(jnverse(i,j), j = 1_pInt,3_pInt*nIntFaceTot) enddo write(6,*)' ' call flush(6) @@ -748,7 +748,7 @@ function homogenization_RGC_updateState(& homogenization_RGC_updateState = (/.true.,.false./) !$OMP CRITICAL (write2out) write(6,'(1x,a,1x,i3,1x,a,1x,i3,1x,a)')'RGC_updateState: ip',ip,'| el',el,'enforces cutback' - write(6,'(1x,a,1x,e14.8)')'due to large relaxation change =',maxval(abs(drelax)) + write(6,'(1x,a,1x,e15.8)')'due to large relaxation change =',maxval(abs(drelax)) call flush(6) !$OMP END CRITICAL (write2out) endif @@ -758,7 +758,7 @@ function homogenization_RGC_updateState(& !$OMP CRITICAL (write2out) write(6,'(1x,a30)')'Returned state: ' do i = 1,3*nIntFaceTot - write(6,'(1x,2(e14.8,1x))')state%p(i) + write(6,'(1x,2(e15.8,1x))')state%p(i) enddo write(6,*)' ' call flush(6) @@ -810,7 +810,7 @@ subroutine homogenization_RGC_averageStressAndItsTangent(& dPdF99 = math_Plain3333to99(dPdF(1:3,1:3,1:3,1:3,iGrain)) write(6,'(1x,a30,1x,i3)')'Stress tangent of grain: ',iGrain do i = 1,9 - write(6,'(1x,(e14.8,1x))') (dPdF99(i,j), j = 1,9) + write(6,'(1x,(e15.8,1x))') (dPdF99(i,j), j = 1,9) enddo write(6,*)' ' enddo @@ -955,7 +955,7 @@ subroutine homogenization_RGC_stressPenalty(& !* Debugging the surface correction factor ! if (ip == 1 .and. el == 1) then ! write(6,'(1x,a20,2(1x,i3))')'Correction factor: ',ip,el -! write(6,'(1x,3(e10.4,1x))')(surfCorr(i), i = 1,3) +! write(6,'(1x,3(e11.4,1x))')(surfCorr(i), i = 1,3) ! endif !* ------------------------------------------------------------------------------------------------------------- @@ -1005,9 +1005,9 @@ subroutine homogenization_RGC_stressPenalty(& ! if (ip == 1 .and. el == 1) then ! write(6,'(1x,a20,i2,1x,a20,1x,i3)')'Mismatch to face: ',intFace(1),'neighbor grain: ',iGNghb ! do i = 1,3 -! write(6,'(1x,3(e10.4,1x))')(nDef(i,j), j = 1,3) +! write(6,'(1x,3(e11.4,1x))')(nDef(i,j), j = 1,3) ! enddo -! write(6,'(1x,a20,e10.4))')'with magnitude: ',nDefNorm +! write(6,'(1x,a20,e11.4))')'with magnitude: ',nDefNorm ! endif !* Compute the stress penalty of all interfaces @@ -1030,7 +1030,7 @@ subroutine homogenization_RGC_stressPenalty(& ! if (ip == 1 .and. el == 1) then ! write(6,'(1x,a20,i2)')'Penalty of grain: ',iGrain ! do i = 1,3 -! write(6,'(1x,3(e10.4,1x))')(rPen(i,j,iGrain), j = 1,3) +! write(6,'(1x,3(e11.4,1x))')(rPen(i,j,iGrain), j = 1,3) ! enddo ! endif @@ -1090,7 +1090,7 @@ subroutine homogenization_RGC_volumePenalty(& ! if (ip == 1 .and. el == 1) then ! write(6,'(1x,a30,i2)')'Volume penalty of grain: ',iGrain ! do i = 1,3 -! write(6,'(1x,3(e10.4,1x))')(vPen(i,j,iGrain), j = 1,3) +! write(6,'(1x,3(e11.4,1x))')(vPen(i,j,iGrain), j = 1,3) ! enddo ! endif @@ -1230,16 +1230,16 @@ function homogenization_RGC_interfaceNormal(& !* Get the normal of the interface, identified from the value of intFace(1) homogenization_RGC_interfaceNormal = 0.0_pReal - nPos = abs(intFace(1)) ! identify the position of the interface in global state array - homogenization_RGC_interfaceNormal(nPos) = intFace(1)/abs(intFace(1)) ! get the normal vector w.r.t. cluster axis + nPos = abs(intFace(1)) ! identify the position of the interface in global state array + homogenization_RGC_interfaceNormal(nPos) = real(intFace(1)/abs(intFace(1)),pReal) ! get the normal vector w.r.t. cluster axis homogenization_RGC_interfaceNormal = & math_mul33x3(homogenization_RGC_orientation(:,:,ip,el),homogenization_RGC_interfaceNormal) - ! map the normal vector into sample coordinate system (basis) + ! map the normal vector into sample coordinate system (basis) ! if (ip == 1 .and. el == 1) then ! write(6,'(1x,a32,3(1x,i3))')'Interface normal: ',intFace(1) -! write(6,'(1x,3(e14.8,1x))')(nVect(i), i = 1,3) +! write(6,'(1x,3(e15.8,1x))')(nVect(i), i = 1,3) ! write(6,*)' ' ! call flush(6) ! endif diff --git a/code/kdtree2.f90 b/code/kdtree2.f90 index 95f069880..30b303cee 100644 --- a/code/kdtree2.f90 +++ b/code/kdtree2.f90 @@ -64,7 +64,7 @@ module kdtree2_priority_queue_module ! There are heap_size active elements. ! Assumes the allocation is always sufficient. Will NOT increase it ! to match. - integer(pInt) :: heap_size = 0 + integer(pInt) :: heap_size = 0_pInt type(kdtree2_result), pointer :: elems(:) end type pq @@ -97,14 +97,14 @@ contains type(pq) :: res ! ! - integer(pInt) :: nalloc + integer :: nalloc nalloc = size(results_in,1) if (nalloc .lt. 1) then write (*,*) 'PQ_CREATE: error, input arrays must be allocated.' end if res%elems => results_in - res%heap_size = 0 + res%heap_size = 0_pInt return end function pq_create @@ -160,8 +160,8 @@ contains i = i_in bigloop: do - l = 2*i ! left(i) - r = l+1 ! right(i) + l = 2_pInt*i ! left(i) + r = l+1_pInt ! right(i) ! ! set 'largest' to the index of either i, l, r ! depending on whose priority is largest. @@ -296,11 +296,11 @@ bigloop: do ! write (*,*) 'PQ_INSERT: error, attempt made to insert element on full PQ' ! stop ! else - a%heap_size = a%heap_size + 1 + a%heap_size = a%heap_size + 1_pInt i = a%heap_size do while (i .gt. 1) - isparent = int(i/2) + isparent = int(i/2_pInt, pInt) !needed casting? parentdis = a%elems(isparent)%dis if (dis .gt. parentdis) then ! move what was in i's parent into i. @@ -339,13 +339,13 @@ bigloop: do e = a%elems(i) parent = i - child = 2*i + child = 2_pInt*i N = a%heap_size do while (child .le. N) if (child .lt. N) then if (a%elems(child)%dis .lt. a%elems(child+1)%dis) then - child = child+1 + child = child+1_pInt endif endif prichild = a%elems(child)%dis @@ -355,7 +355,7 @@ bigloop: do ! move child into parent. a%elems(parent) = a%elems(child) parent = child - child = 2*parent + child = 2_pInt*parent end if end do a%elems(parent) = e @@ -384,8 +384,8 @@ bigloop: do if (.true.) then N=a%heap_size if (N .ge. 1) then - parent =1 - child=2 + parent =1_pInt + child=2_pInt loop: do while (child .le. N) prichild = a%elems(child)%dis @@ -396,9 +396,9 @@ bigloop: do ! if (child .lt. N) then - prichildp1 = a%elems(child+1)%dis + prichildp1 = a%elems(child+1_pInt)%dis if (prichild .lt. prichildp1) then - child = child+1 + child = child+1_pInt prichild = prichildp1 endif endif @@ -411,7 +411,7 @@ bigloop: do ! move child into parent. a%elems(parent) = a%elems(child) parent = child - child = 2*parent + child = 2_pInt*parent end if end do loop a%elems(parent)%dis = dis @@ -449,7 +449,7 @@ bigloop: do ! swap the item to be deleted with the last element ! and shorten heap by one. a%elems(i) = a%elems(a%heap_size) - a%heap_size = a%heap_size - 1 + a%heap_size = a%heap_size - 1_pInt call heapify(a,i) @@ -469,10 +469,10 @@ module kdtree2_module ! This module is identical to 'kd_tree', except that the order ! of subscripts is reversed in the data file. ! In otherwords for an embedding of N D-dimensional vectors, the - ! data file is here, in natural Fortran order data(1:D, 1:N) + ! data file is here, in natural Fortran order myData(1:D, 1:N) ! because Fortran lays out columns first, ! - ! whereas conventionally (C-style) it is data(1:N,1:D) + ! whereas conventionally (C-style) it is myData(1:N,1:D) ! as in the original kd_tree module. ! !-------------DATA TYPE, CREATION, DELETION--------------------- @@ -498,7 +498,7 @@ module kdtree2_module !---------------------------------------------------------------- - integer(pInt), parameter :: bucket_size = 12 + integer(pInt), parameter :: bucket_size = 12_pInt ! The maximum number of points to keep in a terminal node. type interval @@ -525,7 +525,7 @@ module kdtree2_module type :: kdtree2 ! Global information about the tree, one per tree - integer(pInt) :: dimen=0, n=0 + integer(pInt) :: dimen=0_pInt, n=0_pInt ! dimensionality and total # of points real(pReal), pointer :: the_data(:,:) => null() ! pointer to the actual data array @@ -570,7 +570,7 @@ module kdtree2_module integer(pInt) :: dimen integer(pInt) :: nn, nfound real(pReal) :: ballsize - integer(pInt) :: centeridx=999, correltime=9999 + integer(pInt) :: centeridx=999_pInt, correltime=9999_pInt ! exclude points within 'correltime' of 'centeridx', iff centeridx >= 0 integer(pInt) :: nalloc ! how much allocated for results(:)? logical :: rearrange ! are the data rearranged or original? @@ -579,7 +579,7 @@ module kdtree2_module real(pReal), pointer :: qv(:) ! query vector type(kdtree2_result), pointer :: results(:) ! results type(pq) :: pq - real(pReal), pointer :: data(:,:) ! temp pointer to data + real(pReal), pointer :: myData(:,:) ! temp pointer to data integer(pInt), pointer :: ind(:) ! temp pointer to indexes end type tree_search_record @@ -590,7 +590,7 @@ module kdtree2_module contains - function kdtree2_create(input_data,dim,sort,rearrange) result (mr) + function kdtree2_create(input_data,myDim,sort,rearrange) result (mr) ! ! create the actual tree structure, given an input array of data. ! @@ -598,9 +598,9 @@ contains ! THIS IS THE REVERSE OF THE PREVIOUS VERSION OF THIS MODULE. ! The reason for it is cache friendliness, improving performance. ! - ! Optional arguments: If 'dim' is specified, then the tree - ! will only search the first 'dim' components - ! of input_data, otherwise, dim is inferred + ! Optional arguments: If 'myDim' is specified, then the tree + ! will only search the first 'myDim' components + ! of input_data, otherwise, myDim is inferred ! from SIZE(input_data,1). ! ! if sort .eqv. .true. then output results @@ -615,7 +615,7 @@ contains ! ! .. Function Return Cut_value .. type (kdtree2), pointer :: mr - integer(pInt), intent(in), optional :: dim + integer(pInt), intent(in), optional :: myDim logical, intent(in), optional :: sort logical, intent(in), optional :: rearrange ! .. @@ -628,19 +628,19 @@ contains mr%the_data => input_data ! pointer assignment - if (present(dim)) then - mr%dimen = dim + if (present(myDim)) then + mr%dimen = myDim else - mr%dimen = size(input_data,1) + mr%dimen = int(size(input_data,1), pInt) ! size returns default integer end if - mr%n = size(input_data,2) + mr%n = int(size(input_data,2), pInt) ! size returns default integer if (mr%dimen > mr%n) then ! unlikely to be correct write (*,*) 'KD_TREE_TRANS: likely user error.' write (*,*) 'KD_TREE_TRANS: You passed in matrix with D=',mr%dimen write (*,*) 'KD_TREE_TRANS: and N=',mr%n - write (*,*) 'KD_TREE_TRANS: note, that new format is data(1:D,1:N)' + write (*,*) 'KD_TREE_TRANS: note, that new format is myData(1:D,1:N)' write (*,*) 'KD_TREE_TRANS: with usually N >> D. If N =approx= D, then a k-d tree' write (*,*) 'KD_TREE_TRANS: is not an appropriate data structure.' stop @@ -662,7 +662,7 @@ contains if (mr%rearrange) then allocate(mr%rearranged_data(mr%dimen,mr%n)) - do i=1,mr%n + do i=1_pInt,mr%n mr%rearranged_data(:,i) = mr%the_data(:, & mr%ind(i)) enddo @@ -679,7 +679,7 @@ contains type(tree_node), pointer :: dummy => null() ! .. allocate (tp%ind(tp%n)) - forall (j=1:tp%n) + forall (j=1_pInt:tp%n) tp%ind(j) = j end forall tp%root => build_tree_for_range(tp,1_pInt,tp%n, dummy) @@ -729,10 +729,10 @@ contains ! ! always compute true bounding box for terminal nodes. ! - do i=1,dimen + do i=1_pInt,dimen call spread_in_coordinate(tp,i,l,u,res%box(i)) end do - res%cut_dim = 0 + res%cut_dim = 0_pInt res%cut_val = 0.0_pReal res%l = l res%u = u @@ -764,7 +764,7 @@ contains end do - c = maxloc(res%box(1:dimen)%upper-res%box(1:dimen)%lower,1) + c = int(maxloc(res%box(1:dimen)%upper-res%box(1:dimen)%lower,1), pInt) ! ! c is the identity of which coordinate has the greatest spread. ! @@ -867,11 +867,11 @@ contains do while (lb < rb) if ( v(c,ind(lb)) <= alpha ) then ! it is good where it is. - lb = lb+1 + lb = lb+1_pInt else ! swap it with rb. tmp = ind(lb); ind(lb) = ind(rb); ind(rb) = tmp - rb = rb-1 + rb = rb-1_pInt endif end do @@ -879,7 +879,7 @@ contains if (v(c,ind(lb)) <= alpha) then res = lb else - res = lb-1 + res = lb-1_pInt endif end function select_on_coordinate_value @@ -901,9 +901,9 @@ contains do while (l=k) u = m - 1 + if (m<=k) l = m + 1_pInt + if (m>=k) u = m - 1_pInt end do end subroutine select_on_coordinate @@ -944,8 +944,8 @@ contains ulocal = u - do i = l + 2, ulocal, 2 - lmin = v(c,ind(i-1)) + do i = l + 2_pInt, ulocal, 2_pInt + lmin = v(c,ind(i-1_pInt)) lmax = v(c,ind(i)) if (lmin>lmax) then t = lmin @@ -1022,9 +1022,9 @@ contains sr%ballsize = huge(1.0_pReal) sr%qv => qv sr%nn = nn - sr%nfound = 0 - sr%centeridx = -1 - sr%correltime = 0 + sr%nfound = 0_pInt + sr%centeridx = -1_pInt + sr%correltime = 0_pInt sr%overflow = .false. sr%results => results @@ -1034,9 +1034,9 @@ contains sr%ind => tp%ind sr%rearrange = tp%rearrange if (tp%rearrange) then - sr%Data => tp%rearranged_data + sr%myData => tp%rearranged_data else - sr%Data => tp%the_data + sr%myData => tp%the_data endif sr%dimen = tp%dimen @@ -1067,7 +1067,7 @@ contains sr%correltime = correltime sr%nn = nn - sr%nfound = 0 + sr%nfound = 0_pInt sr%dimen = tp%dimen sr%nalloc = nn @@ -1078,9 +1078,9 @@ contains sr%rearrange = tp%rearrange if (sr%rearrange) then - sr%Data => tp%rearranged_data + sr%myData => tp%rearranged_data else - sr%Data => tp%the_data + sr%myData => tp%the_data endif call validate_query_storage(nn) @@ -1116,10 +1116,10 @@ contains ! sr%qv => qv sr%ballsize = r2 - sr%nn = 0 ! flag for fixed ball search - sr%nfound = 0 - sr%centeridx = -1 - sr%correltime = 0 + sr%nn = 0_pInt ! flag for fixed ball search + sr%nfound = 0_pInt + sr%centeridx = -1_pInt + sr%correltime = 0_pInt sr%results => results @@ -1130,9 +1130,9 @@ contains sr%rearrange= tp%rearrange if (tp%rearrange) then - sr%Data => tp%rearranged_data + sr%myData => tp%rearranged_data else - sr%Data => tp%the_data + sr%myData => tp%the_data endif sr%dimen = tp%dimen @@ -1176,8 +1176,8 @@ contains allocate (sr%qv(tp%dimen)) sr%qv = tp%the_data(:,idxin) ! copy the vector sr%ballsize = r2 - sr%nn = 0 ! flag for fixed r search - sr%nfound = 0 + sr%nn = 0_pInt ! flag for fixed r search + sr%nfound = 0_pInt sr%centeridx = idxin sr%correltime = correltime @@ -1195,9 +1195,9 @@ contains sr%rearrange = tp%rearrange if (tp%rearrange) then - sr%Data => tp%rearranged_data + sr%myData => tp%rearranged_data else - sr%Data => tp%the_data + sr%myData => tp%the_data endif sr%rearrange = tp%rearrange sr%dimen = tp%dimen @@ -1236,21 +1236,21 @@ contains sr%qv => qv sr%ballsize = r2 - sr%nn = 0 ! flag for fixed r search - sr%nfound = 0 - sr%centeridx = -1 - sr%correltime = 0 + sr%nn = 0_pInt ! flag for fixed r search + sr%nfound = 0_pInt + sr%centeridx = -1_pInt + sr%correltime = 0_pInt nullify(sr%results) ! for some reason, FTN 95 chokes on '=> null()' - sr%nalloc = 0 ! we do not allocate any storage but that's OK + sr%nalloc = 0_pInt ! we do not allocate any storage but that's OK ! for counting. sr%ind => tp%ind sr%rearrange = tp%rearrange if (tp%rearrange) then - sr%Data => tp%rearranged_data + sr%myData => tp%rearranged_data else - sr%Data => tp%the_data + sr%myData => tp%the_data endif sr%dimen = tp%dimen @@ -1285,22 +1285,22 @@ contains sr%qv = tp%the_data(:,idxin) sr%ballsize = r2 - sr%nn = 0 ! flag for fixed r search - sr%nfound = 0 + sr%nn = 0_pInt ! flag for fixed r search + sr%nfound = 0_pInt sr%centeridx = idxin sr%correltime = correltime nullify(sr%results) - sr%nalloc = 0 ! we do not allocate any storage but that's OK + sr%nalloc = 0_pInt ! we do not allocate any storage but that's OK ! for counting. sr%ind => tp%ind sr%rearrange = tp%rearrange if (sr%rearrange) then - sr%Data => tp%rearranged_data + sr%myData => tp%rearranged_data else - sr%Data => tp%the_data + sr%myData => tp%the_data endif sr%dimen = tp%dimen @@ -1324,7 +1324,7 @@ contains ! integer(pInt), intent(in) :: n - if (size(sr%results,1) .lt. n) then + if (int(size(sr%results,1),pInt) .lt. n) then write (*,*) 'KD_TREE_TRANS: you did not provide enough storage for results(1:n)' stop return @@ -1408,7 +1408,7 @@ contains ! check will also be false. ! box => node%box(1:) - do i=1,sr%dimen + do i=1_pInt,sr%dimen if (i .ne. cut_dim) then dis = dis + dis2_from_bnd(qv(i),box(i)%lower,box(i)%upper) if (dis > ballsize) then @@ -1486,7 +1486,7 @@ contains ! real(pReal), pointer :: qv(:) integer(pInt), pointer :: ind(:) - real(pReal), pointer :: data(:,:) + real(pReal), pointer :: myData(:,:) ! integer(pInt) :: dimen, i, indexofi, k, centeridx, correltime real(pReal) :: ballsize, sd, newpri @@ -1505,7 +1505,7 @@ contains ballsize = sr%ballsize rearrange = sr%rearrange ind => sr%ind(1:) - data => sr%Data(1:,1:) + myData => sr%myData(1:,1:) centeridx = sr%centeridx correltime = sr%correltime @@ -1517,7 +1517,7 @@ contains if (rearrange) then sd = 0.0_pReal do k = 1_pInt,dimen - sd = sd + (data(k,i) - qv(k))**2.0_pReal + sd = sd + (myData(k,i) - qv(k))**2.0_pReal if (sd>ballsize) cycle mainloop end do indexofi = ind(i) ! only read it if we have not broken out @@ -1525,7 +1525,7 @@ contains indexofi = ind(i) sd = 0.0_pReal do k = 1_pInt,dimen - sd = sd + (data(k,indexofi) - qv(k))**2.0_pReal + sd = sd + (myData(k,indexofi) - qv(k))**2.0_pReal if (sd>ballsize) cycle mainloop end do endif @@ -1557,7 +1557,7 @@ contains ! ! add this point unconditionally to fill list. ! - sr%nfound = sr%nfound +1 + sr%nfound = sr%nfound +1_pInt newpri = pq_insert(pqp,sd,indexofi) if (sr%nfound .eq. sr%nn) ballsize = newpri ! we have just filled the working list. @@ -1592,7 +1592,7 @@ contains ! real(pReal), pointer :: qv(:) integer(pInt), pointer :: ind(:) - real(pReal), pointer :: data(:,:) + real(pReal), pointer :: myData(:,:) ! integer(pInt) :: nfound integer(pInt) :: dimen, i, indexofi, k @@ -1608,7 +1608,7 @@ contains ballsize = sr%ballsize rearrange = sr%rearrange ind => sr%ind(1:) - data => sr%Data(1:,1:) + myData => sr%myData(1:,1:) centeridx = sr%centeridx correltime = sr%correltime nn = sr%nn ! number to search for @@ -1640,7 +1640,7 @@ contains if (rearrange) then sd = 0.0_pReal do k = 1_pInt,dimen - sd = sd + (data(k,i) - qv(k))**2.0_pReal + sd = sd + (myData(k,i) - qv(k))**2.0_pReal if (sd>ballsize) cycle mainloop end do indexofi = ind(i) ! only read it if we have not broken out @@ -1648,7 +1648,7 @@ contains indexofi = ind(i) sd = 0.0_pReal do k = 1_pInt,dimen - sd = sd + (data(k,indexofi) - qv(k))**2.0_pReal + sd = sd + (myData(k,indexofi) - qv(k))**2.0_pReal if (sd>ballsize) cycle mainloop end do endif @@ -1698,12 +1698,12 @@ contains do i = 1_pInt, tp%n if (all_distances(i) 1)then - ileft=ileft-1 + if(ileft > 1_pInt)then + ileft=ileft-1_pInt value=a(ileft); ivalue=ind(ileft) else value=a(iright); ivalue=ind(iright) a(iright)=a(1); ind(iright)=ind(1) - iright=iright-1 - if (iright == 1) then + iright=iright-1_pInt + if (iright == 1_pInt) then a(1)=value;ind(1)=ivalue return endif endif i=ileft - j=2*ileft + j=2_pInt*ileft do while (j <= iright) if(j < iright) then - if(a(j) < a(j+1)) j=j+1 + if(a(j) < a(j+1_pInt)) j=j+1_pInt endif if(value < a(j)) then a(i)=a(j); ind(i)=ind(j) i=j j=j+j else - j=iright+1 + j=iright+1_pInt endif end do a(i)=value; ind(i)=ivalue @@ -1842,7 +1843,7 @@ contains integer(pInt) :: i,j integer(pInt) :: ileft,iright - ileft=n/2+1 + ileft=n/2_pInt+1_pInt iright=n ! do i=1,n @@ -1850,33 +1851,33 @@ contains ! Generate initial idum array ! end do - if(n.eq.1) return + if(n.eq.1_pInt) return do - if(ileft > 1)then - ileft=ileft-1 + if(ileft > 1_pInt)then + ileft=ileft-1_pInt value=a(ileft) else value=a(iright) a(iright)=a(1) - iright=iright-1 - if (iright == 1) then + iright=iright-1_pInt + if (iright == 1_pInt) then a(1) = value return endif endif i=ileft - j=2*ileft + j=2_pInt*ileft do while (j <= iright) if(j < iright) then - if(a(j)%dis < a(j+1)%dis) j=j+1 + if(a(j)%dis < a(j+1_pInt)%dis) j=j+1_pInt endif if(value%dis < a(j)%dis) then a(i)=a(j); i=j j=j+j else - j=iright+1 + j=iright+1_pInt endif end do a(i)=value diff --git a/code/lattice.f90 b/code/lattice.f90 index d4b1ecb32..8c42c5949 100644 --- a/code/lattice.f90 +++ b/code/lattice.f90 @@ -39,11 +39,11 @@ implicit none integer(pInt) lattice_Nhexagonal, & ! # of hexagonal lattice structure (from tag CoverA_ratio) lattice_Nstructure ! # of lattice structures (1: fcc,2: bcc,3+: hexagonal) -integer(pInt), parameter :: lattice_maxNslipFamily = 5 ! max # of slip system families over lattice structures -integer(pInt), parameter :: lattice_maxNtwinFamily = 4 ! max # of twin system families over lattice structures -integer(pInt), parameter :: lattice_maxNslip = 54 ! max # of slip systems over lattice structures -integer(pInt), parameter :: lattice_maxNtwin = 24 ! max # of twin systems over lattice structures -integer(pInt), parameter :: lattice_maxNinteraction = 30 ! max # of interaction types (in hardening matrix part) +integer(pInt), parameter :: lattice_maxNslipFamily = 5_pInt ! max # of slip system families over lattice structures +integer(pInt), parameter :: lattice_maxNtwinFamily = 4_pInt ! max # of twin system families over lattice structures +integer(pInt), parameter :: lattice_maxNslip = 54_pInt ! max # of slip systems over lattice structures +integer(pInt), parameter :: lattice_maxNtwin = 24_pInt ! max # of twin systems over lattice structures +integer(pInt), parameter :: lattice_maxNinteraction = 30_pInt ! max # of interaction types (in hardening matrix part) integer(pInt), pointer, dimension(:,:) :: interactionSlipSlip, & interactionSlipTwin, & @@ -81,10 +81,10 @@ integer(pInt), allocatable, dimension(:,:,:) :: lattice_interactionSlipSlip, & !============================== fcc (1) ================================= - integer(pInt), parameter, dimension(lattice_maxNslipFamily) :: lattice_fcc_NslipSystem = (/12, 0, 0, 0, 0/) - integer(pInt), parameter, dimension(lattice_maxNtwinFamily) :: lattice_fcc_NtwinSystem = (/12, 0, 0, 0/) - integer(pInt), parameter :: lattice_fcc_Nslip = 12 ! sum(lattice_fcc_NslipSystem) - integer(pInt), parameter :: lattice_fcc_Ntwin = 12 ! sum(lattice_fcc_NtwinSystem) + integer(pInt), parameter, dimension(lattice_maxNslipFamily) :: lattice_fcc_NslipSystem = int([12, 0, 0, 0, 0],pInt) + integer(pInt), parameter, dimension(lattice_maxNtwinFamily) :: lattice_fcc_NtwinSystem = int([12, 0, 0, 0],pInt) + integer(pInt), parameter :: lattice_fcc_Nslip = 12_pInt ! sum(lattice_fcc_NslipSystem) + integer(pInt), parameter :: lattice_fcc_Ntwin = 12_pInt ! sum(lattice_fcc_NtwinSystem) integer(pInt) :: lattice_fcc_Nstructure = 0_pInt real(pReal), dimension(3+3,lattice_fcc_Nslip), parameter :: lattice_fcc_systemSlip = & @@ -442,10 +442,10 @@ integer(pInt), allocatable, dimension(:,:,:) :: lattice_interactionSlipSlip, & !============================== hex (3+) ================================= - integer(pInt), parameter, dimension(lattice_maxNslipFamily) :: lattice_hex_NslipSystem = (/ 3, 3, 6,12, 6/) - integer(pInt), parameter, dimension(lattice_maxNtwinFamily) :: lattice_hex_NtwinSystem = (/ 6, 6, 6, 6/) - integer(pInt), parameter :: lattice_hex_Nslip = 30 ! sum(lattice_hex_NslipSystem) - integer(pInt), parameter :: lattice_hex_Ntwin = 24 ! sum(lattice_hex_NtwinSystem) + integer(pInt), parameter, dimension(lattice_maxNslipFamily) :: lattice_hex_NslipSystem = int([ 3, 3, 6,12, 6],pInt) + integer(pInt), parameter, dimension(lattice_maxNtwinFamily) :: lattice_hex_NtwinSystem = int([ 6, 6, 6, 6],pInt) + integer(pInt), parameter :: lattice_hex_Nslip = 30_pInt ! sum(lattice_hex_NslipSystem) + integer(pInt), parameter :: lattice_hex_Ntwin = 24_pInt ! sum(lattice_hex_NtwinSystem) integer(pInt) :: lattice_hex_Nstructure = 0_pInt !* sorted by A. Alankar & P. Eisenlohr @@ -824,12 +824,12 @@ function lattice_initializeStructure(struct,CoverA) lattice_fcc_Nstructure = lattice_fcc_Nstructure + 1_pInt ! count fcc instances if (lattice_fcc_Nstructure == 1_pInt) then ! me is first fcc structure processMe = .true. - do i = 1,myNslip ! calculate slip system vectors + do i = 1_pInt,myNslip ! calculate slip system vectors sd(1:3,i) = lattice_fcc_systemSlip(1:3,i)/sqrt(math_mul3x3(lattice_fcc_systemSlip(1:3,i),lattice_fcc_systemSlip(1:3,i))) sn(1:3,i) = lattice_fcc_systemSlip(4:6,i)/sqrt(math_mul3x3(lattice_fcc_systemSlip(4:6,i),lattice_fcc_systemSlip(4:6,i))) st(1:3,i) = math_vectorproduct(sd(1:3,i),sn(1:3,i)) enddo - do i = 1,myNtwin ! calculate twin system vectors and (assign) shears + do i = 1_pInt,myNtwin ! calculate twin system vectors and (assign) shears td(1:3,i) = lattice_fcc_systemTwin(1:3,i)/sqrt(math_mul3x3(lattice_fcc_systemTwin(1:3,i),lattice_fcc_systemTwin(1:3,i))) tn(1:3,i) = lattice_fcc_systemTwin(4:6,i)/sqrt(math_mul3x3(lattice_fcc_systemTwin(4:6,i),lattice_fcc_systemTwin(4:6,i))) tt(1:3,i) = math_vectorproduct(td(1:3,i),tn(1:3,i)) @@ -850,12 +850,12 @@ function lattice_initializeStructure(struct,CoverA) lattice_bcc_Nstructure = lattice_bcc_Nstructure + 1_pInt ! count bcc instances if (lattice_bcc_Nstructure == 1_pInt) then ! me is first bcc structure processMe = .true. - do i = 1,myNslip ! calculate slip system vectors + do i = 1_pInt,myNslip ! calculate slip system vectors sd(1:3,i) = lattice_bcc_systemSlip(1:3,i)/sqrt(math_mul3x3(lattice_bcc_systemSlip(1:3,i),lattice_bcc_systemSlip(1:3,i))) sn(1:3,i) = lattice_bcc_systemSlip(4:6,i)/sqrt(math_mul3x3(lattice_bcc_systemSlip(4:6,i),lattice_bcc_systemSlip(4:6,i))) st(1:3,i) = math_vectorproduct(sd(1:3,i),sn(1:3,i)) enddo - do i = 1,myNtwin ! calculate twin system vectors and (assign) shears + do i = 1_pInt,myNtwin ! calculate twin system vectors and (assign) shears td(1:3,i) = lattice_bcc_systemTwin(1:3,i)/sqrt(math_mul3x3(lattice_bcc_systemTwin(1:3,i),lattice_bcc_systemTwin(1:3,i))) tn(1:3,i) = lattice_bcc_systemTwin(4:6,i)/sqrt(math_mul3x3(lattice_bcc_systemTwin(4:6,i),lattice_bcc_systemTwin(4:6,i))) tt(1:3,i) = math_vectorproduct(td(1:3,i),tn(1:3,i)) @@ -877,7 +877,7 @@ function lattice_initializeStructure(struct,CoverA) myNtwin = lattice_hex_Ntwin ! overall number of twin systems processMe = .true. ! converting from 4 axes coordinate system (a1=a2=a3=c) to ortho-hexgonal system (a, b, c) - do i = 1,myNslip + do i = 1_pInt,myNslip hex_d(1) = lattice_hex_systemSlip(1,i)*1.5_pReal ! direction [uvtw]->[3u/2 (u+2v)*sqrt(3)/2 w*(c/a)] hex_d(2) = (lattice_hex_systemSlip(1,i)+2.0_pReal*lattice_hex_systemSlip(2,i))*(0.5_pReal*sqrt(3.0_pReal)) hex_d(3) = lattice_hex_systemSlip(4,i)*CoverA @@ -889,7 +889,7 @@ function lattice_initializeStructure(struct,CoverA) sn(1:3,i) = hex_n/sqrt(math_mul3x3(hex_n,hex_n)) st(1:3,i) = math_vectorproduct(sd(1:3,i),sn(1:3,i)) enddo - do i = 1,myNtwin + do i = 1_pInt,myNtwin hex_d(1) = lattice_hex_systemTwin(1,i)*1.5_pReal hex_d(2) = (lattice_hex_systemTwin(1,i)+2.0_pReal*lattice_hex_systemTwin(2,i))*(0.5_pReal*sqrt(3.0_pReal)) hex_d(3) = lattice_hex_systemTwin(4,i)*CoverA @@ -923,14 +923,14 @@ function lattice_initializeStructure(struct,CoverA) if (processMe) then if (myStructure > lattice_Nstructure) & call IO_error(666_pInt,0_pInt,0_pInt,0_pInt,'structure index too large') ! check for memory leakage - do i = 1,myNslip ! store slip system vectors and Schmid matrix for my structure + do i = 1_pInt,myNslip ! store slip system vectors and Schmid matrix for my structure lattice_sd(1:3,i,myStructure) = sd(1:3,i) lattice_st(1:3,i,myStructure) = st(1:3,i) lattice_sn(1:3,i,myStructure) = sn(1:3,i) lattice_Sslip(1:3,1:3,i,myStructure) = math_tensorproduct(sd(1:3,i),sn(1:3,i)) lattice_Sslip_v(1:6,i,myStructure) = math_Mandel33to6(math_symmetric33(lattice_Sslip(1:3,1:3,i,myStructure))) enddo - do i = 1,myNtwin ! store twin system vectors and Schmid plus rotation matrix for my structure + do i = 1_pInt,myNtwin ! store twin system vectors and Schmid plus rotation matrix for my structure lattice_td(1:3,i,myStructure) = td(1:3,i) lattice_tt(1:3,i,myStructure) = tt(1:3,i) lattice_tn(1:3,i,myStructure) = tn(1:3,i) diff --git a/code/makefile b/code/makefile index 68d22d4bb..1e848fd46 100644 --- a/code/makefile +++ b/code/makefile @@ -32,42 +32,13 @@ # PREFIX = arbitrary prefix # SUFFIX = arbitrary suffix # STANDARD_CHECK = checking for Fortran 2008, compiler dependend -######################################################################################## -# Here are some useful debugging switches for ifort. Switch on by uncommenting the #SUFFIX line at the end of this section: -# information on http://software.intel.com/en-us/articles/determining-root-cause-of-sigsegv-or-sigbus-errors/ - -# check if an array index is too small (<1) or too large! -DEBUG1 =-check bounds -g - -#will cause a lot of warnings because we create a bunch of temporary arrays -DEBUG2 =-check arg_temp_created - -#check from time to time -DEBUG3 =-fp-stack-check -g -traceback -gen-interfaces -warn interfaces - -#should not be done for OpenMP, but set "ulimit -s unlimited" on shell. Probably it helps also to unlimit other limits -DEBUG4 =-heap-arrays - -#additional warnings -DEBUG5 =-warn all -# or one or more of those: alignments, declarations,general, ignore_loc, uncalled, unuses, usage - -#set precision (check for missing _pInt and _pReal) -DEBUG6= -real-size 32 -integer-size 16 -#or one of those 16/32/64/128 (= 2,4,8,16 bytes) -#SUFFIX =$(DEBUG1) $(DEBUG2) $(DEBUG3) $(DEBUG4) $(DEBUG5) $(DEBUG6) - -# Here are some useful debugging switches for gfortran -# fcheck-bounds: eqv to DEBUG1 of ifort - - ######################################################################################## #auto values will be set by setup_code.py -FFTWROOT := $(DAMASK_ROOT)/lib/fftw -IKMLROOT := -ACMLROOT := /opt/acml4.4.0 -LAPACKROOT := +FFTWROOT :=/$(DAMASK_ROOT)/lib/fftw +IKMLROOT := +ACMLROOT :=/opt/acml4.4.0 +#LAPACKROOT := /usr F90 ?= ifort COMPILERNAME ?= $(F90) @@ -144,38 +115,71 @@ OPTIMIZATION_OFF_ifort :=-O0 OPTIMIZATION_OFF_gfortran :=-O0 OPTIMIZATION_DEFENSIVE_ifort :=-O2 OPTIMIZATION_DEFENSIVE_gfortran :=-O2 -OPTIMIZATION_AGGRESSIVE_ifort :=-O3 $(PORTABLE_SWITCH) -ip -static -fp-model fast=2 -no-prec-div +OPTIMIZATION_AGGRESSIVE_ifort :=-O3 $(PORTABLE_SWITCH) -ipo -static -no-prec-div -fp-model fast=2 OPTIMIZATION_AGGRESSIVE_gfortran :=-O3 $(PORTABLE_SWITCH) -ffast-math -funroll-loops -ftree-vectorize COMPILE_OPTIONS_ifort :=-fpp\ + -implicitnone\ -diag-enable sc3\ - -diag-disable 8291,8290,5268\ + -diag-disable 5268\ -warn declarations\ -warn general\ - -warn usage + -warn usage\ + -warn interfaces\ + -warn ignore_loc\ + -warn alignments\ + -warn unused\ + -warn errors\ + -warn stderrors + +#-fpp: preprocessor +#-fimplicit-none: assume "implicit-none" even if not present in source +#-diag-disable: disables warnings, where +# warning ID 5268: the text exceeds right hand column allowed on the line (we have only comments there) +#-warn: enables warnings, where +# declarations: any undeclared names +# general: warning messages and informational messages are issued by the compiler +# usage: questionable programming practices +# interfaces: checks the interfaces of all SUBROUTINEs called and FUNCTIONs invoked in your compilation against an external set of interface blocks +# ignore_loc: %LOC is stripped from an actual argument +# alignments: data that is not naturally aligned +# unused: declared variables that are never used +# errors: warnings are changed to errors +# stderrors: warnings about Fortran standard violations are changed to errors +# +################################################################################################### +#MORE OPTIONS FOR DEBUGGING DURING COMPILING +#-warn: enables warnings, where +# truncated_source: Determines whether warnings occur when source exceeds the maximum column width in fixed-format files. (too many warnings because we have comments beyond character 132) +# uncalled: Determines whether warnings occur when a statement function is never called +# all: +# +#OPTIONS FOR DEGUBBING DURING RUNTIME +# information on http://software.intel.com/en-us/articles/determining-root-cause-of-sigsegv-or-sigbus-errors/ +#-g: Generate symbolic debugging information in the object file +#-traceback: Generate extra information in the object file to provide source file traceback information when a severe error occurs at run time. +#-gen-interfaces: Generate an interface block for each routine. http://software.intel.com/en-us/blogs/2012/01/05/doctor-fortran-gets-explicit-again/ +#-fp-stack-check: Generate extra code after every function call to ensure that the floating-point (FP) stack is in the expected state. +#-check: checks at runtime, where +# bounds: check if an array index is too small (<1) or too large! +# arg_temp_created: will cause a lot of warnings because we create a bunch of temporary arrays +# format: Checking for the data type of an item being formatted for output. +# output_conversion: Checking for the fit of data items within a designated format descriptor field. +# pointers: Checking for certain disassociated or uninitialized pointers or unallocated allocatable objects. +# uninit: Checking for uninitialized variables. +#-heap-arrays: should not be done for OpenMP, but set "ulimit -s unlimited" on shell. Probably it helps also to unlimit other limits +# +#OPTIONS FOR TYPE DEBUGGING +#-real-size 32: set precision to one of those 32/64/128 (= 4/8/16 bytes) for standard real (=8 for pReal) +#-integer-size 16: set precision to one of those 16/32/64 (= 2/4/8 bytes) for standard integer (=4 for pInt) +################################################################################################### -#alignments: Determines whether warnings occur for data that is not naturally aligned. -#declarations: Determines whether warnings occur for any undeclared names. -#errors: Determines whether warnings are changed to errors. -#general: Determines whether warning messages and informational messages are issued by the compiler. -#ignore_loc: Determines whether warnings occur when %LOC is stripped from an actual argument. -#interfaces: Determines whether the compiler checks the interfaces of all SUBROUTINEs called and FUNCTIONs invoked in your compilation against an external set of interface blocks. -#stderrors: Determines whether warnings about Fortran standard violations are changed to errors. -#truncated_source: Determines whether warnings occur when source exceeds the maximum column width in fixed-format files. -#uncalled: Determines whether warnings occur when a statement function is never called -#unused: Determines whether warnings occur for declared variables that are never used. -#usage: Determines whether warnings occur for questionable programming practices. - -#-fpp: preprocessor -#-diag-disable: disables warnings, where -# warning ID 9291: -# warning ID 8290: -# warning ID 5268: The text exceeds right hand column allowed on the line (we have only comments there) COMPILE_OPTIONS_gfortran :=-xf95-cpp-input\ -ffree-line-length-132\ -fno-range-check\ -fimplicit-none\ + -fall-intrinsics\ -pedantic\ -Warray-bounds\ -Wunused-parameter\ @@ -189,47 +193,50 @@ COMPILE_OPTIONS_gfortran :=-xf95-cpp-input\ -Wunused-value\ -Wunderflow -#-xf95-cpp-input: preprocessor -#-ffree-line-length-132: restrict line length to the standard 132 characters -#-fno-range-check: disables checking if result can be represented by variable. Needs to be set to enable DAMASK_NaN -#-fimplicit-none: assume "implicit-none" even if not present in source -#-pedantic: more strict on standard, enables some of the warnings below -#-Warray-bounds: checks if array reference is out of bounds at compile time. use -fcheck-bounds to also check during runtime -#-Wunused-parameter: find usused variables with "parameter" attribute -#-Wampersand: checks if a character expression is continued proberly by an ampersand at the end of the line and at the beginning of the new line -#-Wno-tabs: do not allow tabs in source -#-Wcharacter-truncation: warn if character expressions (strings) are truncated -#-Wintrinsic-shadow: warn if a user-defined procedure or module procedure has the same name as an intrinsic -#-Waliasing: warn about possible aliasing of dummy arguments. Specifically, it warns if the same actual argument is associated with a dummy argument with "INTENT(IN)" and a dummy argument with "INTENT(OUT)" in a call with an explicit interface. -#-Wconversion: warn about implicit conversions between different type -#-Wsurprising: warn when "suspicious" code constructs are encountered. While technically legal these usually indicate that an error has been made. +#-xf95-cpp-input: preprocessor +#-ffree-line-length-132: restrict line length to the standard 132 characters +#-fno-range-check: disables checking if result can be represented by variable. Needs to be set to enable DAMASK_NaN +#-fimplicit-none: assume "implicit-none" even if not present in source +#-fall-intrinsics: +#-pedantic: more strict on standard, enables some of the warnings below +#-Warray-bounds: checks if array reference is out of bounds at compile time. use -fcheck-bounds to also check during runtime +#-Wunused-parameter: find usused variables with "parameter" attribute +#-Wampersand: checks if a character expression is continued proberly by an ampersand at the end of the line and at the beginning of the new line +#-Wno-tabs: do not allow tabs in source +#-Wcharacter-truncation: warn if character expressions (strings) are truncated +#-Wintrinsic-shadow: warn if a user-defined procedure or module procedure has the same name as an intrinsic +#-Waliasing: warn about possible aliasing of dummy arguments. Specifically, it warns if the same actual argument is associated with a dummy argument with "INTENT(IN)" and a dummy argument with "INTENT(OUT)" in a call with an explicit interface. +#-Wconversion: warn about implicit conversions between different type +#-Wsurprising: warn when "suspicious" code constructs are encountered. While technically legal these usually indicate that an error has been made. #-Wunused-value: -#-Wunderflow: produce a warning when numerical constant expressions are encountered, which yield an UNDERFLOW during compilation - - - -#MORE OPTIONS - -# only for gfortran 4.6: - #-Wsuggest-attribute=const - #-Wsuggest-attribute=noreturn - #-Wsuggest-attribute=pure - -# too many warnings because we have comments beyond character 132: - #-Wline-truncation -# warnings because of "flush" is not longer in the standard, but still an intrinsic fuction of the compilers: - #-Wintrinsic-std -# warnings because we have many temporary arrays (performance issue?): - #-Warray-temporaries - -# -Wimplicit-interface -# -pedantic-errors -# -fmodule-private - +#-Wunderflow: produce a warning when numerical constant expressions are encountered, which yield an UNDERFLOW during compilation +# +################################################################################################### +#OPTIONS FOR GFORTRAN 4.6 +#-Wsuggest-attribute=const: +#-Wsuggest-attribute=noreturn: +#-Wsuggest-attribute=pure: +# +#MORE OPTIONS FOR DEBUGGING DURING COMPILING +#-Wline-truncation: too many warnings because we have comments beyond character 132 +#-Wintrinsic-std: warnings because of "flush" is not longer in the standard, but still an intrinsic fuction of the compilers: +#-Warray-temporarieswarnings: +# because we have many temporary arrays (performance issue?): +#-Wimplicit-interface +#-pedantic-errors +#-fmodule-private +# +#OPTIONS FOR DEGUBBING DURING RUNTIME +#-fcheck-bounds: check if an array index is too small (<1) or too large! +# +#OPTIONS FOR TYPE DEBUGGING +#-fdefault-real-8: set precision to 8 bytes for standard real (=8 for pReal). Will set size of double to 16 bytes as long as -fdefault-double-8 is not set +#-fdefault-integer-8: set precision to 8 bytes for standard integer (=4 for pInt) +################################################################################################## COMPILE =$(OPENMP_FLAG_$(F90)) $(COMPILE_OPTIONS_$(F90)) $(STANDARD_CHECK_$(F90)) $(OPTIMIZATION_$(OPTI)_$(F90)) -c COMPILE_MAXOPTI =$(OPENMP_FLAG_$(F90)) $(COMPILE_OPTIONS_$(F90)) $(STANDARD_CHECK_$(F90)) $(OPTIMIZATION_$(MAXOPTI)_$(F90)) -c - +################################################################################################### DAMASK_spectral.exe: DAMASK_spectral.o CPFEM.a $(PREFIX) $(COMPILERNAME) ${OPENMP_FLAG_${F90}} -o DAMASK_spectral.exe DAMASK_spectral.o CPFEM.a \ diff --git a/code/material.f90 b/code/material.f90 index 30d5a874b..15404db82 100644 --- a/code/material.f90 +++ b/code/material.f90 @@ -114,7 +114,7 @@ subroutine material_init() implicit none !* Definition of variables - integer(pInt), parameter :: fileunit = 200 + integer(pInt), parameter :: fileunit = 200_pInt integer(pInt) i,j !$OMP CRITICAL (write2out) @@ -128,46 +128,46 @@ subroutine material_init() call IO_open_file(fileunit,material_configFile) ! ...open material.config file endif call material_parseHomogenization(fileunit,material_partHomogenization) - if (debug_verbosity > 0) then + if (debug_verbosity > 0_pInt) then !$OMP CRITICAL (write2out) write (6,*) 'Homogenization parsed' !$OMP END CRITICAL (write2out) endif call material_parseMicrostructure(fileunit,material_partMicrostructure) - if (debug_verbosity > 0) then + if (debug_verbosity > 0_pInt) then !$OMP CRITICAL (write2out) write (6,*) 'Microstructure parsed' !$OMP END CRITICAL (write2out) endif call material_parseCrystallite(fileunit,material_partCrystallite) - if (debug_verbosity > 0) then + if (debug_verbosity > 0_pInt) then !$OMP CRITICAL (write2out) write (6,*) 'Crystallite parsed' !$OMP END CRITICAL (write2out) endif call material_parseTexture(fileunit,material_partTexture) - if (debug_verbosity > 0) then + if (debug_verbosity > 0_pInt) then !$OMP CRITICAL (write2out) write (6,*) 'Texture parsed' !$OMP END CRITICAL (write2out) endif call material_parsePhase(fileunit,material_partPhase) - if (debug_verbosity > 0) then + if (debug_verbosity > 0_pInt) then !$OMP CRITICAL (write2out) write (6,*) 'Phase parsed' !$OMP END CRITICAL (write2out) endif close(fileunit) - do i = 1,material_Nmicrostructure - if (microstructure_crystallite(i) < 1 .or. & + do i = 1_pInt,material_Nmicrostructure + if (microstructure_crystallite(i) < 1_pInt .or. & microstructure_crystallite(i) > material_Ncrystallite) call IO_error(150_pInt,i) - if (minval(microstructure_phase(1:microstructure_Nconstituents(i),i)) < 1 .or. & + if (minval(microstructure_phase(1:microstructure_Nconstituents(i),i)) < 1_pInt .or. & maxval(microstructure_phase(1:microstructure_Nconstituents(i),i)) > material_Nphase) call IO_error(151_pInt,i) - if (minval(microstructure_texture(1:microstructure_Nconstituents(i),i)) < 1 .or. & + if (minval(microstructure_texture(1:microstructure_Nconstituents(i),i)) < 1_pInt .or. & maxval(microstructure_texture(1:microstructure_Nconstituents(i),i)) > material_Ntexture) call IO_error(152_pInt,i) if (abs(sum(microstructure_fraction(:,i)) - 1.0_pReal) >= 1.0e-10_pReal) then - if (debug_verbosity > 0) then + if (debug_verbosity > 0_pInt) then !$OMP CRITICAL (write2out) write(6,*)'sum of microstructure fraction = ',sum(microstructure_fraction(:,i)) !$OMP END CRITICAL (write2out) @@ -175,25 +175,25 @@ subroutine material_init() call IO_error(153_pInt,i) endif enddo - if (debug_verbosity > 0) then + if (debug_verbosity > 0_pInt) then !$OMP CRITICAL (write2out) write (6,*) write (6,*) 'MATERIAL configuration' write (6,*) write (6,'(a32,1x,a16,1x,a6)') 'homogenization ','type ','grains' - do i = 1,material_Nhomogenization + do i = 1_pInt,material_Nhomogenization write (6,'(1x,a32,1x,a16,1x,i4)') homogenization_name(i),homogenization_type(i),homogenization_Ngrains(i) enddo write (6,*) write (6,'(a32,1x,a11,1x,a12,1x,a13)') 'microstructure ','crystallite','constituents','homogeneous' - do i = 1,material_Nmicrostructure + do i = 1_pInt,material_Nmicrostructure write (6,'(a32,4x,i4,8x,i4,8x,l1)') microstructure_name(i), & microstructure_crystallite(i), & microstructure_Nconstituents(i), & microstructure_elemhomo(i) if (microstructure_Nconstituents(i) > 0_pInt) then - do j = 1,microstructure_Nconstituents(i) - write (6,'(a1,1x,a32,1x,a32,1x,f6.4)') '>',phase_name(microstructure_phase(j,i)),& + do j = 1_pInt,microstructure_Nconstituents(i) + write (6,'(a1,1x,a32,1x,a32,1x,f7.4)') '>',phase_name(microstructure_phase(j,i)),& texture_name(microstructure_texture(j,i)),& microstructure_fraction(j,i) enddo @@ -209,7 +209,7 @@ endsubroutine !********************************************************************* -subroutine material_parseHomogenization(file,myPart) +subroutine material_parseHomogenization(myFile,myPart) !********************************************************************* use prec, only: pInt @@ -218,14 +218,14 @@ subroutine material_parseHomogenization(file,myPart) implicit none character(len=*), intent(in) :: myPart - integer(pInt), intent(in) :: file + integer(pInt), intent(in) :: myFile integer(pInt), parameter :: maxNchunks = 2_pInt integer(pInt), dimension(1+2*maxNchunks) :: positions integer(pInt) Nsections, section, s character(len=64) tag character(len=1024) line - Nsections = IO_countSections(file,myPart) + Nsections = IO_countSections(myFile,myPart) material_Nhomogenization = Nsections if (Nsections < 1_pInt) call IO_error(160_pInt,ext_msg=myPart) @@ -236,23 +236,23 @@ subroutine material_parseHomogenization(file,myPart) allocate(homogenization_Noutput(Nsections)); homogenization_Noutput = 0_pInt allocate(homogenization_active(Nsections)); homogenization_active = .false. - forall (s = 1:Nsections) homogenization_active(s) = any(mesh_element(3,:) == s) ! current homogenization used in model? Homogenization view, maximum operations depend on maximum number of homog schemes - homogenization_Noutput = IO_countTagInPart(file,myPart,'(output)',Nsections) + forall (s = 1_pInt:Nsections) homogenization_active(s) = any(mesh_element(3,:) == s) ! current homogenization used in model? Homogenization view, maximum operations depend on maximum number of homog schemes + homogenization_Noutput = IO_countTagInPart(myFile,myPart,'(output)',Nsections) - rewind(file) + rewind(myFile) line = '' - section = 0 + section = 0_pInt do while (IO_lc(IO_getTag(line,'<','>')) /= myPart) ! wind forward to myPart - read(file,'(a1024)',END=100) line + read(myFile,'(a1024)',END=100) line enddo do - read(file,'(a1024)',END=100) line + read(myFile,'(a1024)',END=100) line if (IO_isBlank(line)) cycle ! skip empty lines if (IO_getTag(line,'<','>') /= '') exit ! stop at next part if (IO_getTag(line,'[',']') /= '') then ! next section - section = section + 1 + section = section + 1_pInt homogenization_name(section) = IO_getTag(line,'[',']') endif if (section > 0_pInt) then @@ -277,7 +277,7 @@ subroutine material_parseHomogenization(file,myPart) !********************************************************************* -subroutine material_parseMicrostructure(file,myPart) +subroutine material_parseMicrostructure(myFile,myPart) !********************************************************************* use prec, only: pInt @@ -286,14 +286,14 @@ subroutine material_parseMicrostructure(file,myPart) implicit none character(len=*), intent(in) :: myPart - integer(pInt), intent(in) :: file + integer(pInt), intent(in) :: myFile integer(pInt), parameter :: maxNchunks = 7_pInt integer(pInt), dimension(1_pInt+2_pInt*maxNchunks) :: positions integer(pInt) Nsections, section, constituent, e, i character(len=64) tag character(len=1024) line - Nsections = IO_countSections(file,myPart) + Nsections = IO_countSections(myFile,myPart) material_Nmicrostructure = Nsections if (Nsections < 1_pInt) call IO_error(160_pInt,ext_msg=myPart) @@ -303,42 +303,42 @@ subroutine material_parseMicrostructure(file,myPart) allocate(microstructure_active(Nsections)) allocate(microstructure_elemhomo(Nsections)) - forall (e = 1:mesh_NcpElems) microstructure_active(mesh_element(4,e)) = .true. ! current microstructure used in model? Elementwise view, maximum N operations for N elements + forall (e = 1_pInt:mesh_NcpElems) microstructure_active(mesh_element(4,e)) = .true. ! current microstructure used in model? Elementwise view, maximum N operations for N elements - microstructure_Nconstituents = IO_countTagInPart(file,myPart,'(constituent)',Nsections) + microstructure_Nconstituents = IO_countTagInPart(myFile,myPart,'(constituent)',Nsections) microstructure_maxNconstituents = maxval(microstructure_Nconstituents) - microstructure_elemhomo = IO_spotTagInPart(file,myPart,'/elementhomogeneous/',Nsections) + microstructure_elemhomo = IO_spotTagInPart(myFile,myPart,'/elementhomogeneous/',Nsections) allocate(microstructure_phase (microstructure_maxNconstituents,Nsections)); microstructure_phase = 0_pInt allocate(microstructure_texture (microstructure_maxNconstituents,Nsections)); microstructure_texture = 0_pInt allocate(microstructure_fraction(microstructure_maxNconstituents,Nsections)); microstructure_fraction = 0.0_pReal - rewind(file) + rewind(myFile) line = '' - section = 0 + section = 0_pInt do while (IO_lc(IO_getTag(line,'<','>')) /= myPart) ! wind forward to myPart - read(file,'(a1024)',END=100) line + read(myFile,'(a1024)',END=100) line enddo do - read(file,'(a1024)',END=100) line + read(myFile,'(a1024)',END=100) line if (IO_isBlank(line)) cycle ! skip empty lines if (IO_getTag(line,'<','>') /= '') exit ! stop at next part if (IO_getTag(line,'[',']') /= '') then ! next section - section = section + 1 - constituent = 0 + section = section + 1_pInt + constituent = 0_pInt microstructure_name(section) = IO_getTag(line,'[',']') endif - if (section > 0) then + if (section > 0_pInt) then positions = IO_stringPos(line,maxNchunks) tag = IO_lc(IO_stringValue(line,positions,1_pInt)) ! extract key select case(tag) case ('crystallite') microstructure_crystallite(section) = IO_intValue(line,positions,2_pInt) case ('(constituent)') - constituent = constituent + 1 - do i=2,6,2 + constituent = constituent + 1_pInt + do i=2_pInt,6_pInt,2_pInt tag = IO_lc(IO_stringValue(line,positions,i)) select case (tag) case('phase') @@ -357,7 +357,7 @@ subroutine material_parseMicrostructure(file,myPart) !********************************************************************* -subroutine material_parseCrystallite(file,myPart) +subroutine material_parseCrystallite(myFile,myPart) !********************************************************************* use prec, only: pInt @@ -366,33 +366,33 @@ subroutine material_parseCrystallite(file,myPart) implicit none character(len=*), intent(in) :: myPart - integer(pInt), intent(in) :: file + integer(pInt), intent(in) :: myFile integer(pInt) Nsections, section character(len=1024) line - Nsections = IO_countSections(file,myPart) + Nsections = IO_countSections(myFile,myPart) material_Ncrystallite = Nsections if (Nsections < 1_pInt) call IO_error(160_pInt,ext_msg=myPart) allocate(crystallite_name(Nsections)); crystallite_name = '' allocate(crystallite_Noutput(Nsections)); crystallite_Noutput = 0_pInt - crystallite_Noutput = IO_countTagInPart(file,myPart,'(output)',Nsections) + crystallite_Noutput = IO_countTagInPart(myFile,myPart,'(output)',Nsections) - rewind(file) + rewind(myFile) line = '' - section = 0 + section = 0_pInt do while (IO_lc(IO_getTag(line,'<','>')) /= myPart) ! wind forward to myPart - read(file,'(a1024)',END=100) line + read(myFile,'(a1024)',END=100) line enddo do - read(file,'(a1024)',END=100) line + read(myFile,'(a1024)',END=100) line if (IO_isBlank(line)) cycle ! skip empty lines if (IO_getTag(line,'<','>') /= '') exit ! stop at next part if (IO_getTag(line,'[',']') /= '') then ! next section - section = section + 1 + section = section + 1_pInt crystallite_name(section) = IO_getTag(line,'[',']') endif enddo @@ -401,7 +401,7 @@ subroutine material_parseCrystallite(file,myPart) !********************************************************************* -subroutine material_parsePhase(file,myPart) +subroutine material_parsePhase(myFile,myPart) !********************************************************************* use prec, only: pInt @@ -409,14 +409,14 @@ subroutine material_parsePhase(file,myPart) implicit none character(len=*), intent(in) :: myPart - integer(pInt), intent(in) :: file - integer(pInt), parameter :: maxNchunks = 2 + integer(pInt), intent(in) :: myFile + integer(pInt), parameter :: maxNchunks = 2_pInt integer(pInt), dimension(1+2*maxNchunks) :: positions integer(pInt) Nsections, section, s character(len=64) tag character(len=1024) line - Nsections = IO_countSections(file,myPart) + Nsections = IO_countSections(myFile,myPart) material_Nphase = Nsections if (Nsections < 1_pInt) call IO_error(160_pInt,ext_msg=myPart) @@ -426,23 +426,23 @@ subroutine material_parsePhase(file,myPart) allocate(phase_Noutput(Nsections)) allocate(phase_localConstitution(Nsections)) - phase_Noutput = IO_countTagInPart(file,myPart,'(output)',Nsections) - phase_localConstitution = .not. IO_spotTagInPart(file,myPart,'/nonlocal/',Nsections) + phase_Noutput = IO_countTagInPart(myFile,myPart,'(output)',Nsections) + phase_localConstitution = .not. IO_spotTagInPart(myFile,myPart,'/nonlocal/',Nsections) - rewind(file) + rewind(myFile) line = '' - section = 0 + section = 0_pInt do while (IO_lc(IO_getTag(line,'<','>')) /= myPart) ! wind forward to myPart - read(file,'(a1024)',END=100) line + read(myFile,'(a1024)',END=100) line enddo do - read(file,'(a1024)',END=100) line + read(myFile,'(a1024)',END=100) line if (IO_isBlank(line)) cycle ! skip empty lines if (IO_getTag(line,'<','>') /= '') exit ! stop at next part if (IO_getTag(line,'[',']') /= '') then ! next section - section = section + 1 + section = section + 1_pInt phase_name(section) = IO_getTag(line,'[',']') endif if (section > 0_pInt) then @@ -451,9 +451,9 @@ subroutine material_parsePhase(file,myPart) select case(tag) case ('constitution') phase_constitution(section) = IO_lc(IO_stringValue(line,positions,2_pInt)) - do s = 1,section + do s = 1_pInt,section if (phase_constitution(s) == phase_constitution(section)) & - phase_constitutionInstance(section) = phase_constitutionInstance(section) + 1 ! count instances + phase_constitutionInstance(section) = phase_constitutionInstance(section) + 1_pInt ! count instances enddo end select endif @@ -463,7 +463,7 @@ subroutine material_parsePhase(file,myPart) !********************************************************************* -subroutine material_parseTexture(file,myPart) +subroutine material_parseTexture(myFile,myPart) !********************************************************************* use prec, only: pInt, pReal @@ -472,15 +472,15 @@ subroutine material_parseTexture(file,myPart) implicit none character(len=*), intent(in) :: myPart - integer(pInt), intent(in) :: file - integer(pInt), parameter :: maxNchunks = 13 + integer(pInt), intent(in) :: myFile + integer(pInt), parameter :: maxNchunks = 13_pInt integer(pInt), dimension(1+2*maxNchunks) :: positions integer(pInt) Nsections, section, gauss, fiber, i character(len=64) tag character(len=1024) line - Nsections = IO_countSections(file,myPart) + Nsections = IO_countSections(myFile,myPart) material_Ntexture = Nsections if (Nsections < 1_pInt) call IO_error(160_pInt,ext_msg=myPart) @@ -490,33 +490,33 @@ subroutine material_parseTexture(file,myPart) allocate(texture_Ngauss(Nsections)); texture_Ngauss = 0_pInt allocate(texture_Nfiber(Nsections)); texture_Nfiber = 0_pInt - texture_Ngauss = IO_countTagInPart(file,myPart,'(gauss)', Nsections) + & - IO_countTagInPart(file,myPart,'(random)',Nsections) - texture_Nfiber = IO_countTagInPart(file,myPart,'(fiber)', Nsections) + texture_Ngauss = IO_countTagInPart(myFile,myPart,'(gauss)', Nsections) + & + IO_countTagInPart(myFile,myPart,'(random)',Nsections) + texture_Nfiber = IO_countTagInPart(myFile,myPart,'(fiber)', Nsections) texture_maxNgauss = maxval(texture_Ngauss) texture_maxNfiber = maxval(texture_Nfiber) allocate(texture_Gauss (5,texture_maxNgauss,Nsections)); texture_Gauss = 0.0_pReal allocate(texture_Fiber (6,texture_maxNfiber,Nsections)); texture_Fiber = 0.0_pReal - rewind(file) + rewind(myFile) line = '' - section = 0 + section = 0_pInt do while (IO_lc(IO_getTag(line,'<','>')) /= myPart) ! wind forward to myPart - read(file,'(a1024)',END=100) line + read(myFile,'(a1024)',END=100) line enddo do - read(file,'(a1024)',END=100) line + read(myFile,'(a1024)',END=100) line if (IO_isBlank(line)) cycle ! skip empty lines if (IO_getTag(line,'<','>') /= '') exit ! stop at next part if (IO_getTag(line,'[',']') /= '') then ! next section - section = section + 1 - gauss = 0 - fiber = 0 + section = section + 1_pInt + gauss = 0_pInt + fiber = 0_pInt texture_name(section) = IO_getTag(line,'[',']') endif - if (section > 0) then + if (section > 0_pInt) then positions = IO_stringPos(line,maxNchunks) tag = IO_lc(IO_stringValue(line,positions,1_pInt)) ! extract key select case(tag) @@ -528,17 +528,17 @@ subroutine material_parseTexture(file,myPart) tag = IO_lc(IO_stringValue(line,positions,2_pInt)) select case (tag) case('orthotropic') - texture_symmetry(section) = 4 + texture_symmetry(section) = 4_pInt case('monoclinic') - texture_symmetry(section) = 2 + texture_symmetry(section) = 2_pInt case default - texture_symmetry(section) = 1 + texture_symmetry(section) = 1_pInt end select case ('(random)') - gauss = gauss + 1 + gauss = gauss + 1_pInt texture_Gauss(1:3,gauss,section) = math_sampleRandomOri() - do i = 2,4,2 + do i = 2_pInt,4_pInt,2_pInt tag = IO_lc(IO_stringValue(line,positions,i)) select case (tag) case('scatter') @@ -549,8 +549,8 @@ subroutine material_parseTexture(file,myPart) enddo case ('(gauss)') - gauss = gauss + 1 - do i = 2,10,2 + gauss = gauss + 1_pInt + do i = 2_pInt,10_pInt,2_pInt tag = IO_lc(IO_stringValue(line,positions,i)) select case (tag) case('phi1') @@ -567,8 +567,8 @@ subroutine material_parseTexture(file,myPart) enddo case ('(fiber)') - fiber = fiber + 1 - do i = 2,12,2 + fiber = fiber + 1_pInt + do i = 2_pInt,12_pInt,2_pInt tag = IO_lc(IO_stringValue(line,positions,i)) select case (tag) case('alpha1') @@ -629,7 +629,7 @@ subroutine material_populateGrains() allocate(Nelems(material_Nhomogenization,material_Nmicrostructure)); Nelems = 0_pInt ! precounting of elements for each homog/micro pair - do e = 1, mesh_NcpElems + do e = 1_pInt, mesh_NcpElems homog = mesh_element(3,e) micro = mesh_element(4,e) Nelems(homog,micro) = Nelems(homog,micro) + 1_pInt @@ -641,12 +641,12 @@ subroutine material_populateGrains() Nelems = 0_pInt ! reuse as counter ! identify maximum grain count per IP (from element) and find grains per homog/micro pair - do e = 1,mesh_NcpElems + do e = 1_pInt,mesh_NcpElems homog = mesh_element(3,e) micro = mesh_element(4,e) - if (homog < 1 .or. homog > material_Nhomogenization) & ! out of bounds + if (homog < 1_pInt .or. homog > material_Nhomogenization) & ! out of bounds call IO_error(154_pInt,e,0_pInt,0_pInt) - if (micro < 1 .or. micro > material_Nmicrostructure) & ! out of bounds + if (micro < 1_pInt .or. micro > material_Nmicrostructure) & ! out of bounds call IO_error(155_pInt,e,0_pInt,0_pInt) if (microstructure_elemhomo(micro)) then dGrains = homogenization_Ngrains(homog) @@ -664,7 +664,7 @@ subroutine material_populateGrains() allocate(textureOfGrain(maxval(Ngrains))) ! reserve memory for maximum case allocate(orientationOfGrain(3,maxval(Ngrains))) ! reserve memory for maximum case - if (debug_verbosity > 0) then + if (debug_verbosity > 0_pInt) then !$OMP CRITICAL (write2out) write (6,*) write (6,*) 'MATERIAL grain population' @@ -672,12 +672,12 @@ subroutine material_populateGrains() write (6,'(a32,1x,a32,1x,a6)') 'homogenization_name','microstructure_name','grain#' !$OMP END CRITICAL (write2out) endif - do homog = 1,material_Nhomogenization ! loop over homogenizations + do homog = 1_pInt,material_Nhomogenization ! loop over homogenizations dGrains = homogenization_Ngrains(homog) ! grain number per material point - do micro = 1,material_Nmicrostructure ! all pairs of homog and micro - if (Ngrains(homog,micro) > 0) then ! an active pair of homog and micro + do micro = 1_pInt,material_Nmicrostructure ! all pairs of homog and micro + if (Ngrains(homog,micro) > 0_pInt) then ! an active pair of homog and micro myNgrains = Ngrains(homog,micro) ! assign short name for total number of grains to populate - if (debug_verbosity > 0) then + if (debug_verbosity > 0_pInt) then !$OMP CRITICAL (write2out) write (6,*) write (6,'(a32,1x,a32,1x,i6)') homogenization_name(homog),microstructure_name(micro),myNgrains @@ -690,12 +690,12 @@ subroutine material_populateGrains() do hme = 1_pInt, Nelems(homog,micro) e = elemsOfHomogMicro(hme,homog,micro) ! my combination of homog and micro, only perform calculations for elements with homog, micro combinations which is indexed in cpElemsindex if (microstructure_elemhomo(micro)) then ! homogeneous distribution of grains over each element's IPs - volumeOfGrain(grain+1:grain+dGrains) = sum(mesh_ipVolume(1:FE_Nips(mesh_element(2,e)),e))/& + volumeOfGrain(grain+1_pInt:grain+dGrains) = sum(mesh_ipVolume(1:FE_Nips(mesh_element(2,e)),e))/& real(dGrains,pReal) grain = grain + dGrains ! wind forward by NgrainsPerIP else - forall (i = 1:FE_Nips(mesh_element(2,e))) & ! loop over IPs - volumeOfGrain(grain+(i-1)*dGrains+1:grain+i*dGrains) = & + forall (i = 1_pInt:FE_Nips(mesh_element(2,e))) & ! loop over IPs + volumeOfGrain(grain+(i-1)*dGrains+1_pInt:grain+i*dGrains) = & mesh_ipVolume(i,e)/dGrains ! assign IPvolume/Ngrains to all grains of IP grain = grain + FE_Nips(mesh_element(2,e)) * dGrains ! wind forward by Nips*NgrainsPerIP endif @@ -703,13 +703,13 @@ subroutine material_populateGrains() ! ---------------------------------------------------------------------------- divide myNgrains as best over constituents NgrainsOfConstituent = 0_pInt - forall (i = 1:microstructure_Nconstituents(micro)) & + forall (i = 1_pInt:microstructure_Nconstituents(micro)) & NgrainsOfConstituent(i) = nint(microstructure_fraction(i,micro) * myNgrains, pInt) ! do rounding integer conversion do while (sum(NgrainsOfConstituent) /= myNgrains) ! total grain count over constituents wrong? sgn = sign(1_pInt, myNgrains - sum(NgrainsOfConstituent)) ! direction of required change extreme = 0.0_pReal t = 0_pInt - do i = 1,microstructure_Nconstituents(micro) ! find largest deviator + do i = 1_pInt,microstructure_Nconstituents(micro) ! find largest deviator if (real(sgn,pReal)*log(NgrainsOfConstituent(i)/myNgrains/microstructure_fraction(i,micro)) > extreme) then extreme = real(sgn,pReal)*log(NgrainsOfConstituent(i)/myNgrains/microstructure_fraction(i,micro)) t = i @@ -723,21 +723,21 @@ subroutine material_populateGrains() orientationOfGrain = 0.0_pReal grain = 0_pInt ! reset microstructure grain index - do i = 1,microstructure_Nconstituents(micro) ! loop over constituents + do i = 1_pInt,microstructure_Nconstituents(micro) ! loop over constituents phaseID = microstructure_phase(i,micro) textureID = microstructure_texture(i,micro) - phaseOfGrain(grain+1:grain+NgrainsOfConstituent(i)) = phaseID ! assign resp. phase - textureOfGrain(grain+1:grain+NgrainsOfConstituent(i)) = textureID ! assign resp. texture + phaseOfGrain(grain+1_pInt:grain+NgrainsOfConstituent(i)) = phaseID ! assign resp. phase + textureOfGrain(grain+1_pInt:grain+NgrainsOfConstituent(i)) = textureID ! assign resp. texture myNorientations = ceiling(real(NgrainsOfConstituent(i),pReal)/& - real(texture_symmetry(textureID),pReal)) ! max number of unique orientations (excl. symmetry) + real(texture_symmetry(textureID),pReal),pInt) ! max number of unique orientations (excl. symmetry) constituentGrain = 0_pInt ! constituent grain index ! --------- if (texture_ODFfile(textureID) == '') then ! dealing with texture components ! --------- - do t = 1,texture_Ngauss(textureID) ! loop over Gauss components - do g = 1,int(myNorientations*texture_Gauss(5,t,textureID)) ! loop over required grain count + do t = 1_pInt,texture_Ngauss(textureID) ! loop over Gauss components + do g = 1_pInt,int(myNorientations*texture_Gauss(5,t,textureID),pInt) ! loop over required grain count orientationOfGrain(:,grain+constituentGrain+g) = & math_sampleGaussOri(texture_Gauss(1:3,t,textureID),& texture_Gauss( 4,t,textureID)) @@ -745,17 +745,17 @@ subroutine material_populateGrains() constituentGrain = constituentGrain + int(myNorientations*texture_Gauss(5,t,textureID)) enddo - do t = 1,texture_Nfiber(textureID) ! loop over fiber components - do g = 1,int(myNorientations*texture_Fiber(6,t,textureID)) ! loop over required grain count + do t = 1_pInt,texture_Nfiber(textureID) ! loop over fiber components + do g = 1_pInt,int(myNorientations*texture_Fiber(6,t,textureID),pInt) ! loop over required grain count orientationOfGrain(:,grain+constituentGrain+g) = & math_sampleFiberOri(texture_Fiber(1:2,t,textureID),& texture_Fiber(3:4,t,textureID),& texture_Fiber( 5,t,textureID)) enddo - constituentGrain = constituentGrain + int(myNorientations*texture_fiber(6,t,textureID)) + constituentGrain = constituentGrain + int(myNorientations*texture_fiber(6,t,textureID),pInt) enddo - do j = constituentGrain+1,myNorientations ! fill remainder with random + do j = constituentGrain+1_pInt,myNorientations ! fill remainder with random orientationOfGrain(:,grain+j) = math_sampleRandomOri() enddo ! --------- @@ -770,12 +770,12 @@ subroutine material_populateGrains() symExtension = texture_symmetry(textureID) - 1_pInt if (symExtension > 0_pInt) then ! sample symmetry constituentGrain = NgrainsOfConstituent(i)-myNorientations ! calc remainder of array - do j = 1,myNorientations ! loop over each "real" orientation + do j = 1_pInt,myNorientations ! loop over each "real" orientation symOrientation = math_symmetricEulers(texture_symmetry(textureID),orientationOfGrain(:,j)) ! get symmetric equivalents e = min(symExtension,constituentGrain) ! are we at end of constituent grain array? if (e > 0_pInt) then - orientationOfGrain(:,grain+myNorientations+1+(j-1)*symExtension:& - grain+myNorientations+e+(j-1)*symExtension) = & + orientationOfGrain(:,grain+myNorientations+1+(j-1_pInt)*symExtension:& + grain+myNorientations+e+(j-1_pInt)*symExtension) = & symOrientation(:,1:e) constituentGrain = constituentGrain - e ! remainder shrinks by e endif @@ -787,7 +787,7 @@ subroutine material_populateGrains() ! ---------------------------------------------------------------------------- if (.not. microstructure_elemhomo(micro)) then ! unless element homogeneous, reshuffle grains - do i=1,myNgrains-1 ! walk thru grains + do i=1_pInt,myNgrains-1_pInt ! walk thru grains call random_number(rnd) t = nint(rnd*(myNgrains-i)+i+0.5_pReal,pInt) ! select a grain in remaining list m = phaseOfGrain(t) ! exchange current with random @@ -809,7 +809,7 @@ subroutine material_populateGrains() do hme = 1_pInt, Nelems(homog,micro) e = elemsOfHomogMicro(hme,homog,micro) ! only perform calculations for elements with homog, micro combinations which is indexed in cpElemsindex if (microstructure_elemhomo(micro)) then ! homogeneous distribution of grains over each element's IPs - forall (i = 1:FE_Nips(mesh_element(2,e)), g = 1:dGrains) ! loop over IPs and grains + forall (i = 1_pInt:FE_Nips(mesh_element(2,e)), g = 1_pInt:dGrains) ! loop over IPs and grains material_volume(g,i,e) = volumeOfGrain(grain+g) material_phase(g,i,e) = phaseOfGrain(grain+g) material_texture(g,i,e) = textureOfGrain(grain+g) @@ -818,11 +818,11 @@ subroutine material_populateGrains() FEsolving_execIP(2,e) = 1_pInt ! restrict calculation to first IP only, since all other results are to be copied from this grain = grain + dGrains ! wind forward by NgrainsPerIP else - forall (i = 1:FE_Nips(mesh_element(2,e)), g = 1:dGrains) ! loop over IPs and grains - material_volume(g,i,e) = volumeOfGrain(grain+(i-1)*dGrains+g) - material_phase(g,i,e) = phaseOfGrain(grain+(i-1)*dGrains+g) - material_texture(g,i,e) = textureOfGrain(grain+(i-1)*dGrains+g) - material_EulerAngles(:,g,i,e) = orientationOfGrain(:,grain+(i-1)*dGrains+g) + forall (i = 1_pInt:FE_Nips(mesh_element(2,e)), g = 1_pInt:dGrains) ! loop over IPs and grains + material_volume(g,i,e) = volumeOfGrain(grain+(i-1_pInt)*dGrains+g) + material_phase(g,i,e) = phaseOfGrain(grain+(i-1_pInt)*dGrains+g) + material_texture(g,i,e) = textureOfGrain(grain+(i-1_pInt)*dGrains+g) + material_EulerAngles(:,g,i,e) = orientationOfGrain(:,grain+(i-1_pInt)*dGrains+g) end forall grain = grain + FE_Nips(mesh_element(2,e)) * dGrains ! wind forward by Nips*NgrainsPerIP endif diff --git a/code/math.f90 b/code/math.f90 index fcfd86947..56b72c64d 100644 --- a/code/math.f90 +++ b/code/math.f90 @@ -33,7 +33,7 @@ real(pReal), parameter :: pi = 3.14159265358979323846264338327950288419716939937510_pReal real(pReal), parameter :: inDeg = 180.0_pReal/pi real(pReal), parameter :: inRad = pi/180.0_pReal - complex(pReal), parameter :: two_pi_img = cmplx(0.0_pReal,2.0_pReal* pi, pReal) + complex(pReal), parameter :: two_pi_img = (0.0_pReal,2.0_pReal)* pi ! *** 3x3 Identity *** real(pReal), dimension(3,3), parameter :: math_I3 = & @@ -265,7 +265,7 @@ real(pReal), dimension(4,36), parameter :: math_symOperations = & integer(pInt), intent(in) :: istart,iend integer(pInt) :: d,i,j,k,x,tmp - d = size(a,1_pInt) ! number of linked data + d = int(size(a,1_pInt), pInt) ! number of linked data ! set the starting and ending points, and the pivot point i = istart @@ -3433,7 +3433,7 @@ subroutine curl_fft(res,geomdim,vec_tens,field,curl) xi(i,j,k,1:3) = real(k_s, pReal)/geomdim enddo; enddo; enddo - do k = 1_pInt, res(3); do j = 1_pInt, res(2);do i = 1_pInt, res1_red + do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res1_red do l = 1_pInt, vec_tens curl_fourier(i,j,k,l,1) = ( field_fourier(i,j,k,l,3)*xi(i,j,k,2)& -field_fourier(i,j,k,l,2)*xi(i,j,k,3) )*two_pi_img diff --git a/code/mesh.f90 b/code/mesh.f90 index a69a1e858..00816fe39 100644 --- a/code/mesh.f90 +++ b/code/mesh.f90 @@ -93,75 +93,75 @@ integer(pInt) :: hypoelasticTableStyle integer(pInt) :: initialcondTableStyle - integer(pInt), parameter :: FE_Nelemtypes = 10 - integer(pInt), parameter :: FE_maxNnodes = 8 - integer(pInt), parameter :: FE_maxNsubNodes = 56 - integer(pInt), parameter :: FE_maxNips = 27 - integer(pInt), parameter :: FE_maxNipNeighbors = 6 - integer(pInt), parameter :: FE_maxmaxNnodesAtIP = 8 - integer(pInt), parameter :: FE_NipFaceNodes = 4 + integer(pInt), parameter :: FE_Nelemtypes = 10_pInt + integer(pInt), parameter :: FE_maxNnodes = 8_pInt + integer(pInt), parameter :: FE_maxNsubNodes = 56_pInt + integer(pInt), parameter :: FE_maxNips = 27_pInt + integer(pInt), parameter :: FE_maxNipNeighbors = 6_pInt + integer(pInt), parameter :: FE_maxmaxNnodesAtIP = 8_pInt + integer(pInt), parameter :: FE_NipFaceNodes = 4_pInt integer(pInt), dimension(FE_Nelemtypes), parameter :: FE_Nnodes = & - (/8, & ! element 7 - 4, & ! element 134 - 4, & ! element 11 - 4, & ! element 27 - 4, & ! element 157 - 6, & ! element 136 - 8, & ! element 21 - 8, & ! element 117 - 8, & ! element 57 (c3d20r == c3d8 --> copy of 7) - 3 & ! element 155, 125, 128 - /) + int([8, & ! element 7 + 4, & ! element 134 + 4, & ! element 11 + 4, & ! element 27 + 4, & ! element 157 + 6, & ! element 136 + 8, & ! element 21 + 8, & ! element 117 + 8, & ! element 57 (c3d20r == c3d8 --> copy of 7) + 3 & ! element 155, 125, 128 + ],pInt) integer(pInt), dimension(FE_Nelemtypes), parameter :: FE_NoriginalNodes = & - (/8, & ! element 7 - 4, & ! element 134 - 4, & ! element 11 - 8, & ! element 27 - 4, & ! element 157 - 6, & ! element 136 - 20,& ! element 21 - 8, & ! element 117 - 20,& ! element 57 (c3d20r == c3d8 --> copy of 7) - 6 & ! element 155, 125, 128 - /) + int([8, & ! element 7 + 4, & ! element 134 + 4, & ! element 11 + 8, & ! element 27 + 4, & ! element 157 + 6, & ! element 136 + 20,& ! element 21 + 8, & ! element 117 + 20,& ! element 57 (c3d20r == c3d8 --> copy of 7) + 6 & ! element 155, 125, 128 + ],pInt) integer(pInt), dimension(FE_Nelemtypes), parameter :: FE_Nips = & - (/8, & ! element 7 - 1, & ! element 134 - 4, & ! element 11 - 9, & ! element 27 - 4, & ! element 157 - 6, & ! element 136 - 27,& ! element 21 - 1, & ! element 117 - 8, & ! element 57 (c3d20r == c3d8 --> copy of 7) - 3 & ! element 155, 125, 128 - /) + int([8, & ! element 7 + 1, & ! element 134 + 4, & ! element 11 + 9, & ! element 27 + 4, & ! element 157 + 6, & ! element 136 + 27,& ! element 21 + 1, & ! element 117 + 8, & ! element 57 (c3d20r == c3d8 --> copy of 7) + 3 & ! element 155, 125, 128 + ],pInt) integer(pInt), dimension(FE_Nelemtypes), parameter :: FE_NipNeighbors = & - (/6, & ! element 7 - 4, & ! element 134 - 4, & ! element 11 - 4, & ! element 27 - 6, & ! element 157 - 6, & ! element 136 - 6, & ! element 21 - 6, & ! element 117 - 6, & ! element 57 (c3d20r == c3d8 --> copy of 7) - 4 & ! element 155, 125, 128 - /) + int([6, & ! element 7 + 4, & ! element 134 + 4, & ! element 11 + 4, & ! element 27 + 6, & ! element 157 + 6, & ! element 136 + 6, & ! element 21 + 6, & ! element 117 + 6, & ! element 57 (c3d20r == c3d8 --> copy of 7) + 4 & ! element 155, 125, 128 + ],pInt) integer(pInt), dimension(FE_Nelemtypes), parameter :: FE_NsubNodes = & - (/19,& ! element 7 - 0, & ! element 134 - 5, & ! element 11 - 12,& ! element 27 - 0, & ! element 157 - 15,& ! element 136 - 56,& ! element 21 - 0, & ! element 117 - 19,& ! element 57 (c3d20r == c3d8 --> copy of 7) - 4 & ! element 155, 125, 128 - /) + int([19,& ! element 7 + 0, & ! element 134 + 5, & ! element 11 + 12,& ! element 27 + 0, & ! element 157 + 15,& ! element 136 + 56,& ! element 21 + 0, & ! element 117 + 19,& ! element 57 (c3d20r == c3d8 --> copy of 7) + 4 & ! element 155, 125, 128 + ],pInt) integer(pInt), dimension(FE_maxNipNeighbors,FE_Nelemtypes), parameter :: FE_NfaceNodes = & - reshape((/& + reshape(int([& 4,4,4,4,4,4, & ! element 7 3,3,3,3,0,0, & ! element 134 2,2,2,2,0,0, & ! element 11 @@ -172,21 +172,21 @@ 4,4,4,4,4,4, & ! element 117 4,4,4,4,4,4, & ! element 57 (c3d20r == c3d8 --> copy of 7) 2,2,2,0,0,0 & ! element 155, 125, 128 - /),(/FE_maxNipNeighbors,FE_Nelemtypes/)) + ],pInt),(/FE_maxNipNeighbors,FE_Nelemtypes/)) integer(pInt), dimension(FE_Nelemtypes), parameter :: FE_maxNnodesAtIP = & - (/1, & ! element 7 - 4, & ! element 134 - 1, & ! element 11 - 2, & ! element 27 - 1, & ! element 157 - 1, & ! element 136 - 4, & ! element 21 - 8, & ! element 117 - 1, & ! element 57 (c3d20r == c3d8 --> copy of 7) - 1 & ! element 155, 125, 128 - /) + int([1, & ! element 7 + 4, & ! element 134 + 1, & ! element 11 + 2, & ! element 27 + 1, & ! element 157 + 1, & ! element 136 + 4, & ! element 21 + 8, & ! element 117 + 1, & ! element 57 (c3d20r == c3d8 --> copy of 7) + 1 & ! element 155, 125, 128 + ],pInt) integer(pInt), dimension(FE_NipFaceNodes,FE_maxNipNeighbors,FE_Nelemtypes), parameter :: FE_nodeOnFace = & - reshape((/& + reshape(int([& 1,2,3,4 , & ! element 7 2,1,5,6 , & 3,2,6,7 , & @@ -247,7 +247,7 @@ 0,0,0,0 , & 0,0,0,0 , & 0,0,0,0 & - /),(/FE_NipFaceNodes,FE_maxNipNeighbors,FE_Nelemtypes/)) + ],pInt),(/FE_NipFaceNodes,FE_maxNipNeighbors,FE_Nelemtypes/)) CONTAINS ! --------------------------- @@ -263,13 +263,14 @@ subroutine mesh_init (ip,element) use DAMASK_interface + use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) use prec, only: pInt use IO, only: IO_error,IO_open_InputFile,IO_abaqus_hasNoPart use FEsolving, only: parallelExecution, FEsolving_execElem, FEsolving_execIP, calcMode, lastMode, FEmodelGeometry implicit none - integer(pInt), parameter :: fileUnit = 222 + integer(pInt), parameter :: fileUnit = 222_pInt integer(pInt) e,element,ip !$OMP CRITICAL (write2out) @@ -333,7 +334,7 @@ parallelExecution = (parallelExecution .and. (mesh_Nelems == mesh_NcpElems)) ! plus potential killer from non-local constitutive FEsolving_execElem = [ 1_pInt,mesh_NcpElems] - allocate(FEsolving_execIP(2,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(mesh_element(2,e)) allocate(calcMode(mesh_maxNips,mesh_NcpElems)) @@ -426,23 +427,23 @@ upper = int(size(lookupMap,2_pInt),pInt) ! check at bounds QUESTION is it valid to extend bounds by 1 and just do binary search w/o init check at bounds? - if (lookupMap(1,lower) == id) then - mesh_FEasCP = lookupMap(2,lower) + if (lookupMap(1_pInt,lower) == id) then + mesh_FEasCP = lookupMap(2_pInt,lower) return - elseif (lookupMap(1,upper) == id) then - mesh_FEasCP = lookupMap(2,upper) + elseif (lookupMap(1_pInt,upper) == id) then + mesh_FEasCP = lookupMap(2_pInt,upper) return endif ! binary search in between bounds - do while (upper-lower > 1) - center = (lower+upper)/2 - if (lookupMap(1,center) < id) then + do while (upper-lower > 1_pInt) + center = (lower+upper)/2_pInt + if (lookupMap(1_pInt,center) < id) then lower = center - elseif (lookupMap(1,center) > id) then + elseif (lookupMap(1_pInt,center) > id) then upper = center else - mesh_FEasCP = lookupMap(2,center) + mesh_FEasCP = lookupMap(2_pInt,center) exit endif enddo @@ -487,11 +488,11 @@ logical checkTwins minNsharedElems = mesh_maxNsharedElems + 1_pInt ! init to worst case matchingFace = 0_pInt matchingElem = 0_pInt ! intialize to "no match found" -myType = mesh_element(2,elem) ! figure elemType +myType = mesh_element(2_pInt,elem) ! figure elemType -do n = 1,FE_NfaceNodes(face,myType) ! loop over nodes on face - myFaceNodes(n) = mesh_FEasCP('node',mesh_element(4+FE_nodeOnFace(n,face,myType),elem)) ! CP id of face node - NsharedElems = mesh_sharedElem(1,myFaceNodes(n)) ! figure # shared elements for this node +do n = 1_pInt,FE_NfaceNodes(face,myType) ! loop over nodes on face + myFaceNodes(n) = mesh_FEasCP('node',mesh_element(4_pInt+FE_nodeOnFace(n,face,myType),elem)) ! CP id of face node + NsharedElems = mesh_sharedElem(1_pInt,myFaceNodes(n)) ! figure # shared elements for this node if (NsharedElems < minNsharedElems) then minNsharedElems = NsharedElems ! remember min # shared elems lonelyNode = n ! remember most lonely node @@ -501,30 +502,30 @@ enddo allocate(element_seen(minNsharedElems)) element_seen = 0_pInt -checkCandidate: do i = 1,minNsharedElems ! iterate over lonelyNode's shared elements - candidateElem = mesh_sharedElem(1+i,myFaceNodes(lonelyNode)) ! present candidate elem +checkCandidate: do i = 1_pInt,minNsharedElems ! iterate over lonelyNode's shared elements + candidateElem = mesh_sharedElem(1_pInt+i,myFaceNodes(lonelyNode)) ! present candidate elem if (all(element_seen /= candidateElem)) then ! element seen for the first time? element_seen(i) = candidateElem - candidateType = mesh_element(2,candidateElem) ! figure elemType of candidate -checkCandidateFace: do candidateFace = 1,FE_maxNipNeighbors ! check each face of candidate + candidateType = mesh_element(2_pInt,candidateElem) ! figure elemType of candidate +checkCandidateFace: do candidateFace = 1_pInt,FE_maxNipNeighbors ! check each face of candidate if (FE_NfaceNodes(candidateFace,candidateType) /= FE_NfaceNodes(face,myType) & ! incompatible face .or. (candidateElem == elem .and. candidateFace == face)) then ! this is my face cycle checkCandidateFace endif checkTwins = .false. - do n = 1,FE_NfaceNodes(candidateFace,candidateType) ! loop through nodes on face - candidateFaceNode = mesh_FEasCP('node', mesh_element(4+FE_nodeOnFace(n,candidateFace,candidateType),candidateElem)) + do n = 1_pInt,FE_NfaceNodes(candidateFace,candidateType) ! loop through nodes on face + candidateFaceNode = mesh_FEasCP('node', mesh_element(4_pInt+FE_nodeOnFace(n,candidateFace,candidateType),candidateElem)) if (all(myFaceNodes /= candidateFaceNode)) then ! candidate node does not match any of my face nodes checkTwins = .true. ! perhaps the twin nodes do match exit endif enddo if(checkTwins) then -checkCandidateFaceTwins: do dir = 1,3 - do n = 1,FE_NfaceNodes(candidateFace,candidateType) ! loop through nodes on face +checkCandidateFaceTwins: do dir = 1_pInt,3_pInt + do n = 1_pInt,FE_NfaceNodes(candidateFace,candidateType) ! loop through nodes on face candidateFaceNode = mesh_FEasCP('node', mesh_element(4+FE_nodeOnFace(n,candidateFace,candidateType),candidateElem)) if (all(myFaceNodes /= mesh_nodeTwins(dir,candidateFaceNode))) then ! node twin does not match either - if (dir == 3) then + if (dir == 3_pInt) then cycle checkCandidateFace else cycle checkCandidateFaceTwins ! try twins in next dimension @@ -564,7 +565,7 @@ endsubroutine ! fill FE_nodesAtIP with data FE_nodesAtIP(:,:FE_Nips(1),1) = & ! element 7 - reshape((/& + reshape(int([& 1, & 2, & 4, & @@ -573,20 +574,20 @@ endsubroutine 6, & 8, & 7 & - /),(/FE_maxNnodesAtIP(1),FE_Nips(1)/)) + ],pInt),(/FE_maxNnodesAtIP(1),FE_Nips(1)/)) FE_nodesAtIP(:,:FE_Nips(2),2) = & ! element 134 - reshape((/& + reshape(int([& 1,2,3,4 & - /),(/FE_maxNnodesAtIP(2),FE_Nips(2)/)) + ],pInt),(/FE_maxNnodesAtIP(2),FE_Nips(2)/)) FE_nodesAtIP(:,:FE_Nips(3),3) = & ! element 11 - reshape((/& + reshape(int([& 1, & 2, & 4, & 3 & - /),(/FE_maxNnodesAtIP(3),FE_Nips(3)/)) + ],pInt),(/FE_maxNnodesAtIP(3),FE_Nips(3)/)) FE_nodesAtIP(:,:FE_Nips(4),4) = & ! element 27 - reshape((/& + reshape(int([& 1,0, & 1,2, & 2,0, & @@ -596,25 +597,25 @@ endsubroutine 4,0, & 3,4, & 3,0 & - /),(/FE_maxNnodesAtIP(4),FE_Nips(4)/)) + ],pInt),(/FE_maxNnodesAtIP(4),FE_Nips(4)/)) FE_nodesAtIP(:,:FE_Nips(5),5) = & ! element 157 - reshape((/& + reshape(int([& 1, & 2, & 3, & 4 & - /),(/FE_maxNnodesAtIP(5),FE_Nips(5)/)) + ],pInt),(/FE_maxNnodesAtIP(5),FE_Nips(5)/)) FE_nodesAtIP(:,:FE_Nips(6),6) = & ! element 136 - reshape((/& + reshape(int([& 1, & 2, & 3, & 4, & 5, & 6 & - /),(/FE_maxNnodesAtIP(6),FE_Nips(6)/)) + ],pInt),(/FE_maxNnodesAtIP(6),FE_Nips(6)/)) FE_nodesAtIP(:,:FE_Nips(7),7) = & ! element 21 - reshape((/& + reshape(int([& 1,0, 0,0, & 1,2, 0,0, & 2,0, 0,0, & @@ -642,13 +643,13 @@ endsubroutine 8,0, 0,0, & 7,8, 0,0, & 7,0, 0,0 & - /),(/FE_maxNnodesAtIP(7),FE_Nips(7)/)) + ],pInt),(/FE_maxNnodesAtIP(7),FE_Nips(7)/)) ! FE_nodesAtIP(:,:FE_Nips(8),8) = & ! element 117 (c3d8r --> single IP per element, so no need for this mapping) ! reshape((/& ! 1,2,3,4,5,6,7,8 & ! /),(/FE_maxNnodesAtIP(8),FE_Nips(8)/)) FE_nodesAtIP(:,:FE_Nips(9),9) = & ! element 57 (c3d20r == c3d8 --> copy of 7) - reshape((/& + reshape(int([& 1, & 2, & 4, & @@ -657,13 +658,13 @@ endsubroutine 6, & 8, & 7 & - /),(/FE_maxNnodesAtIP(9),FE_Nips(9)/)) + ],pInt),(/FE_maxNnodesAtIP(9),FE_Nips(9)/)) FE_nodesAtIP(:,:FE_Nips(10),10) = & ! element 155, 125, 128 - reshape((/& + reshape(int([& 1, & 2, & 3 & - /),(/FE_maxNnodesAtIP(10),FE_Nips(10)/)) + ],pInt),(/FE_maxNnodesAtIP(10),FE_Nips(10)/)) ! *** FE_ipNeighbor *** ! is a list of the neighborhood of each IP. @@ -672,7 +673,7 @@ endsubroutine ! Negative integers denote the interface behind which the neighboring (extra-FE) IP will be located. FE_ipNeighbor(:FE_NipNeighbors(1),:FE_Nips(1),1) = & ! element 7 - reshape((/& + reshape(int([& 2,-5, 3,-2, 5,-1, & -3, 1, 4,-2, 6,-1, & 4,-5,-4, 1, 7,-1, & @@ -681,20 +682,20 @@ endsubroutine -3, 5, 8,-2,-6, 2, & 8,-5,-4, 5,-6, 3, & -3, 7,-4, 6,-6, 4 & - /),(/FE_NipNeighbors(1),FE_Nips(1)/)) + ],pInt),(/FE_NipNeighbors(1),FE_Nips(1)/)) FE_ipNeighbor(:FE_NipNeighbors(2),:FE_Nips(2),2) = & ! element 134 - reshape((/& + reshape(int([& -1,-2,-3,-4 & - /),(/FE_NipNeighbors(2),FE_Nips(2)/)) + ],pInt),(/FE_NipNeighbors(2),FE_Nips(2)/)) FE_ipNeighbor(:FE_NipNeighbors(3),:FE_Nips(3),3) = & ! element 11 - reshape((/& + reshape(int([& 2,-4, 3,-1, & -2, 1, 4,-1, & 4,-4,-3, 1, & -2, 3,-3, 2 & - /),(/FE_NipNeighbors(3),FE_Nips(3)/)) + ],pInt),(/FE_NipNeighbors(3),FE_Nips(3)/)) FE_ipNeighbor(:FE_NipNeighbors(4),:FE_Nips(4),4) = & ! element 27 - reshape((/& + reshape(int([& 2,-4, 4,-1, & 3, 1, 5,-1, & -2, 2, 6,-1, & @@ -704,25 +705,25 @@ endsubroutine 8,-4,-3, 4, & 9, 7,-3, 5, & -2, 8,-3, 6 & - /),(/FE_NipNeighbors(4),FE_Nips(4)/)) + ],pInt),(/FE_NipNeighbors(4),FE_Nips(4)/)) FE_ipNeighbor(:FE_NipNeighbors(5),:FE_Nips(5),5) = & ! element 157 - reshape((/& + reshape(int([& 2,-4, 3,-2, 4,-1, & 3,-2, 1,-3, 4,-1, & 1,-3, 2,-4, 4,-1, & 1,-3, 2,-4, 3,-2 & - /),(/FE_NipNeighbors(5),FE_Nips(5)/)) + ],pInt),(/FE_NipNeighbors(5),FE_Nips(5)/)) FE_ipNeighbor(:FE_NipNeighbors(6),:FE_Nips(6),6) = & ! element 136 - reshape((/& + reshape(int([& 2,-4, 3,-2, 4,-1, & -3, 1, 3,-2, 5,-1, & 2,-4,-3, 1, 6,-1, & 5,-4, 6,-2,-5, 1, & -3, 4, 6,-2,-5, 2, & 5,-4,-3, 4,-5, 3 & - /),(/FE_NipNeighbors(6),FE_Nips(6)/)) + ],pInt),(/FE_NipNeighbors(6),FE_Nips(6)/)) FE_ipNeighbor(:FE_NipNeighbors(7),:FE_Nips(7),7) = & ! element 21 - reshape((/& + reshape(int([& 2,-5, 4,-2,10,-1, & 3, 1, 5,-2,11,-1, & -3, 2, 6,-2,12,-1, & @@ -750,13 +751,13 @@ endsubroutine 26,-5,-4,22,-6,16, & 27,25,-4,23,-6,17, & -3,26,-4,24,-6,18 & - /),(/FE_NipNeighbors(7),FE_Nips(7)/)) + ],pInt),(/FE_NipNeighbors(7),FE_Nips(7)/)) FE_ipNeighbor(:FE_NipNeighbors(8),:FE_Nips(8),8) = & ! element 117 - reshape((/& + reshape(int([& -3,-5,-4,-2,-6,-1 & - /),(/FE_NipNeighbors(8),FE_Nips(8)/)) + ],pInt),(/FE_NipNeighbors(8),FE_Nips(8)/)) FE_ipNeighbor(:FE_NipNeighbors(9),:FE_Nips(9),9) = & ! element 57 (c3d20r == c3d8 --> copy of 7) - reshape((/& + reshape(int([& 2,-5, 3,-2, 5,-1, & -3, 1, 4,-2, 6,-1, & 4,-5,-4, 1, 7,-1, & @@ -765,13 +766,13 @@ FE_ipNeighbor(:FE_NipNeighbors(8),:FE_Nips(8),8) = & ! element 117 -3, 5, 8,-2,-6, 2, & 8,-5,-4, 5,-6, 3, & -3, 7,-4, 6,-6, 4 & - /),(/FE_NipNeighbors(9),FE_Nips(9)/)) + ],pInt),(/FE_NipNeighbors(9),FE_Nips(9)/)) FE_ipNeighbor(:FE_NipNeighbors(10),:FE_Nips(10),10) = & ! element 155, 125, 128 - reshape((/& + reshape(int([& 2,-3, 3,-1, & -2, 1, 3,-1, & 2,-3,-2, 1 & - /),(/FE_NipNeighbors(10),FE_Nips(10)/)) + ],pInt),(/FE_NipNeighbors(10),FE_Nips(10)/)) ! *** FE_subNodeParent *** ! lists the group of nodes for which the center of gravity @@ -781,7 +782,7 @@ FE_ipNeighbor(:FE_NipNeighbors(8),:FE_Nips(8),8) = & ! element 117 ! 1, 2, 3, 4, 0, 0, 0, 0 FE_subNodeParent(:FE_Nips(1),:FE_NsubNodes(1),1) = & ! element 7 - reshape((/& + reshape(int([& 1, 2, 0, 0, 0, 0, 0, 0, & 2, 3, 0, 0, 0, 0, 0, 0, & 3, 4, 0, 0, 0, 0, 0, 0, & @@ -801,18 +802,18 @@ FE_ipNeighbor(:FE_NipNeighbors(8),:FE_Nips(8),8) = & ! element 117 1, 4, 8, 5, 0, 0, 0, 0, & 5, 6, 7, 8, 0, 0, 0, 0, & 1, 2, 3, 4, 5, 6, 7, 8 & - /),(/FE_Nips(1),FE_NsubNodes(1)/)) + ],pInt),(/FE_Nips(1),FE_NsubNodes(1)/)) !FE_subNodeParent(:FE_Nips(2),:FE_NsubNodes(2),2) ! element 134 has no subnodes FE_subNodeParent(:FE_Nips(3),:FE_NsubNodes(3),3) = & ! element 11 - reshape((/& + reshape(int([& 1, 2, 0, 0, & 2, 3, 0, 0, & 3, 4, 0, 0, & 4, 1, 0, 0, & 1, 2, 3, 4 & - /),(/FE_Nips(3),FE_NsubNodes(3)/)) + ],pInt),(/FE_Nips(3),FE_NsubNodes(3)/)) FE_subNodeParent(:FE_Nips(4),:FE_NsubNodes(4),4) = & ! element 27 - reshape((/& + reshape(int([& 1, 1, 2, 0, 0, 0, 0, 0, 0, & 1, 2, 2, 0, 0, 0, 0, 0, 0, & 2, 2, 3, 0, 0, 0, 0, 0, 0, & @@ -825,13 +826,13 @@ FE_ipNeighbor(:FE_NipNeighbors(8),:FE_Nips(8),8) = & ! element 117 2, 2, 2, 2, 1, 1, 3, 3, 4, & 3, 3, 3, 3, 2, 2, 4, 4, 1, & 4, 4, 4, 4, 1, 1, 3, 3, 2 & - /),(/FE_Nips(4),FE_NsubNodes(4)/)) + ],pInt),(/FE_Nips(4),FE_NsubNodes(4)/)) !FE_subNodeParent(:FE_Nips(5),:FE_NsubNodes(5),5) = & ! element 157 ! reshape((/& ! *still to be defined* - ! /),(/FE_Nips(5),FE_NsubNodes(5)/)) + ! ],pInt),(/FE_Nips(5),FE_NsubNodes(5)/)) FE_subNodeParent(:FE_Nips(6),:FE_NsubNodes(6),6) = & ! element 136 - reshape((/& + reshape(int([& 1, 2, 0, 0, 0, 0, & 2, 3, 0, 0, 0, 0, & 3, 1, 0, 0, 0, 0, & @@ -847,9 +848,9 @@ FE_ipNeighbor(:FE_NipNeighbors(8),:FE_Nips(8),8) = & ! element 117 1, 3, 4, 6, 0, 0, & 4, 5, 6, 0, 0, 0, & 1, 2, 3, 4, 5, 6 & - /),(/FE_Nips(6),FE_NsubNodes(6)/)) + ],pInt),(/FE_Nips(6),FE_NsubNodes(6)/)) FE_subNodeParent(:FE_Nips(7),:FE_NsubNodes(7),7) = & ! element 21 - reshape((/& + reshape(int([& 1, 1, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & 1, 2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & 2, 2, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & @@ -906,10 +907,10 @@ FE_ipNeighbor(:FE_NipNeighbors(8),:FE_Nips(8),8) = & ! element 117 6, 6, 6, 6, 6, 6, 6, 6, 2, 2, 2, 2, 5, 5, 5, 5, 7, 7, 7, 7, 1, 1, 3, 3, 8, 8, 4, & 7, 7, 7, 7, 7, 7, 7, 7, 3, 3, 3, 3, 6, 6, 6, 6, 8, 8, 8, 8, 2, 2, 4, 4, 5, 5, 1, & 8, 8, 8, 8, 8, 8, 8, 8, 4, 4, 4, 4, 5, 5, 5, 5, 7, 7, 7, 7, 1, 1, 3, 3, 6, 6, 2 & - /),(/FE_Nips(7),FE_NsubNodes(7)/)) + ],pInt),(/FE_Nips(7),FE_NsubNodes(7)/)) !FE_subNodeParent(:FE_Nips(8),:FE_NsubNodes(8),8) ! element 117 has no subnodes FE_subNodeParent(:FE_Nips(9),:FE_NsubNodes(9),9) = & ! element 57 (c3d20r == c3d8 --> copy of 7) - reshape((/& + reshape(int([& 1, 2, 0, 0, 0, 0, 0, 0, & 2, 3, 0, 0, 0, 0, 0, 0, & 3, 4, 0, 0, 0, 0, 0, 0, & @@ -929,14 +930,14 @@ FE_ipNeighbor(:FE_NipNeighbors(8),:FE_Nips(8),8) = & ! element 117 1, 4, 8, 5, 0, 0, 0, 0, & 5, 6, 7, 8, 0, 0, 0, 0, & 1, 2, 3, 4, 5, 6, 7, 8 & - /),(/FE_Nips(9),FE_NsubNodes(9)/)) + ],pInt),(/FE_Nips(9),FE_NsubNodes(9)/)) FE_subNodeParent(:FE_Nips(10),:FE_NsubNodes(10),10) = & ! element 155, 125, 128 - reshape((/& + reshape(int([& 1, 2, 0, & 2, 3, 0, & 3, 1, 0, & 1, 2, 3 & - /),(/FE_Nips(10),FE_NsubNodes(10)/)) + ],pInt),(/FE_Nips(10),FE_NsubNodes(10)/)) ! *** FE_subNodeOnIPFace *** ! indicates which subnodes make up the interfaces enclosing the IP volume. @@ -950,7 +951,7 @@ FE_ipNeighbor(:FE_NipNeighbors(8),:FE_Nips(8),8) = & ! element 117 ! currently not foreseen in 2D cases. FE_subNodeOnIPFace(:FE_NipFaceNodes,:FE_NipNeighbors(1),:FE_Nips(1),1) = & ! element 7 - reshape((/& + reshape(int([& 9,21,27,22, & ! 1 1,13,25,12, & 12,25,27,21, & @@ -999,16 +1000,16 @@ FE_ipNeighbor(:FE_NipNeighbors(8),:FE_Nips(8),8) = & ! element 117 18,26,27,23, & 7,19,26,18, & 15,23,27,24 & - /),(/FE_NipFaceNodes,FE_NipNeighbors(1),FE_Nips(1)/)) + ],pInt),(/FE_NipFaceNodes,FE_NipNeighbors(1),FE_Nips(1)/)) FE_subNodeOnIPFace(:FE_NipFaceNodes,:FE_NipNeighbors(2),:FE_Nips(2),2) = & ! element 134 - reshape((/& + reshape(int([& 1, 1, 3, 2, & ! 1 1, 1, 2, 4, & 2, 2, 3, 4, & 1, 1, 4, 3 & - /),(/FE_NipFaceNodes,FE_NipNeighbors(2),FE_Nips(2)/)) + ],pInt),(/FE_NipFaceNodes,FE_NipNeighbors(2),FE_Nips(2)/)) FE_subNodeOnIPFace(:FE_NipFaceNodes,:FE_NipNeighbors(3),:FE_Nips(3),3) = & ! element 11 - reshape((/& + reshape(int([& 5, 9, 9, 5 , & ! 1 1, 8, 8, 1 , & 8, 9, 9, 8 , & @@ -1025,9 +1026,9 @@ FE_ipNeighbor(:FE_NipNeighbors(8),:FE_Nips(8),8) = & ! element 117 4, 8, 8, 4 , & 4, 7, 7, 4 , & 8, 9, 9, 8 & - /),(/FE_NipFaceNodes,FE_NipNeighbors(3),FE_Nips(3)/)) + ],pInt),(/FE_NipFaceNodes,FE_NipNeighbors(3),FE_Nips(3)/)) FE_subNodeOnIPFace(:FE_NipFaceNodes,:FE_NipNeighbors(4),:FE_Nips(4),4) = & ! element 27 - reshape((/& + reshape(int([& 9,17,17, 9 , & ! 1 1,16,16, 1 , & 16,17,17,16 , & @@ -1064,13 +1065,13 @@ FE_ipNeighbor(:FE_NipNeighbors(8),:FE_Nips(8),8) = & ! element 117 13,19,19,13 , & 3,13,13, 3 , & 12,19,19,12 & - /),(/FE_NipFaceNodes,FE_NipNeighbors(4),FE_Nips(4)/)) + ],pInt),(/FE_NipFaceNodes,FE_NipNeighbors(4),FE_Nips(4)/)) !FE_subNodeOnIPFace(:FE_NipFaceNodes,:FE_NipNeighbors(5),:FE_Nips(5),5) = & ! element 157 ! reshape((/& ! *still to be defined* ! /),(/FE_NipFaceNodes,FE_NipNeighbors(5),FE_Nips(5)/)) FE_subNodeOnIPFace(:FE_NipFaceNodes,:FE_NipNeighbors(6),:FE_Nips(6),6) = & ! element 136 - reshape((/& + reshape(int([& 7,16,21,17, & ! 1 1,10,19, 9, & 9,19,21,16, & @@ -1107,9 +1108,9 @@ FE_ipNeighbor(:FE_NipNeighbors(8),:FE_Nips(8),8) = & ! element 117 15,19,21,20, & 6,15,20,14, & 12,18,21,19 & - /),(/FE_NipFaceNodes,FE_NipNeighbors(6),FE_Nips(6)/)) + ],pInt),(/FE_NipFaceNodes,FE_NipNeighbors(6),FE_Nips(6)/)) FE_subNodeOnIPFace(:FE_NipFaceNodes,:FE_NipNeighbors(7),:FE_Nips(7),7) = & ! element 21 - reshape((/& + reshape(int([& 9,33,57,37, & ! 1 1,17,44,16, & 33,16,44,57, & @@ -1272,18 +1273,18 @@ FE_ipNeighbor(:FE_NipNeighbors(8),:FE_Nips(8),8) = & ! element 117 63,48,28,55, & 55,28, 7,29, & 63,49,23,48 & - /),(/FE_NipFaceNodes,FE_NipNeighbors(7),FE_Nips(7)/)) + ],pInt),(/FE_NipFaceNodes,FE_NipNeighbors(7),FE_Nips(7)/)) FE_subNodeOnIPFace(:FE_NipFaceNodes,:FE_NipNeighbors(8),:FE_Nips(8),8) = & ! element 117 - reshape((/& + reshape(int([& 2, 3, 7, 6, & ! 1 1, 5, 8, 4, & 3, 4, 8, 7, & 1, 2, 6, 5, & 5, 6, 7, 8, & 1, 4, 3, 2 & - /),(/FE_NipFaceNodes,FE_NipNeighbors(8),FE_Nips(8)/)) + ],pInt),(/FE_NipFaceNodes,FE_NipNeighbors(8),FE_Nips(8)/)) FE_subNodeOnIPFace(:FE_NipFaceNodes,:FE_NipNeighbors(9),:FE_Nips(9),9) = & ! element 57 (c3d20r == c3d8 --> copy of 7) - reshape((/& + reshape(int([& 9,21,27,22, & ! 1 1,13,25,12, & 12,25,27,21, & @@ -1332,9 +1333,9 @@ FE_ipNeighbor(:FE_NipNeighbors(8),:FE_Nips(8),8) = & ! element 117 18,26,27,23, & 7,19,26,18, & 15,23,27,24 & - /),(/FE_NipFaceNodes,FE_NipNeighbors(9),FE_Nips(9)/)) + ],pInt),(/FE_NipFaceNodes,FE_NipNeighbors(9),FE_Nips(9)/)) FE_subNodeOnIPFace(:FE_NipFaceNodes,:FE_NipNeighbors(10),:FE_Nips(10),10) = & ! element 155, 125, 128 - reshape((/& + reshape(int([& 4, 7, 7, 4 , & ! 1 1, 6, 6, 1 , & 6, 7, 7, 6 , & @@ -1347,7 +1348,7 @@ FE_ipNeighbor(:FE_NipNeighbors(8),:FE_Nips(8),8) = & ! element 117 3, 6, 6, 3 , & 3, 5, 5, 3 , & 6, 7, 7, 6 & - /),(/FE_NipFaceNodes,FE_NipNeighbors(10),FE_Nips(10)/)) + ],pInt),(/FE_NipFaceNodes,FE_NipNeighbors(10),FE_Nips(10)/)) endsubroutine @@ -1363,7 +1364,7 @@ FE_ipNeighbor(:FE_NipNeighbors(8),:FE_Nips(8),8) = & ! element 117 use IO implicit none - integer(pInt), parameter :: maxNchunks = 6 + integer(pInt), parameter :: maxNchunks = 6_pInt integer(pInt), dimension (1+2*maxNchunks) :: myPos integer(pInt) myUnit @@ -1403,7 +1404,7 @@ implicit none integer(pInt), intent(in) :: myUnit -integer(pInt), parameter :: maxNchunks = 5 +integer(pInt), parameter :: maxNchunks = 5_pInt integer(pInt), dimension (1+2*maxNchunks) :: myPos integer(pInt) chunk, Nchunks character(len=300) line, keyword, damaskOption, v @@ -1467,7 +1468,7 @@ enddo use IO implicit none - integer(pInt), parameter :: maxNchunks = 7 + integer(pInt), parameter :: maxNchunks = 7_pInt integer(pInt), dimension (1+2*maxNchunks) :: myPos integer(pInt) a,b,c,i,j,headerLength @@ -1488,11 +1489,11 @@ enddo endif rewind(myUnit) - do i = 1, headerLength + do i = 1_pInt, headerLength read(myUnit,'(a1024)') line myPos = IO_stringPos(line,maxNchunks) if ( IO_lc(IO_StringValue(line,myPos,1_pInt)) == 'resolution') then - do j = 2,6,2 + do j = 2_pInt,6_pInt,2_pInt select case (IO_lc(IO_stringValue(line,myPos,j))) case('a') a = IO_intValue(line,myPos,j+1_pInt) @@ -1503,7 +1504,7 @@ enddo end select enddo mesh_Nelems = a * b * c - mesh_Nnodes = (1 + a)*(1 + b)*(1 + c) + mesh_Nnodes = (1_pInt + a)*(1_pInt + b)*(1_pInt + c) endif enddo @@ -1520,7 +1521,7 @@ enddo use IO implicit none - integer(pInt), parameter :: maxNchunks = 4 + integer(pInt), parameter :: maxNchunks = 4_pInt integer(pInt), dimension (1+2*maxNchunks) :: myPos integer(pInt) myUnit @@ -1556,15 +1557,15 @@ enddo use IO implicit none - integer(pInt), parameter :: maxNchunks = 2 + integer(pInt), parameter :: maxNchunks = 2_pInt integer(pInt), dimension (1+2*maxNchunks) :: myPos character(len=300) line integer(pInt) myUnit logical inPart - mesh_Nnodes = 0 - mesh_Nelems = 0 + mesh_Nnodes = 0_pInt + mesh_Nelems = 0_pInt 610 FORMAT(A300) @@ -1599,8 +1600,8 @@ enddo endif enddo -620 if (mesh_Nnodes < 2) call IO_error(error_ID=900_pInt) - if (mesh_Nelems == 0) call IO_error(error_ID=901_pInt) +620 if (mesh_Nnodes < 2_pInt) call IO_error(error_ID=900_pInt) + if (mesh_Nelems == 0_pInt) call IO_error(error_ID=901_pInt) endsubroutine @@ -1616,7 +1617,7 @@ enddo use IO implicit none - integer(pInt), parameter :: maxNchunks = 2 + integer(pInt), parameter :: maxNchunks = 2_pInt integer(pInt), dimension (1+2*maxNchunks) :: myPos integer(pInt) myUnit @@ -1654,7 +1655,7 @@ enddo use IO implicit none - integer(pInt), parameter :: maxNchunks = 2 + integer(pInt), parameter :: maxNchunks = 2_pInt integer(pInt), dimension (1+2*maxNchunks) :: myPos character(len=300) line @@ -1696,8 +1697,8 @@ enddo use IO implicit none - integer(pInt), parameter :: maxNchunks = 2 - integer(pInt), dimension (1+2*maxNchunks) :: myPos + integer(pInt), parameter :: maxNchunks = 2_pInt + integer(pInt), dimension (1_pInt+2_pInt*maxNchunks) :: myPos character(len=300) line integer(pInt) myUnit @@ -1722,7 +1723,7 @@ enddo mesh_Nmaterials = mesh_Nmaterials + 1_pInt enddo -620 if (mesh_Nmaterials == 0) call IO_error(error_ID=903_pInt) +620 if (mesh_Nmaterials == 0_pInt) call IO_error(error_ID=903_pInt) endsubroutine @@ -1753,7 +1754,7 @@ enddo use IO implicit none - integer(pInt), parameter :: maxNchunks = 1 + integer(pInt), parameter :: maxNchunks = 1_pInt integer(pInt), dimension (1+2*maxNchunks) :: myPos integer(pInt) myUnit,i @@ -1769,7 +1770,7 @@ enddo myPos = IO_stringPos(line,maxNchunks) if ( IO_lc(IO_stringValue(line,myPos,1_pInt)) == 'hypoelastic') then - do i=1,3+hypoelasticTableStyle ! Skip 3 or 4 lines + do i=1_pInt,3_pInt+hypoelasticTableStyle ! Skip 3 or 4 lines read (myUnit,610,END=620) line enddo mesh_NcpElems = mesh_NcpElems + IO_countContinousIntValues(myUnit) @@ -1791,7 +1792,7 @@ enddo use IO implicit none - integer(pInt), parameter :: maxNchunks = 2 + integer(pInt), parameter :: maxNchunks = 2_pInt integer(pInt), dimension (1+2*maxNchunks) :: myPos character(len=300) line @@ -1799,7 +1800,7 @@ enddo logical materialFound character (len=64) materialName,elemSetName - mesh_NcpElems = 0 + mesh_NcpElems = 0_pInt 610 FORMAT(A300) @@ -1813,10 +1814,10 @@ enddo materialFound = materialName /= '' ! valid name? case('*user') if (IO_lc(IO_StringValue(line,myPos,2_pInt)) == 'material' .and. materialFound) then - do i = 1,mesh_Nmaterials ! look thru material names + do i = 1_pInt,mesh_Nmaterials ! look thru material names if (materialName == mesh_nameMaterial(i)) then ! found one elemSetName = mesh_mapMaterial(i) ! take corresponding elemSet - do k = 1,mesh_NelemSets ! look thru all elemSet definitions + do k = 1_pInt,mesh_NelemSets ! look thru all elemSet definitions if (elemSetName == mesh_nameElemSet(k)) & ! matched? mesh_NcpElems = mesh_NcpElems + mesh_mapElemSet(1,k) ! add those elem count enddo @@ -1827,7 +1828,7 @@ enddo endselect enddo -620 if (mesh_NcpElems == 0) call IO_error(error_ID=906_pInt) +620 if (mesh_NcpElems == 0_pInt) call IO_error(error_ID=906_pInt) endsubroutine @@ -1844,7 +1845,7 @@ enddo implicit none - integer(pInt), parameter :: maxNchunks = 4 + integer(pInt), parameter :: maxNchunks = 4_pInt integer(pInt), dimension (1+2*maxNchunks) :: myPos character(len=300) line @@ -1862,7 +1863,7 @@ enddo myPos = IO_stringPos(line,maxNchunks) if( (IO_lc(IO_stringValue(line,myPos,1_pInt)) == 'define' ) .and. & (IO_lc(IO_stringValue(line,myPos,2_pInt)) == 'element' ) ) then - elemSet = elemSet+1 + elemSet = elemSet+1_pInt mesh_nameElemSet(elemSet) = trim(IO_stringValue(line,myPos,4_pInt)) mesh_mapElemSet(:,elemSet) = IO_continousIntValues(myUnit,mesh_maxNelemInSet,mesh_nameElemSet,mesh_mapElemSet,mesh_NelemSets) endif @@ -1933,11 +1934,11 @@ enddo use IO implicit none - integer(pInt), parameter :: maxNchunks = 20 - integer(pInt), dimension (1+2*maxNchunks) :: myPos + integer(pInt), parameter :: maxNchunks = 20_pInt + integer(pInt), dimension (1_pInt+2_pInt*maxNchunks) :: myPos character(len=300) line - integer(pInt) myUnit,i,count + integer(pInt) myUnit,i,c logical inPart character(len=64) elemSetName,materialName @@ -1946,7 +1947,7 @@ enddo 610 FORMAT(A300) - count = 0 + c = 0_pInt inPart = .false. rewind(myUnit) do @@ -1963,7 +1964,7 @@ enddo elemSetName = '' materialName = '' - do i = 3,myPos(1_pInt) + do i = 3_pInt,myPos(1_pInt) if (IO_extractValue(IO_lc(IO_stringValue(line,myPos,i)),'elset') /= '') & elemSetName = trim(IO_extractValue(IO_lc(IO_stringValue(line,myPos,i)),'elset')) if (IO_extractValue(IO_lc(IO_stringValue(line,myPos,i)),'material') /= '') & @@ -1971,15 +1972,15 @@ enddo enddo if (elemSetName /= '' .and. materialName /= '') then - count = count + 1_pInt - mesh_nameMaterial(count) = materialName ! name of material used for this section - mesh_mapMaterial(count) = elemSetName ! mapped to respective element set + c = c + 1_pInt + mesh_nameMaterial(c) = materialName ! name of material used for this section + mesh_mapMaterial(c) = elemSetName ! mapped to respective element set endif endif enddo -620 if (count==0) call IO_error(error_ID=905_pInt) - do i=1,count +620 if (c==0_pInt) call IO_error(error_ID=905_pInt) + do i=1_pInt,c ! write(6,*)'name of materials: ',i,mesh_nameMaterial(i) ! write(6,*)'name of elemSets: ',i,mesh_mapMaterial(i) if (mesh_nameMaterial(i)=='' .or. mesh_mapMaterial(i)=='') call IO_error(error_ID=905_pInt) @@ -2003,7 +2004,7 @@ enddo allocate (mesh_mapFEtoCPnode(2_pInt,mesh_Nnodes)) ; mesh_mapFEtoCPnode = 0_pInt - forall (i = 1:mesh_Nnodes) & + forall (i = 1_pInt:mesh_Nnodes) & mesh_mapFEtoCPnode(1:2,i) = i endsubroutine @@ -2042,7 +2043,7 @@ enddo myPos = IO_stringPos(line,maxNchunks) if( IO_lc(IO_stringValue(line,myPos,1_pInt)) == 'coordinates' ) then read (myUnit,610,END=650) line ! skip crap line - do i = 1,mesh_Nnodes + do i = 1_pInt,mesh_Nnodes read (myUnit,610,END=650) line mesh_mapFEtoCPnode(1_pInt,i) = IO_fixedIntValue (line,[ 0_pInt,10_pInt],1_pInt) mesh_mapFEtoCPnode(2_pInt,i) = i @@ -2070,11 +2071,11 @@ enddo implicit none - integer(pInt), parameter :: maxNchunks = 2 - integer(pInt), dimension (1+2*maxNchunks) :: myPos + integer(pInt), parameter :: maxNchunks = 2_pInt + integer(pInt), dimension (1_pInt+2_pInt*maxNchunks) :: myPos character(len=300) line - integer(pInt) myUnit,i,count,cpNode + integer(pInt) myUnit,i,c,cpNode logical inPart allocate (mesh_mapFEtoCPnode(2_pInt,mesh_Nnodes)) ; mesh_mapFEtoCPnode = 0_pInt @@ -2098,11 +2099,11 @@ enddo IO_lc(IO_stringValue(line,myPos,2_pInt)) /= 'file' .and. & IO_lc(IO_stringValue(line,myPos,2_pInt)) /= 'response' ) & ) then - count = IO_countDataLines(myUnit) - do i = 1,count + c = IO_countDataLines(myUnit) + do i = 1_pInt,c backspace(myUnit) enddo - do i = 1,count + do i = 1_pInt,c read (myUnit,610,END=650) line myPos = IO_stringPos(line,maxNchunks) cpNode = cpNode + 1_pInt @@ -2133,7 +2134,7 @@ enddo allocate (mesh_mapFEtoCPelem(2_pInt,mesh_NcpElems)) ; mesh_mapFEtoCPelem = 0_pInt - forall (i = 1:mesh_NcpElems) & + forall (i = 1_pInt:mesh_NcpElems) & mesh_mapFEtoCPelem(1:2,i) = i endsubroutine @@ -2153,11 +2154,11 @@ enddo implicit none - integer(pInt), parameter :: maxNchunks = 1 - integer(pInt), dimension (1+2*maxNchunks) :: myPos + integer(pInt), parameter :: maxNchunks = 1_pInt + integer(pInt), dimension (1_pInt+2_pInt*maxNchunks) :: myPos character(len=300) line - integer(pInt), dimension (1+mesh_NcpElems) :: contInts + integer(pInt), dimension (1_pInt+mesh_NcpElems) :: contInts integer(pInt) myUnit,i,cpElem allocate (mesh_mapFEtoCPelem(2,mesh_NcpElems)) ; mesh_mapFEtoCPelem = 0_pInt @@ -2170,13 +2171,13 @@ enddo read (myUnit,610,END=660) line myPos = IO_stringPos(line,maxNchunks) if( IO_lc(IO_stringValue(line,myPos,1_pInt)) == 'hypoelastic' ) then - do i=1,3+hypoelasticTableStyle ! skip three (or four if new table style!) lines + do i=1_pInt,3_pInt+hypoelasticTableStyle ! skip three (or four if new table style!) lines read (myUnit,610,END=660) line enddo contInts = IO_continousIntValues(myUnit,mesh_NcpElems,mesh_nameElemSet,mesh_mapElemSet,mesh_NelemSets) - do i = 1,contInts(1) - cpElem = cpElem+1 - mesh_mapFEtoCPelem(1,cpElem) = contInts(1+i) + do i = 1_pInt,contInts(1) + cpElem = cpElem+1_pInt + mesh_mapFEtoCPelem(1,cpElem) = contInts(1_pInt+i) mesh_mapFEtoCPelem(2,cpElem) = cpElem enddo endif @@ -2200,8 +2201,8 @@ enddo implicit none - integer(pInt), parameter :: maxNchunks = 2 - integer(pInt), dimension (1+2*maxNchunks) :: myPos + integer(pInt), parameter :: maxNchunks = 2_pInt + integer(pInt), dimension (1_pInt+2_pInt*maxNchunks) :: myPos character(len=300) line integer(pInt) myUnit,i,j,k,cpElem @@ -2223,14 +2224,14 @@ enddo materialFound = materialName /= '' ! valid name? case('*user') if (IO_lc(IO_stringValue(line,myPos,2_pInt)) == 'material' .and. materialFound) then - do i = 1,mesh_Nmaterials ! look thru material names + do i = 1_pInt,mesh_Nmaterials ! look thru material names if (materialName == mesh_nameMaterial(i)) then ! found one elemSetName = mesh_mapMaterial(i) ! take corresponding elemSet - do k = 1,mesh_NelemSets ! look thru all elemSet definitions + do k = 1_pInt,mesh_NelemSets ! look thru all elemSet definitions if (elemSetName == mesh_nameElemSet(k)) then ! matched? - do j = 1,mesh_mapElemSet(1,k) + do j = 1_pInt,mesh_mapElemSet(1,k) cpElem = cpElem + 1_pInt - mesh_mapFEtoCPelem(1,cpElem) = mesh_mapElemSet(1+j,k) ! store FE id + mesh_mapFEtoCPelem(1,cpElem) = mesh_mapElemSet(1_pInt+j,k) ! store FE id mesh_mapFEtoCPelem(2,cpElem) = cpElem ! store our id enddo endif @@ -2284,8 +2285,8 @@ subroutine mesh_marc_count_cpSizes (myUnit) use IO implicit none - integer(pInt), parameter :: maxNchunks = 2 - integer(pInt), dimension (1+2*maxNchunks) :: myPos + integer(pInt), parameter :: maxNchunks = 2_pInt + integer(pInt), dimension (1_pInt+2_pInt*maxNchunks) :: myPos character(len=300) line integer(pInt) myUnit,i,t,e @@ -2302,11 +2303,11 @@ subroutine mesh_marc_count_cpSizes (myUnit) myPos = IO_stringPos(line,maxNchunks) if( IO_lc(IO_stringValue(line,myPos,1_pInt)) == 'connectivity' ) then read (myUnit,610,END=630) line ! Garbage line - do i=1,mesh_Nelems ! read all elements + do i=1_pInt,mesh_Nelems ! read all elements read (myUnit,610,END=630) line myPos = IO_stringPos(line,maxNchunks) ! limit to id and type e = mesh_FEasCP('elem',IO_intValue(line,myPos,1_pInt)) - if (e /= 0) then + if (e /= 0_pInt) then t = FE_mapElemtype(IO_stringValue(line,myPos,2_pInt)) mesh_maxNnodes = max(mesh_maxNnodes,FE_Nnodes(t)) mesh_maxNips = max(mesh_maxNips,FE_Nips(t)) @@ -2334,11 +2335,11 @@ subroutine mesh_marc_count_cpSizes (myUnit) use IO implicit none - integer(pInt), parameter :: maxNchunks = 2 - integer(pInt), dimension (1+2*maxNchunks) :: myPos + integer(pInt), parameter :: maxNchunks = 2_pInt + integer(pInt), dimension (1_pInt+2_pInt*maxNchunks) :: myPos character(len=300) line - integer(pInt) myUnit,i,count,t + integer(pInt) myUnit,i,c,t logical inPart mesh_maxNnodes = 0_pInt @@ -2364,15 +2365,15 @@ subroutine mesh_marc_count_cpSizes (myUnit) IO_lc(IO_stringValue(line,myPos,2_pInt)) /= 'response' ) & ) then t = FE_mapElemtype(IO_extractValue(IO_lc(IO_stringValue(line,myPos,2_pInt)),'type')) ! remember elem type - if (t==0) call IO_error(error_ID=910_pInt,ext_msg='mesh_abaqus_count_cpSizes') - count = IO_countDataLines(myUnit) - do i = 1,count + if (t==0_pInt) call IO_error(error_ID=910_pInt,ext_msg='mesh_abaqus_count_cpSizes') + c = IO_countDataLines(myUnit) + do i = 1_pInt,c backspace(myUnit) enddo - do i = 1,count + do i = 1_pInt,c read (myUnit,610,END=620) line myPos = IO_stringPos(line,maxNchunks) ! limit to 64 nodes max - if (mesh_FEasCP('elem',IO_intValue(line,myPos,1_pInt)) /= 0) then ! disregard non CP elems + if (mesh_FEasCP('elem',IO_intValue(line,myPos,1_pInt)) /= 0_pInt) then ! disregard non CP elems mesh_maxNnodes = max(mesh_maxNnodes,FE_Nnodes(t)) mesh_maxNips = max(mesh_maxNips,FE_Nips(t)) mesh_maxNipNeighbors = max(mesh_maxNipNeighbors,FE_NipNeighbors(t)) @@ -2397,8 +2398,8 @@ subroutine mesh_marc_count_cpSizes (myUnit) use IO implicit none - integer(pInt), parameter :: maxNchunks = 7 - integer(pInt), dimension (1+2*maxNchunks) :: myPos + integer(pInt), parameter :: maxNchunks = 7_pInt + integer(pInt), dimension (1_pInt+2_pInt*maxNchunks) :: myPos integer(pInt) a,b,c,n,i,j,headerLength real(pReal) x,y,z logical gotResolution,gotDimension @@ -2429,13 +2430,13 @@ subroutine mesh_marc_count_cpSizes (myUnit) endif rewind(myUnit) - do i = 1, headerLength + do i = 1_pInt, headerLength read(myUnit,'(a1024)') line myPos = IO_stringPos(line,maxNchunks) select case ( IO_lc(IO_StringValue(line,myPos,1_pInt)) ) case ('dimension') gotDimension = .true. - do j = 2,6,2 + do j = 2_pInt,6_pInt,2_pInt select case (IO_lc(IO_stringValue(line,myPos,j))) case('x') x = IO_floatValue(line,myPos,j+1_pInt) @@ -2447,7 +2448,7 @@ subroutine mesh_marc_count_cpSizes (myUnit) enddo case ('resolution') gotResolution = .true. - do j = 2,6,2 + do j = 2_pInt,6_pInt,2_pInt select case (IO_lc(IO_stringValue(line,myPos,j))) case('a') a = 1_pInt + IO_intValue(line,myPos,j+1_pInt) @@ -2463,13 +2464,13 @@ subroutine mesh_marc_count_cpSizes (myUnit) ! --- sanity checks --- if ((.not. gotDimension) .or. (.not. gotResolution)) call IO_error(error_ID=842_pInt) - if ((a < 1) .or. (b < 1) .or. (c < 0)) call IO_error(error_ID=843_pInt) ! 1_pInt is already added + if ((a < 1_pInt) .or. (b < 1_pInt) .or. (c < 0_pInt)) call IO_error(error_ID=843_pInt) ! 1_pInt is already added if ((x <= 0.0_pReal) .or. (y <= 0.0_pReal) .or. (z <= 0.0_pReal)) call IO_error(error_ID=844_pInt) - forall (n = 0:mesh_Nnodes-1) - mesh_node0(1,n+1) = x * dble(mod(n,a)) / dble(a-1_pInt) - mesh_node0(2,n+1) = y * dble(mod(n/a,b)) / dble(b-1_pInt) - mesh_node0(3,n+1) = z * dble(mod(n/a/b,c)) / dble(c-1_pInt) + forall (n = 0_pInt:mesh_Nnodes-1_pInt) + mesh_node0(1,n+1_pInt) = x * real(mod(n,a),pReal) / real(a-1_pInt,pReal) + mesh_node0(2,n+1_pInt) = y * real(mod(n/a,b),pReal) / real(b-1_pInt,pReal) + mesh_node0(3,n+1_pInt) = z * real(mod(n/a/b,c),pReal) / real(c-1_pInt,pReal) end forall mesh_node = mesh_node0 !why? @@ -2489,9 +2490,9 @@ subroutine mesh_marc_count_cpSizes (myUnit) use IO implicit none - integer(pInt), dimension(5), parameter :: node_ends = (/0,10,30,50,70/) - integer(pInt), parameter :: maxNchunks = 1 - integer(pInt), dimension (1+2*maxNchunks) :: myPos + integer(pInt), dimension(5), parameter :: node_ends = int([0,10,30,50,70],pInt) + integer(pInt), parameter :: maxNchunks = 1_pInt + integer(pInt), dimension (1_pInt+2_pInt*maxNchunks) :: myPos character(len=300) line integer(pInt) myUnit,i,j,m @@ -2507,7 +2508,7 @@ subroutine mesh_marc_count_cpSizes (myUnit) myPos = IO_stringPos(line,maxNchunks) if( IO_lc(IO_stringValue(line,myPos,1_pInt)) == 'coordinates' ) then read (myUnit,610,END=670) line ! skip crap line - do i=1,mesh_Nnodes + do i=1_pInt,mesh_Nnodes read (myUnit,610,END=670) line m = mesh_FEasCP('node',IO_fixedIntValue(line,node_ends,1_pInt)) forall (j = 1_pInt:3_pInt) mesh_node0(j,m) = IO_fixedNoEFloatValue(line,node_ends,j+1_pInt) @@ -2533,11 +2534,11 @@ subroutine mesh_marc_count_cpSizes (myUnit) use IO implicit none - integer(pInt), parameter :: maxNchunks = 4 - integer(pInt), dimension (1+2*maxNchunks) :: myPos + integer(pInt), parameter :: maxNchunks = 4_pInt + integer(pInt), dimension (1_pInt+2_pInt*maxNchunks) :: myPos character(len=300) line - integer(pInt) myUnit,i,j,m,count + integer(pInt) myUnit,i,j,m,c logical inPart allocate ( mesh_node0 (3,mesh_Nnodes) ); mesh_node0 = 0.0_pReal @@ -2561,15 +2562,15 @@ subroutine mesh_marc_count_cpSizes (myUnit) IO_lc(IO_stringValue(line,myPos,2_pInt)) /= 'file' .and. & IO_lc(IO_stringValue(line,myPos,2_pInt)) /= 'response' ) & ) then - count = IO_countDataLines(myUnit) ! how many nodes are defined here? - do i = 1,count + c = IO_countDataLines(myUnit) ! how many nodes are defined here? + do i = 1_pInt,c backspace(myUnit) ! rewind to first entry enddo - do i = 1,count + do i = 1_pInt,c read (myUnit,610,END=670) line myPos = IO_stringPos(line,maxNchunks) m = mesh_FEasCP('node',IO_intValue(line,myPos,1_pInt)) - forall (j=1:3) mesh_node0(j,m) = IO_floatValue(line,myPos,j+1_pInt) + forall (j=1_pInt:3_pInt) mesh_node0(j,m) = IO_floatValue(line,myPos,j+1_pInt) enddo endif enddo @@ -2592,8 +2593,8 @@ subroutine mesh_marc_count_cpSizes (myUnit) use IO implicit none - integer(pInt), parameter :: maxNchunks = 7 - integer(pInt), dimension (1+2*maxNchunks) :: myPos + integer(pInt), parameter :: maxNchunks = 7_pInt + integer(pInt), dimension (1_pInt+2_pInt*maxNchunks) :: myPos integer(pInt) a,b,c,e,i,j,homog,headerLength,maxIntCount integer(pInt), dimension(:), allocatable :: microstructures integer(pInt), dimension(1,1) :: dummySet = 0_pInt @@ -2617,12 +2618,12 @@ subroutine mesh_marc_count_cpSizes (myUnit) endif rewind(myUnit) - do i = 1, headerLength + do i = 1_pInt, headerLength read(myUnit,'(a65536)') line myPos = IO_stringPos(line,maxNchunks) select case ( IO_lc(IO_StringValue(line,myPos,1_pInt)) ) case ('resolution') - do j = 2,6,2 + do j = 2_pInt,6_pInt,2_pInt select case (IO_lc(IO_stringValue(line,myPos,j))) case('a') a = 1_pInt + IO_intValue(line,myPos,j+1_pInt) @@ -2646,18 +2647,18 @@ subroutine mesh_marc_count_cpSizes (myUnit) enddo rewind (myUnit) - do i=1,headerLength ! skip header + do i=1_pInt,headerLength ! skip header read(myUnit,'(a65536)') line enddo -100 allocate (mesh_element (4+mesh_maxNnodes,mesh_NcpElems)) ; mesh_element = 0_pInt +100 allocate (mesh_element (4_pInt+mesh_maxNnodes,mesh_NcpElems)) ; mesh_element = 0_pInt allocate (microstructures (1_pInt+maxIntCount)) ; microstructures = 2_pInt e = 0_pInt do while (e < mesh_NcpElems .and. microstructures(1) > 0_pInt) ! fill expected number of elements, stop at end of data (or blank line!) microstructures = IO_continousIntValues(myUnit,maxIntCount,dummyName,dummySet,0_pInt) ! get affected elements - do i = 1,microstructures(1_pInt) - e = e+1 ! valid element entry + do i = 1_pInt,microstructures(1_pInt) + e = e+1_pInt ! valid element entry mesh_element( 1,e) = e ! FE id mesh_element( 2,e) = FE_mapElemtype('C3D8R') ! elem type mesh_element( 3,e) = homog ! homogenization @@ -2694,14 +2695,14 @@ subroutine mesh_marc_count_cpSizes (myUnit) use IO implicit none - integer(pInt), parameter :: maxNchunks = 66 - integer(pInt), dimension (1+2*maxNchunks) :: myPos + integer(pInt), parameter :: maxNchunks = 66_pInt + integer(pInt), dimension (1_pInt+2_pInt*maxNchunks) :: myPos character(len=300) line - integer(pInt), dimension(1+mesh_NcpElems) :: contInts + integer(pInt), dimension(1_pInt+mesh_NcpElems) :: contInts integer(pInt) myUnit,i,j,sv,val,e - allocate (mesh_element (4+mesh_maxNnodes,mesh_NcpElems)) ; mesh_element = 0_pInt + allocate (mesh_element (4_pInt+mesh_maxNnodes,mesh_NcpElems)) ; mesh_element = 0_pInt 610 FORMAT(A300) @@ -2711,15 +2712,15 @@ subroutine mesh_marc_count_cpSizes (myUnit) myPos(1:1+2*1) = IO_stringPos(line,1_pInt) if( IO_lc(IO_stringValue(line,myPos,1_pInt)) == 'connectivity' ) then read (myUnit,610,END=620) line ! Garbage line - do i = 1,mesh_Nelems + do i = 1_pInt,mesh_Nelems read (myUnit,610,END=620) line myPos = IO_stringPos(line,maxNchunks) ! limit to 64 nodes max (plus ID, type) e = mesh_FEasCP('elem',IO_intValue(line,myPos,1_pInt)) - if (e /= 0) then ! disregard non CP elems + if (e /= 0_pInt) then ! disregard non CP elems mesh_element(1,e) = IO_IntValue (line,myPos,1_pInt) ! FE id mesh_element(2,e) = FE_mapElemtype(IO_StringValue(line,myPos,2_pInt)) ! elem type - forall (j = 1:FE_Nnodes(mesh_element(2,e))) & - mesh_element(j+4,e) = IO_IntValue(line,myPos,j+2_pInt) ! copy FE ids of nodes + forall (j = 1_pInt:FE_Nnodes(mesh_element(2,e))) & + mesh_element(j+4_pInt,e) = IO_IntValue(line,myPos,j+2_pInt) ! copy FE ids of nodes call IO_skipChunks(myUnit,FE_NoriginalNodes(mesh_element(2_pInt,e))-(myPos(1_pInt)-2_pInt)) ! read on if FE_Nnodes exceeds node count present on current line endif enddo @@ -2740,19 +2741,19 @@ subroutine mesh_marc_count_cpSizes (myUnit) if( (sv == 2_pInt).or.(sv == 3_pInt) ) then ! only state vars 2 and 3 of interest read (myUnit,610,END=620) line ! read line with value of state var myPos(1:1+2*1) = IO_stringPos(line,1_pInt) - do while (scan(IO_stringValue(line,myPos,1_pInt),'+-',back=.true.)>1_pInt) ! is noEfloat value? - val = NINT(IO_fixedNoEFloatValue(line,(/0_pInt,20_pInt/),1_pInt)) ! state var's value - mesh_maxValStateVar(sv-1) = max(val,mesh_maxValStateVar(sv-1_pInt)) ! remember max val of homogenization and microstructure index - if (initialcondTableStyle == 2) then + do while (scan(IO_stringValue(line,myPos,1_pInt),'+-',back=.true.)>1) ! is noEfloat value? + val = nint(IO_fixedNoEFloatValue(line,[0_pInt,20_pInt],1_pInt),pInt) ! state var's value + mesh_maxValStateVar(sv-1_pInt) = max(val,mesh_maxValStateVar(sv-1_pInt)) ! remember max val of homogenization and microstructure index + if (initialcondTableStyle == 2_pInt) then read (myUnit,610,END=630) line ! read extra line read (myUnit,610,END=630) line ! read extra line endif contInts = IO_continousIntValues(myUnit,mesh_Nelems,mesh_nameElemSet,mesh_mapElemSet,mesh_NelemSets) ! get affected elements - do i = 1,contInts(1) + do i = 1_pInt,contInts(1) e = mesh_FEasCP('elem',contInts(1_pInt+i)) mesh_element(1_pInt+sv,e) = val enddo - if (initialcondTableStyle == 0) read (myUnit,610,END=620) line ! ignore IP range for old table style + if (initialcondTableStyle == 0_pInt) read (myUnit,610,END=620) line ! ignore IP range for old table style read (myUnit,610,END=630) line myPos(1:1+2*1) = IO_stringPos(line,1_pInt) enddo @@ -2776,15 +2777,15 @@ subroutine mesh_marc_count_cpSizes (myUnit) use IO implicit none - integer(pInt), parameter :: maxNchunks = 65 - integer(pInt), dimension (1+2*maxNchunks) :: myPos + integer(pInt), parameter :: maxNchunks = 65_pInt + integer(pInt), dimension (1_pInt+2_pInt*maxNchunks) :: myPos - integer(pInt) myUnit,i,j,k,count,e,t,homog,micro + integer(pInt) myUnit,i,j,k,c,e,t,homog,micro logical inPart,materialFound character (len=64) materialName,elemSetName character(len=300) line - allocate (mesh_element (4+mesh_maxNnodes,mesh_NcpElems)) ; mesh_element = 0_pInt + allocate (mesh_element (4_pInt+mesh_maxNnodes,mesh_NcpElems)) ; mesh_element = 0_pInt 610 FORMAT(A300) @@ -2804,19 +2805,19 @@ subroutine mesh_marc_count_cpSizes (myUnit) IO_lc(IO_stringValue(line,myPos,2_pInt)) /= 'response' ) & ) then t = FE_mapElemtype(IO_extractValue(IO_lc(IO_stringValue(line,myPos,2_pInt)),'type')) ! remember elem type - if (t==0) call IO_error(error_ID=910_pInt,ext_msg='mesh_abaqus_build_elements') - count = IO_countDataLines(myUnit) - do i = 1,count + if (t==0_pInt) call IO_error(error_ID=910_pInt,ext_msg='mesh_abaqus_build_elements') + c = IO_countDataLines(myUnit) + do i = 1_pInt,c backspace(myUnit) enddo - do i = 1,count + do i = 1_pInt,c read (myUnit,610,END=620) line myPos = IO_stringPos(line,maxNchunks) ! limit to 64 nodes max e = mesh_FEasCP('elem',IO_intValue(line,myPos,1_pInt)) - if (e /= 0) then ! disregard non CP elems + if (e /= 0_pInt) then ! disregard non CP elems mesh_element(1,e) = IO_intValue(line,myPos,1_pInt) ! FE id mesh_element(2,e) = t ! elem type - forall (j=1:FE_Nnodes(t)) & + forall (j=1_pInt:FE_Nnodes(t)) & mesh_element(4_pInt+j,e) = IO_intValue(line,myPos,1_pInt+j) ! copy FE ids of nodes to position 5: call IO_skipChunks(myUnit,FE_NoriginalNodes(t)-(myPos(1_pInt)-1_pInt)) ! read on (even multiple lines) if FE_NoriginalNodes exceeds required node count endif @@ -2840,14 +2841,14 @@ subroutine mesh_marc_count_cpSizes (myUnit) materialFound ) then read (myUnit,610,END=630) line ! read homogenization and microstructure myPos(1:1+2*2) = IO_stringPos(line,2_pInt) - homog = NINT(IO_floatValue(line,myPos,1_pInt)) - micro = NINT(IO_floatValue(line,myPos,2_pInt)) - do i = 1,mesh_Nmaterials ! look thru material names + homog = nint(IO_floatValue(line,myPos,1_pInt),pInt) + micro = nint(IO_floatValue(line,myPos,2_pInt),pInt) + do i = 1_pInt,mesh_Nmaterials ! look thru material names if (materialName == mesh_nameMaterial(i)) then ! found one elemSetName = mesh_mapMaterial(i) ! take corresponding elemSet - do k = 1,mesh_NelemSets ! look thru all elemSet definitions + do k = 1_pInt,mesh_NelemSets ! look thru all elemSet definitions if (elemSetName == mesh_nameElemSet(k)) then ! matched? - do j = 1,mesh_mapElemSet(1,k) + do j = 1_pInt,mesh_mapElemSet(1,k) e = mesh_FEasCP('elem',mesh_mapElemSet(1+j,k)) mesh_element(3,e) = homog ! store homogenization mesh_element(4,e) = micro ! store microstructure @@ -2882,7 +2883,7 @@ integer(pint) e, & ! element index t, & ! element type node, & ! CP node index j, & ! node index per element - dim, & ! dimension index + myDim, & ! dimension index nodeTwin ! node twin in the specified dimension integer(pInt), dimension (mesh_Nnodes) :: node_count integer(pInt), dimension (:), allocatable :: node_seen @@ -2892,16 +2893,16 @@ allocate(node_seen(maxval(FE_Nnodes))) node_count = 0_pInt -do e = 1,mesh_NcpElems +do e = 1_pInt,mesh_NcpElems t = mesh_element(2,e) ! get element type node_seen = 0_pInt ! reset node duplicates - do j = 1,FE_Nnodes(t) ! check each node of element + do j = 1_pInt,FE_Nnodes(t) ! check each node of element node = mesh_FEasCP('node',mesh_element(4+j,e)) ! translate to internal (consecutive) numbering if (all(node_seen /= node)) then node_count(node) = node_count(node) + 1_pInt ! if FE node not yet encountered -> count it - do dim = 1,3 ! check in each dimension... - nodeTwin = mesh_nodeTwins(dim,node) + do myDim = 1_pInt,3_pInt ! check in each dimension... + nodeTwin = mesh_nodeTwins(myDim,node) if (nodeTwin > 0_pInt) & ! if I am a twin of some node... node_count(nodeTwin) = node_count(nodeTwin) + 1_pInt ! -> count me again for the twin node enddo @@ -2910,22 +2911,22 @@ do e = 1,mesh_NcpElems enddo enddo -mesh_maxNsharedElems = maxval(node_count) ! most shared node +mesh_maxNsharedElems = int(maxval(node_count),pInt) ! most shared node allocate(mesh_sharedElem(1+mesh_maxNsharedElems,mesh_Nnodes)) mesh_sharedElem = 0_pInt -do e = 1,mesh_NcpElems +do e = 1_pInt,mesh_NcpElems t = mesh_element(2,e) node_seen = 0_pInt - do j = 1,FE_Nnodes(t) - node = mesh_FEasCP('node',mesh_element(4+j,e)) + do j = 1_pInt,FE_Nnodes(t) + node = mesh_FEasCP('node',mesh_element(4_pInt+j,e)) if (all(node_seen /= node)) then mesh_sharedElem(1,node) = mesh_sharedElem(1,node) + 1_pInt ! count for each node the connected elements - mesh_sharedElem(mesh_sharedElem(1,node)+1,node) = e ! store the respective element id - do dim = 1,3 ! check in each dimension... - nodeTwin = mesh_nodeTwins(dim,node) - if (nodeTwin > 0) then ! if i am a twin of some node... + mesh_sharedElem(mesh_sharedElem(1,node)+1_pInt,node) = e ! store the respective element id + do myDim = 1_pInt,3_pInt ! check in each dimension... + nodeTwin = mesh_nodeTwins(myDim,node) + if (nodeTwin > 0_pInt) then ! if i am a twin of some node... mesh_sharedElem(1,nodeTwin) = mesh_sharedElem(1,nodeTwin) + 1_pInt ! ...count me again for the twin mesh_sharedElem(mesh_sharedElem(1,nodeTwin)+1,nodeTwin) = e ! store the respective element id endif @@ -2976,17 +2977,17 @@ mesh_ipNeighborhood = 0_pInt linkedNodes = 0_pInt -do myElem = 1,mesh_NcpElems ! loop over cpElems +do myElem = 1_pInt,mesh_NcpElems ! loop over cpElems myType = mesh_element(2,myElem) ! get elemType - do myIP = 1,FE_Nips(myType) ! loop over IPs of elem + do myIP = 1_pInt,FE_Nips(myType) ! loop over IPs of elem - do neighbor = 1,FE_NipNeighbors(myType) ! loop over neighbors of IP + do neighbor = 1_pInt,FE_NipNeighbors(myType) ! loop over neighbors of IP neighboringIPkey = FE_ipNeighbor(neighbor,myIP,myType) !*** if the key is positive, the neighbor is inside the element !*** that means, we have already found our neighboring IP - if (neighboringIPkey > 0) then + if (neighboringIPkey > 0_pInt) then mesh_ipNeighborhood(1,neighbor,myIP,myElem) = myElem mesh_ipNeighborhood(2,neighbor,myIP,myElem) = neighboringIPkey @@ -2994,7 +2995,7 @@ do myElem = 1,mesh_NcpElems !*** if the key is negative, the neighbor resides in a neighboring element !*** that means, we have to look through the face indicated by the key and see which element is behind that face - elseif (neighboringIPkey < 0) then ! neighboring element's IP + elseif (neighboringIPkey < 0_pInt) then ! neighboring element's IP myFace = -neighboringIPkey call mesh_faceMatch(myElem, myFace, matchingElem, matchingFace) ! get face and CP elem id of face match if (matchingElem > 0_pInt) then ! found match? @@ -3012,13 +3013,13 @@ do myElem = 1,mesh_NcpElems NlinkedNodes = 0_pInt linkedNodes = 0_pInt - do a = 1,FE_maxNnodesAtIP(myType) ! figure my anchor nodes on connecting face + do a = 1_pInt,FE_maxNnodesAtIP(myType) ! figure my anchor nodes on connecting face anchor = FE_nodesAtIP(a,myIP,myType) if (anchor /= 0_pInt) then ! valid anchor node if (any(FE_nodeOnFace(:,myFace,myType) == anchor)) then ! ip anchor sits on face? NlinkedNodes = NlinkedNodes + 1_pInt linkedNodes(NlinkedNodes) = & - mesh_FEasCP('node',mesh_element(4+anchor,myElem)) ! CP id of anchor node + mesh_FEasCP('node',mesh_element(4_pInt+anchor,myElem)) ! CP id of anchor node else ! something went wrong with the linkage, since not all anchors sit on my face NlinkedNodes = 0_pInt linkedNodes = 0_pInt @@ -3031,10 +3032,10 @@ do myElem = 1,mesh_NcpElems !*** and try to find an ip with matching nodes !*** also try to match with node twins -checkCandidateIP: do candidateIP = 1,FE_Nips(neighboringType) +checkCandidateIP: do candidateIP = 1_pInt,FE_Nips(neighboringType) NmatchingNodes = 0_pInt matchingNodes = 0_pInt - do a = 1,FE_maxNnodesAtIP(neighboringType) ! check each anchor node of that ip + do a = 1_pInt,FE_maxNnodesAtIP(neighboringType) ! check each anchor node of that ip anchor = FE_nodesAtIP(a,candidateIP,neighboringType) if (anchor /= 0_pInt) then ! valid anchor node if (any(FE_nodeOnFace(:,matchingFace,neighboringType) == anchor)) then ! sits on matching face? @@ -3055,7 +3056,7 @@ checkCandidateIP: do candidateIP = 1,FE_Nips(neighboringType) !*** check "normal" nodes whether they match or not checkTwins = .false. - do a = 1,NlinkedNodes + do a = 1_pInt,NlinkedNodes if (all(matchingNodes /= linkedNodes(a))) then ! this linkedNode does not match any matchingNode checkTwins = .true. exit ! no need to search further @@ -3065,8 +3066,8 @@ checkCandidateIP: do candidateIP = 1,FE_Nips(neighboringType) !*** if no match found, then also check node twins if(checkTwins) then - dir = maxloc(abs(mesh_ipAreaNormal(1:3,neighbor,myIP,myElem)),1) ! check for twins only in direction of the surface normal - do a = 1,NlinkedNodes + dir = int(maxloc(abs(mesh_ipAreaNormal(1:3,neighbor,myIP,myElem)),1),pInt) ! check for twins only in direction of the surface normal + do a = 1_pInt,NlinkedNodes twin_of_linkedNode = mesh_nodeTwins(dir,linkedNodes(a)) if (twin_of_linkedNode == 0_pInt .or. & ! twin of linkedNode does not exist... all(matchingNodes /= twin_of_linkedNode)) then ! ... or it does not match any matchingNode @@ -3109,19 +3110,19 @@ endsubroutine endif mesh_subNodeCoord = 0.0_pReal - do e = 1,mesh_NcpElems ! loop over cpElems + do e = 1_pInt,mesh_NcpElems ! loop over cpElems t = mesh_element(2,e) ! get elemType - do n = 1,FE_Nnodes(t) - mesh_subNodeCoord(1:3,n,e) = mesh_node(1:3,mesh_FEasCP('node',mesh_element(4+n,e))) ! loop over nodes of this element type + do n = 1_pInt,FE_Nnodes(t) + mesh_subNodeCoord(1:3,n,e) = mesh_node(1:3,mesh_FEasCP('node',mesh_element(4_pInt+n,e))) ! loop over nodes of this element type enddo - do n = 1,FE_NsubNodes(t) ! now for the true subnodes - do p = 1,FE_Nips(t) ! loop through possible parent nodes - if (FE_subNodeParent(p,n,t) > 0) & ! valid parent node + do n = 1_pInt,FE_NsubNodes(t) ! now for the true subnodes + do p = 1_pInt,FE_Nips(t) ! loop through possible parent nodes + if (FE_subNodeParent(p,n,t) > 0_pInt) & ! valid parent node mesh_subNodeCoord(1:3,FE_Nnodes(t)+n,e) = mesh_subNodeCoord(1:3,FE_Nnodes(t)+n,e) & - + mesh_node(1:3,mesh_FEasCP('node',mesh_element(4+FE_subNodeParent(p,n,t),e))) ! add up parents + + mesh_node(1:3,mesh_FEasCP('node',mesh_element(4_pInt+FE_subNodeParent(p,n,t),e))) ! add up parents enddo mesh_subNodeCoord(1:3,n+FE_Nnodes(t),e) = mesh_subNodeCoord(1:3,n+FE_Nnodes(t),e) & - /real(count(FE_subNodeParent(:,n,t) > 0),pReal) + /real(count(FE_subNodeParent(:,n,t) > 0_pInt),pReal) enddo enddo @@ -3140,7 +3141,7 @@ endsubroutine implicit none integer(pInt) e,f,t,i,j,k,n - logical(pInt), dimension(mesh_maxNnodes+mesh_maxNsubNodes) :: gravityNode ! flagList to find subnodes determining center of grav + logical, dimension(mesh_maxNnodes+mesh_maxNsubNodes) :: gravityNode ! flagList to find subnodes determining center of grav real(pReal), dimension(3,mesh_maxNnodes+mesh_maxNsubNodes) :: gravityNodePos ! coordinates of subnodes determining center of grav real(pReal), dimension(3) :: centerOfGravity @@ -3148,21 +3149,21 @@ endsubroutine allocate(mesh_ipCenterOfGravity(3,mesh_maxNips,mesh_NcpElems)) endif - do e = 1,mesh_NcpElems ! loop over cpElems + do e = 1_pInt,mesh_NcpElems ! loop over cpElems t = mesh_element(2,e) ! get elemType - do i = 1,FE_Nips(t) ! loop over IPs of elem + do i = 1_pInt,FE_Nips(t) ! loop over IPs of elem gravityNode = .false. ! reset flagList gravityNodePos = 0.0_pReal ! reset coordinates - do f = 1,FE_NipNeighbors(t) ! loop over interfaces of IP - do n = 1,FE_NipFaceNodes ! loop over nodes on interface + do f = 1_pInt,FE_NipNeighbors(t) ! loop over interfaces of IP + do n = 1_pInt,FE_NipFaceNodes ! loop over nodes on interface gravityNode(FE_subNodeOnIPFace(n,f,i,t)) = .true. gravityNodePos(:,FE_subNodeOnIPFace(n,f,i,t)) = mesh_subNodeCoord(:,FE_subNodeOnIPFace(n,f,i,t),e) enddo enddo - do j = 1,mesh_maxNnodes+mesh_maxNsubNodes-1 ! walk through entire flagList except last + do j = 1_pInt,mesh_maxNnodes+mesh_maxNsubNodes-1_pInt ! walk through entire flagList except last if (gravityNode(j)) then ! valid node index - do k = j+1,mesh_maxNnodes+mesh_maxNsubNodes ! walk through remainder of list + do k = j+1_pInt,mesh_maxNnodes+mesh_maxNsubNodes ! walk through remainder of list if (gravityNode(k) .and. all(abs(gravityNodePos(:,j) - gravityNodePos(:,k)) < tol_gravityNodePos)) then ! found duplicate gravityNode(j) = .false. ! delete first instance gravityNodePos(:,j) = 0.0_pReal @@ -3192,7 +3193,7 @@ endsubroutine implicit none integer(pInt) e,f,t,i,j,n - integer(pInt), parameter :: Ntriangles = FE_NipFaceNodes-2 ! each interface is made up of this many triangles + integer(pInt), parameter :: Ntriangles = FE_NipFaceNodes-2_pInt ! each interface is made up of this many triangles real(pReal), dimension(3,FE_NipFaceNodes) :: nPos ! coordinates of nodes on IP face real(pReal), dimension(Ntriangles,FE_NipFaceNodes) :: volume ! volumes of possible tetrahedra @@ -3201,16 +3202,16 @@ endsubroutine endif mesh_ipVolume = 0.0_pReal - do e = 1,mesh_NcpElems ! loop over cpElems + do e = 1_pInt,mesh_NcpElems ! loop over cpElems t = mesh_element(2_pInt,e) ! get elemType - do i = 1,FE_Nips(t) ! loop over IPs of elem - do f = 1,FE_NipNeighbors(t) ! loop over interfaces of IP and add tetrahedra which connect to CoG - forall (n = 1:FE_NipFaceNodes) & + do i = 1_pInt,FE_Nips(t) ! loop over IPs of elem + do f = 1_pInt,FE_NipNeighbors(t) ! loop over interfaces of IP and add tetrahedra which connect to CoG + forall (n = 1_pInt:FE_NipFaceNodes) & nPos(:,n) = mesh_subNodeCoord(:,FE_subNodeOnIPFace(n,f,i,t),e) - forall (n = 1:FE_NipFaceNodes, j = 1:Ntriangles) & ! start at each interface node and build valid triangles to cover interface + forall (n = 1_pInt:FE_NipFaceNodes, j = 1_pInt:Ntriangles) & ! start at each interface node and build valid triangles to cover interface volume(j,n) = math_volTetrahedron(nPos(:,n), & ! calc volume of respective tetrahedron to CoG nPos(:,1_pInt+mod(n-1_pInt +j ,FE_NipFaceNodes)), & ! start at offset j - nPos(:,1+mod(n-1_pInt +j+1_pInt,FE_NipFaceNodes)), & ! and take j's neighbor + nPos(:,1_pInt+mod(n-1_pInt +j+1_pInt,FE_NipFaceNodes)), & ! and take j's neighbor mesh_ipCenterOfGravity(:,i,e)) mesh_ipVolume(i,e) = mesh_ipVolume(i,e) + sum(volume) ! add contribution from this interface enddo @@ -3234,28 +3235,28 @@ endsubroutine implicit none integer(pInt) e,f,t,i,j,n - integer(pInt), parameter :: Ntriangles = FE_NipFaceNodes-2 ! each interface is made up of this many triangles + integer(pInt), parameter :: Ntriangles = FE_NipFaceNodes-2_pInt ! each interface is made up of this many triangles real(pReal), dimension (3,FE_NipFaceNodes) :: nPos ! coordinates of nodes on IP face real(pReal), dimension(3,Ntriangles,FE_NipFaceNodes) :: normal real(pReal), dimension(Ntriangles,FE_NipFaceNodes) :: area allocate(mesh_ipArea(mesh_maxNipNeighbors,mesh_maxNips,mesh_NcpElems)) ; mesh_ipArea = 0.0_pReal allocate(mesh_ipAreaNormal(3_pInt,mesh_maxNipNeighbors,mesh_maxNips,mesh_NcpElems)) ; mesh_ipAreaNormal = 0.0_pReal - do e = 1,mesh_NcpElems ! loop over cpElems + do e = 1_pInt,mesh_NcpElems ! loop over cpElems t = mesh_element(2,e) ! get elemType - do i = 1,FE_Nips(t) ! loop over IPs of elem - do f = 1,FE_NipNeighbors(t) ! loop over interfaces of IP - forall (n = 1:FE_NipFaceNodes) nPos(:,n) = mesh_subNodeCoord(:,FE_subNodeOnIPFace(n,f,i,t),e) - forall (n = 1:FE_NipFaceNodes, j = 1:Ntriangles) ! start at each interface node and build valid triangles to cover interface + do i = 1_pInt,FE_Nips(t) ! loop over IPs of elem + do f = 1_pInt,FE_NipNeighbors(t) ! loop over interfaces of IP + forall (n = 1_pInt:FE_NipFaceNodes) nPos(:,n) = mesh_subNodeCoord(:,FE_subNodeOnIPFace(n,f,i,t),e) + forall (n = 1_pInt:FE_NipFaceNodes, j = 1_pInt:Ntriangles) ! start at each interface node and build valid triangles to cover interface normal(:,j,n) = math_vectorproduct(nPos(:,1_pInt+mod(n+j-1_pInt,FE_NipFaceNodes)) - nPos(:,n), & ! calc their normal vectors nPos(:,1_pInt+mod(n+j-0_pInt,FE_NipFaceNodes)) - nPos(:,n)) area(j,n) = sqrt(sum(normal(:,j,n)*normal(:,j,n))) ! and area end forall - forall (n = 1:FE_NipFaceNodes, j = 1:Ntriangles, area(j,n) > 0.0_pReal) & + forall (n = 1_pInt:FE_NipFaceNodes, j = 1_pInt:Ntriangles, area(j,n) > 0.0_pReal) & normal(1:3,j,n) = normal(1:3,j,n) / area(j,n) ! make myUnit normal mesh_ipArea(f,i,e) = sum(area) / (FE_NipFaceNodes*2.0_pReal) ! area of parallelograms instead of triangles - mesh_ipAreaNormal(:,f,i,e) = sum(sum(normal,3),2)/& ! average of all valid normals + mesh_ipAreaNormal(:,f,i,e) = sum(sum(normal,3),2_pInt)/& ! average of all valid normals real(count(area > 0.0_pReal),pReal) enddo enddo @@ -3292,7 +3293,7 @@ mesh_nodeTwins = 0_pInt tolerance = 0.001_pReal * minval(mesh_ipVolume) ** 0.333_pReal -do dir = 1,3 ! check periodicity in directions of x,y,z +do dir = 1_pInt,3_pInt ! check periodicity in directions of x,y,z if (mesh_periodicSurface(dir)) then ! only if periodicity is requested @@ -3303,13 +3304,13 @@ do dir = 1,3 ! check periodicity in direction maximumNodes = 0_pInt minCoord = minval(mesh_node0(dir,:)) maxCoord = maxval(mesh_node0(dir,:)) - do node = 1,mesh_Nnodes ! loop through all nodes and find surface nodes + do node = 1_pInt,mesh_Nnodes ! loop through all nodes and find surface nodes if (abs(mesh_node0(dir,node) - minCoord) <= tolerance) then minimumNodes(1) = minimumNodes(1) + 1_pInt - minimumNodes(minimumNodes(1)+1) = node + minimumNodes(minimumNodes(1)+1_pInt) = node elseif (abs(mesh_node0(dir,node) - maxCoord) <= tolerance) then maximumNodes(1) = maximumNodes(1) + 1_pInt - maximumNodes(maximumNodes(1)+1) = node + maximumNodes(maximumNodes(1)+1_pInt) = node endif enddo @@ -3317,11 +3318,11 @@ do dir = 1,3 ! check periodicity in direction !*** find the corresponding node on the other side with the same position in this dimension unpaired = .true. - do n1 = 1,minimumNodes(1) - minimumNode = minimumNodes(n1+1) + do n1 = 1_pInt,minimumNodes(1) + minimumNode = minimumNodes(n1+1_pInt) if (unpaired(minimumNode)) then - do n2 = 1,maximumNodes(1) - maximumNode = maximumNodes(n2+1) + do n2 = 1_pInt,maximumNodes(1) + maximumNode = maximumNodes(n2+1_pInt) distance = abs(mesh_node0(:,minimumNode) - mesh_node0(:,maximumNode)) if (sum(distance) - distance(dir) <= tolerance) then ! minimum possible distance (within tolerance) mesh_nodeTwins(dir,minimumNode) = maximumNode @@ -3366,14 +3367,14 @@ if (mesh_maxValStateVar(1) < 1_pInt) call IO_error(error_ID=170_pInt) ! no homog if (mesh_maxValStateVar(2) < 1_pInt) call IO_error(error_ID=180_pInt) ! no microstructure specified allocate (mesh_HomogMicro(mesh_maxValStateVar(1),mesh_maxValStateVar(2))); mesh_HomogMicro = 0_pInt -do e = 1,mesh_NcpElems +do e = 1_pInt,mesh_NcpElems if (mesh_element(3,e) < 1_pInt) call IO_error(error_ID=170_pInt,e=e) ! no homogenization specified if (mesh_element(4,e) < 1_pInt) call IO_error(error_ID=180_pInt,e=e) ! no microstructure specified mesh_HomogMicro(mesh_element(3,e),mesh_element(4,e)) = & - mesh_HomogMicro(mesh_element(3,e),mesh_element(4,e)) + 1 ! count combinations of homogenization and microstructure + mesh_HomogMicro(mesh_element(3,e),mesh_element(4,e)) + 1_pInt ! count combinations of homogenization and microstructure enddo -if (debug_verbosity > 0) then +if (debug_verbosity > 0_pInt) then !$OMP CRITICAL (write2out) write (6,*) write (6,*) 'Input Parser: STATISTICS' @@ -3395,7 +3396,7 @@ if (debug_verbosity > 0) then write (fmt,'(a,i32.32,a)') '(9x,a2,1x,',mesh_maxValStateVar(2),'(i8))' write (6,fmt) '+-',math_range(mesh_maxValStateVar(2)) write (fmt,'(a,i32.32,a)') '(i8,1x,a2,1x,',mesh_maxValStateVar(2),'(i8))' - do i=1,mesh_maxValStateVar(1) ! loop over all (possibly assigned) homogenizations + do i=1_pInt,mesh_maxValStateVar(1) ! loop over all (possibly assigned) homogenizations write (6,fmt) i,'| ',mesh_HomogMicro(i,:) ! loop over all (possibly assigned) microstructures enddo write(6,*) @@ -3413,13 +3414,13 @@ if (debug_verbosity > 1) then write (6,*) 'Input Parser: SUBNODE COORDINATES' write (6,*) write(6,'(a8,1x,a5,1x,a15,1x,a15,1x,a20,3(1x,a12))') 'elem','IP','IP neighbor','IPFaceNodes','subNodeOnIPFace','x','y','z' - do e = 1,mesh_NcpElems ! loop over cpElems + do e = 1_pInt,mesh_NcpElems ! loop over cpElems if (debug_selectiveDebugger .and. debug_e /= e) cycle t = mesh_element(2,e) ! get elemType - do i = 1,FE_Nips(t) ! loop over IPs of elem + do i = 1_pInt,FE_Nips(t) ! loop over IPs of elem if (debug_selectiveDebugger .and. debug_i /= i) cycle - do f = 1,FE_NipNeighbors(t) ! loop over interfaces of IP - do n = 1,FE_NipFaceNodes ! loop over nodes on interface + do f = 1_pInt,FE_NipNeighbors(t) ! loop over interfaces of IP + do n = 1_pInt,FE_NipFaceNodes ! loop over nodes on interface write(6,'(i8,1x,i5,1x,i15,1x,i15,1x,i20,3(1x,f12.8))') e,i,f,n,FE_subNodeOnIPFace(n,f,i,t),& mesh_subNodeCoord(1,FE_subNodeOnIPFace(n,f,i,t),e),& mesh_subNodeCoord(2,FE_subNodeOnIPFace(n,f,i,t),e),& @@ -3431,9 +3432,9 @@ if (debug_verbosity > 1) then write(6,*) write(6,*) 'Input Parser: IP COORDINATES' write(6,'(a8,1x,a5,3(1x,a12))') 'elem','IP','x','y','z' - do e = 1,mesh_NcpElems + do e = 1_pInt,mesh_NcpElems if (debug_selectiveDebugger .and. debug_e /= e) cycle - do i = 1,FE_Nips(mesh_element(2,e)) + do i = 1_pInt,FE_Nips(mesh_element(2,e)) if (debug_selectiveDebugger .and. debug_i /= i) cycle write (6,'(i8,1x,i5,3(1x,f12.8))') e, i, mesh_ipCenterOfGravity(:,i,e) enddo @@ -3444,12 +3445,12 @@ if (debug_verbosity > 1) then write (6,'(a13,1x,e15.8)') 'total volume', sum(mesh_ipVolume) write (6,*) write (6,'(a8,1x,a5,1x,a15,1x,a5,1x,a15,1x,a16)') 'elem','IP','volume','face','area','-- normal --' - do e = 1,mesh_NcpElems + do e = 1_pInt,mesh_NcpElems if (debug_selectiveDebugger .and. debug_e /= e) cycle - do i = 1,FE_Nips(mesh_element(2,e)) + do i = 1_pInt,FE_Nips(mesh_element(2,e)) if (debug_selectiveDebugger .and. debug_i /= i) cycle write (6,'(i8,1x,i5,1x,e15.8)') e,i,mesh_IPvolume(i,e) - do f = 1,FE_NipNeighbors(mesh_element(2,e)) + do f = 1_pInt,FE_NipNeighbors(mesh_element(2,e)) write (6,'(i33,1x,e15.8,1x,3(f6.3,1x))') f,mesh_ipArea(f,i,e),mesh_ipAreaNormal(:,f,i,e) enddo enddo @@ -3458,7 +3459,7 @@ if (debug_verbosity > 1) then write (6,*) 'Input Parser: NODE TWINS' write (6,*) write(6,'(a6,3(3x,a6))') ' node','twin_x','twin_y','twin_z' - do n = 1,mesh_Nnodes ! loop over cpNodes + do n = 1_pInt,mesh_Nnodes ! loop over cpNodes if (debug_e <= mesh_NcpElems) then if (any(mesh_element(5:,debug_e) == n)) then write(6,'(i6,3(3x,i6))') n, mesh_nodeTwins(1:3,n) @@ -3469,12 +3470,12 @@ if (debug_verbosity > 1) then write(6,*) 'Input Parser: IP NEIGHBORHOOD' write(6,*) write(6,'(a8,1x,a10,1x,a10,1x,a3,1x,a13,1x,a13)') 'elem','IP','neighbor','','elemNeighbor','ipNeighbor' - do e = 1,mesh_NcpElems ! loop over cpElems + do e = 1_pInt,mesh_NcpElems ! loop over cpElems if (debug_selectiveDebugger .and. debug_e /= e) cycle t = mesh_element(2,e) ! get elemType - do i = 1,FE_Nips(t) ! loop over IPs of elem + do i = 1_pInt,FE_Nips(t) ! loop over IPs of elem if (debug_selectiveDebugger .and. debug_i /= i) cycle - do n = 1,FE_NipNeighbors(t) ! loop over neighbors of IP + do n = 1_pInt,FE_NipNeighbors(t) ! loop over neighbors of IP write (6,'(i8,1x,i10,1x,i10,1x,a3,1x,i13,1x,i13)') e,i,n,'-->',mesh_ipNeighborhood(1,n,i,e),mesh_ipNeighborhood(2,n,i,e) enddo enddo diff --git a/code/numerics.f90 b/code/numerics.f90 index 786ff4c82..84c99a397 100644 --- a/code/numerics.f90 +++ b/code/numerics.f90 @@ -69,18 +69,18 @@ real(pReal) :: relevantStrain = 1.0e-7_pReal, & err_stress_tolrel = 0.01_pReal , & ! relative tolerance for fullfillment of stress BC, Default: 0.01 allowing deviation of 1% of maximum stress fftw_timelimit = -1.0_pReal, & ! sets the timelimit of plan creation for FFTW, see manual on www.fftw.org, Default -1.0: disable timelimit rotation_tol = 1.0e-12_pReal ! tolerance of rotation specified in loadcase, Default 1.0e-12: first guess -character(len=64) :: fftw_planner_string = 'FFTW_PATIENT' ! reads the planing-rigor flag, see manual on www.fftw.org, Default FFTW_PATIENT: use patiant planner flag -integer(pInt) :: fftw_planner_flag = -1_pInt ! conversion of fftw_planner_string to integer, basically what is usually done in the include file of fftw -logical :: memory_efficient = .true. ,& ! for fast execution (pre calculation of gamma_hat), Default .true.: do not precalculate - divergence_correction = .false. ,& ! correct divergence calculation in fourier space, Default .false.: no correction - update_gamma = .false.,& ! update gamma operator with current stiffness, Default .false.: use initial stiffness +character(len=64) :: fftw_plan_mode = 'FFTW_PATIENT' ! reads the planing-rigor flag, see manual on www.fftw.org, Default FFTW_PATIENT: use patiant planner flag +integer(pInt) :: fftw_planner_flag = -1_pInt, & ! conversion of fftw_plan_mode to integer, basically what is usually done in the include file of fftw + itmax = 20_pInt ! maximum number of iterations +logical :: memory_efficient = .true., & ! for fast execution (pre calculation of gamma_hat), Default .true.: do not precalculate + divergence_correction = .false., & ! correct divergence calculation in fourier space, Default .false.: no correction + update_gamma = .false., & ! update gamma operator with current stiffness, Default .false.: use initial stiffness simplified_algorithm = .true. ! use short algorithm without fluctuation field, Default .true.: use simplified algorithm -real(pReal) :: cut_off_value = 0.0_pReal ! percentage of frequencies to cut away, Default 0.0: use all frequencies -integer(pInt) :: itmax = 20_pInt , & ! maximum number of iterations + !* Random seeding parameters - fixedSeed = 0_pInt ! fixed seeding for pseudo-random number generator, Default 0: use random seed +integer(pInt) :: fixedSeed = 0_pInt ! fixed seeding for pseudo-random number generator, Default 0: use random seed !* OpenMP variable integer(pInt) :: DAMASK_NumThreadsInt = 0_pInt ! value stored in environment variable DAMASK_NUM_THREADS, set to zero if no OpenMP directive @@ -111,7 +111,7 @@ subroutine numerics_init() !*** local variables ***! integer(pInt), parameter :: fileunit = 300_pInt integer(pInt), parameter :: maxNchunks = 2_pInt - integer(pInt) :: gotDAMASK_NUM_THREADS = 1_pInt + integer :: gotDAMASK_NUM_THREADS = 1 integer(pInt), dimension(1+2*maxNchunks) :: positions character(len=64) :: tag character(len=1024) :: line @@ -127,9 +127,9 @@ subroutine numerics_init() !$OMP END CRITICAL (write2out) !$ call GET_ENVIRONMENT_VARIABLE(NAME='DAMASK_NUM_THREADS',VALUE=DAMASK_NumThreadsString,STATUS=gotDAMASK_NUM_THREADS) ! get environment variable DAMASK_NUM_THREADS... -!$ if(gotDAMASK_NUM_THREADS /= 0_pInt) call IO_warning(47_pInt,ext_msg=DAMASK_NumThreadsString) +!$ if(gotDAMASK_NUM_THREADS /= 0) call IO_warning(47_pInt,ext_msg=DAMASK_NumThreadsString) !$ read(DAMASK_NumThreadsString,'(i6)') DAMASK_NumThreadsInt ! ...convert it to integer... -!$ if (DAMASK_NumThreadsInt < 1) DAMASK_NumThreadsInt = 1 ! ...ensure that its at least one... +!$ if (DAMASK_NumThreadsInt < 1_pInt) DAMASK_NumThreadsInt = 1_pInt ! ...ensure that its at least one... !$ call omp_set_num_threads(DAMASK_NumThreadsInt) ! ...and use it as number of threads for parallel execution ! try to open the config file @@ -238,8 +238,8 @@ subroutine numerics_init() memory_efficient = IO_intValue(line,positions,2_pInt) > 0_pInt case ('fftw_timelimit') fftw_timelimit = IO_floatValue(line,positions,2_pInt) - case ('fftw_planner_string') - fftw_planner_string = IO_stringValue(line,positions,2_pInt) + case ('fftw_plan_mode') + fftw_plan_mode = IO_stringValue(line,positions,2_pInt) case ('rotation_tol') rotation_tol = IO_floatValue(line,positions,2_pInt) case ('divergence_correction') @@ -248,8 +248,6 @@ subroutine numerics_init() update_gamma = IO_intValue(line,positions,2_pInt) > 0_pInt case ('simplified_algorithm') simplified_algorithm = IO_intValue(line,positions,2_pInt) > 0_pInt - case ('cut_off_value') - cut_off_value = IO_floatValue(line,positions,2_pInt) !* Random seeding parameters @@ -272,7 +270,7 @@ subroutine numerics_init() endif - select case(IO_lc(fftw_planner_string)) ! setting parameters for the plan creation of FFTW. Basically a translation from fftw3.f + select case(IO_lc(fftw_plan_mode)) ! setting parameters for the plan creation of FFTW. Basically a translation from fftw3.f case('estimate','fftw_estimate') ! ordered from slow execution (but fast plan creation) to fast execution fftw_planner_flag = 64_pInt case('measure','fftw_measure') @@ -282,7 +280,7 @@ subroutine numerics_init() case('exhaustive','fftw_exhaustive') fftw_planner_flag = 8_pInt case default - call IO_warning(warning_ID=47_pInt,ext_msg=trim(IO_lc(fftw_planner_string))) + call IO_warning(warning_ID=47_pInt,ext_msg=trim(IO_lc(fftw_plan_mode))) fftw_planner_flag = 32_pInt end select @@ -341,13 +339,12 @@ subroutine numerics_init() else write(6,'(a24,1x,e8.1)') ' fftw_timelimit: ',fftw_timelimit endif - write(6,'(a24,1x,a)') ' fftw_planner_string: ',trim(fftw_planner_string) + write(6,'(a24,1x,a)') ' fftw_plan_mode: ',trim(fftw_plan_mode) write(6,'(a24,1x,i8)') ' fftw_planner_flag: ',fftw_planner_flag write(6,'(a24,1x,e8.1)') ' rotation_tol: ',rotation_tol write(6,'(a24,1x,L8,/)') ' divergence_correction: ',divergence_correction write(6,'(a24,1x,L8,/)') ' update_gamma: ',update_gamma write(6,'(a24,1x,L8,/)') ' simplified_algorithm: ',simplified_algorithm - write(6,'(a24,1x,e8.1)') ' cut_off_value: ',cut_off_value !* Random seeding parameters @@ -411,7 +408,7 @@ subroutine numerics_init() if (fixedSeed <= 0_pInt) then !$OMP CRITICAL (write2out) - write(6,'(a)') 'Random is random!' + write(6,'(a)') ' Random is random!' !$OMP END CRITICAL (write2out) endif endsubroutine diff --git a/code/prec.f90 b/code/prec.f90 index 515ac0ec7..e8ef60a08 100644 --- a/code/prec.f90 +++ b/code/prec.f90 @@ -35,13 +35,13 @@ real(pReal), parameter, public :: tol_gravityNodePos = 1.0e-100_pReal ! from http://www.hpc.unimelb.edu.au/doc/f90lrm/dfum_035.html ! copy can be found in documentation/Code/Fortran #ifdef __INTEL_COMPILER -#if __INTEL_COMPILER<12000 - real(pReal), parameter, public :: DAMASK_NaN = Z'7FF0000000000001' +#if __INTEL_COMPILER<1200 + real(pReal), parameter, public :: DAMASK_NaN = Z'7FF0000000000001' #else - real(pReal), parameter, public :: DAMASK_NaN = real(Z'7FF0000000000001', pReal) + real(pReal), parameter, public :: DAMASK_NaN = real(Z'7FF0000000000001', pReal) #endif #else -real(pReal), parameter, public :: DAMASK_NaN = real(Z'7FF0000000000001', pReal) + real(pReal), parameter, public :: DAMASK_NaN = real(Z'7FF0000000000001', pReal) #endif type :: p_vec real(pReal), dimension(:), pointer :: p diff --git a/code/prec_single.f90 b/code/prec_single.f90 index 9deb3a4e2..cc3ae07a1 100644 --- a/code/prec_single.f90 +++ b/code/prec_single.f90 @@ -35,13 +35,13 @@ real(pReal), parameter, public :: tol_gravityNodePos = 1.0e-36_pReal ! from http://www.hpc.unimelb.edu.au/doc/f90lrm/dfum_035.html ! copy can be found in documentation/Code/Fortran #ifdef __INTEL_COMPILER -#if __INTEL_COMPILER<12000 - real(pReal), parameter, public :: DAMASK_NaN = Z'Z'7F800001', pReal' +#if __INTEL_COMPILER<1200 + real(pReal), parameter, public :: DAMASK_NaN = Z'Z'7F800001', pReal' #else - real(pReal), parameter, public :: DAMASK_NaN = real(Z'7F800001', pReal) + real(pReal), parameter, public :: DAMASK_NaN = real(Z'7F800001', pReal) #endif #else -real(pReal), parameter, public :: DAMASK_NaN = real(Z'7F800001', pReal) + real(pReal), parameter, public :: DAMASK_NaN = real(Z'7F800001', pReal) #endif type :: p_vec real(pReal), dimension(:), pointer :: p