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

did a lot of polishing: 

- removed unnecessary "return" before end of subroutine or function:
- changed undetermined array length (:) to (1:3)

To prevent problems with some code analysing tools:
- "3D oneliner loops" (with ";) only for "do" and "enddo" at the same time
- removed line continuation in OMP statements

made the makefile more flexible, removed heap-arrays switch
This commit is contained in:
Martin Diehl 2011-09-13 15:54:06 +00:00
parent d235de1aa5
commit c35ea33f8e
9 changed files with 128 additions and 206 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

14
code/DAMASK_spectral.f90
View File

@ -152,7 +152,7 @@ program DAMASK_spectral
print '(a,/,a)', 'Workingdir: ',trim(getSolverWorkingDirectoryName()) print '(a,/,a)', 'Workingdir: ',trim(getSolverWorkingDirectoryName())
print '(a,/,a)', 'SolverJobName: ',trim(getSolverJobName()) print '(a,/,a)', 'SolverJobName: ',trim(getSolverJobName())
if (.not. IO_open_file(unit,path)) call IO_error(30,ext_msg = path) if (.not. IO_open_file(unit,path)) call IO_error(30,ext_msg = trim(path))
rewind(unit) rewind(unit)
do do
@ -181,7 +181,7 @@ program DAMASK_spectral
101 N_Loadcases = N_n 101 N_Loadcases = N_n
if ((N_l + N_Fdot /= N_n) .or. (N_n /= N_t)) & ! sanity check if ((N_l + N_Fdot /= N_n) .or. (N_n /= N_t)) & ! sanity check
call IO_error(31,ext_msg = path) ! error message for incomplete loadcase call IO_error(31,ext_msg = trim(path)) ! error message for incomplete loadcase
! allocate memory depending on lines in input file ! allocate memory depending on lines in input file
allocate (bc_deformation(3,3,N_Loadcases)); bc_deformation = 0.0_pReal allocate (bc_deformation(3,3,N_Loadcases)); bc_deformation = 0.0_pReal
@ -476,7 +476,7 @@ program DAMASK_spectral
time = time + timeinc time = time + timeinc
if (velGradApplied(loadcase)) & ! calculate fDot from given L and current F if (velGradApplied(loadcase)) & ! calculate fDot from given L and current F
fDot = math_mul33x33(bc_deformation(:,:,loadcase), defgradAim) fDot = math_mul33x33(bc_deformation(1:3,1:3,loadcase), defgradAim)
!winding forward of deformation aim !winding forward of deformation aim
temp33_Real = defgradAim temp33_Real = defgradAim
@ -489,11 +489,11 @@ program DAMASK_spectral
do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1) do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1)
temp33_Real = defgrad(i,j,k,:,:) temp33_Real = defgrad(i,j,k,:,:)
if (velGradApplied(loadcase)) & ! use velocity gradient to calculate new deformation gradient (if not guessing) if (velGradApplied(loadcase)) & ! use velocity gradient to calculate new deformation gradient (if not guessing)
fDot = math_mul33x33(bc_deformation(:,:,loadcase),defgradold(i,j,k,:,:)) fDot = math_mul33x33(bc_deformation(1:3,1:3,loadcase),defgradold(i,j,k,1:3,1:3))
defgrad(i,j,k,:,:) = defgrad(i,j,k,:,:) & ! decide if guessing along former trajectory or apply homogeneous addon defgrad(i,j,k,1:3,1:3) = defgrad(i,j,k,1:3,1:3) & ! decide if guessing along former trajectory or apply homogeneous addon
+ guessmode * (defgrad(i,j,k,:,:) - defgradold(i,j,k,:,:))& ! guessing... + guessmode * (defgrad(i,j,k,1:3,1:3) - defgradold(i,j,k,1:3,1:3))& ! guessing...
+ (1.0_pReal-guessmode) * mask_defgrad * fDot *timeinc ! apply the prescribed value where deformation is given if not guessing + (1.0_pReal-guessmode) * mask_defgrad * fDot *timeinc ! apply the prescribed value where deformation is given if not guessing
defgradold(i,j,k,:,:) = temp33_Real defgradold(i,j,k,1:3,1:3) = temp33_Real
enddo; enddo; enddo enddo; enddo; enddo
guessmode = 1.0_pReal ! keep guessing along former trajectory during same loadcase guessmode = 1.0_pReal ! keep guessing along former trajectory during same loadcase

2
code/DAMASK_spectral_interface.f90
View File

@ -190,7 +190,7 @@ function rectifyPath(path)
!remove ./ from path !remove ./ from path
l = len_trim(path) l = len_trim(path)
rectifyPath = path rectifyPath = path
do i = l,2,-1 do i = l,3,-1
if ( rectifyPath(i-1:i) == './' .and. rectifyPath(i-2:i-2) /= '.' ) & if ( rectifyPath(i-1:i) == './' .and. rectifyPath(i-2:i-2) /= '.' ) &
rectifyPath(i-1:l) = rectifyPath(i+1:l)//' ' rectifyPath(i-1:l) = rectifyPath(i+1:l)//' '
enddo enddo

6
code/IO.f90
View File

@ -58,7 +58,6 @@ subroutine IO_init ()
call flush(6) call flush(6)
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
return
endsubroutine endsubroutine
@ -115,7 +114,6 @@ recursive function IO_abaqus_assembleInputFile(unit1,unit2) result(createSuccess
return return
200 createSuccess =.false. 200 createSuccess =.false.
return
end function end function
@ -536,7 +534,7 @@ end function
!******************************************************************** !********************************************************************
! get tagged content of line ! get tagged content of line
!******************************************************************** !********************************************************************
pure function IO_getTag (line,openChar,closechar) pure function IO_getTag (line,openChar,closeChar)
use prec, only: pInt use prec, only: pInt
implicit none implicit none
@ -553,8 +551,6 @@ end function
if (left == verify(line,sep) .and. right > left) & ! openChar is first and closeChar occurs if (left == verify(line,sep) .and. right > left) & ! openChar is first and closeChar occurs
IO_getTag = line(left+1:right-1) IO_getTag = line(left+1:right-1)
return
endfunction endfunction

15
code/crystallite.f90
View File

@ -805,7 +805,9 @@ if(updateJaco) then
mySizeDotState = constitutive_sizeDotState(g,i,e) mySizeDotState = constitutive_sizeDotState(g,i,e)
constitutive_state(g,i,e)%p(1:mySizeState) = constitutive_state_backup(g,i,e)%p(1:mySizeState) constitutive_state(g,i,e)%p(1:mySizeState) = constitutive_state_backup(g,i,e)%p(1:mySizeState)
constitutive_dotState(g,i,e)%p(1:mySizeDotState) = constitutive_dotState_backup(g,i,e)%p(1:mySizeDotState) constitutive_dotState(g,i,e)%p(1:mySizeDotState) = constitutive_dotState_backup(g,i,e)%p(1:mySizeDotState)
enddo; enddo; enddo enddo
enddo
enddo
!OMP END PARALLEL DO !OMP END PARALLEL DO
crystallite_Temperature = Temperature_backup crystallite_Temperature = Temperature_backup
crystallite_subF = F_backup crystallite_subF = F_backup
@ -842,8 +844,10 @@ if(updateJaco) then
elseif (crystallite_requested(g,i,e) .and. .not. crystallite_converged(g,i,e)) then ! central solution did not converge elseif (crystallite_requested(g,i,e) .and. .not. crystallite_converged(g,i,e)) then ! central solution did not converge
crystallite_dPdF(1:3,1:3,1:3,1:3,g,i,e) = crystallite_fallbackdPdF(1:3,1:3,1:3,1:3,g,i,e) ! use (elastic) fallback crystallite_dPdF(1:3,1:3,1:3,1:3,g,i,e) = crystallite_fallbackdPdF(1:3,1:3,1:3,1:3,g,i,e) ! use (elastic) fallback
endif endif
enddo; enddo; enddo enddo
!OMP END PARALLEL DO enddo
enddo
!$OMP END PARALLEL DO
endif ! jacobian calculation endif ! jacobian calculation
@ -2996,8 +3000,7 @@ logical error
! --- UPDATE SOME ADDITIONAL VARIABLES THAT ARE NEEDED FOR NONLOCAL MATERIAL --- ! --- UPDATE SOME ADDITIONAL VARIABLES THAT ARE NEEDED FOR NONLOCAL MATERIAL ---
! --- we use crystallite_orientation from above, so need a seperate loop ! --- we use crystallite_orientation from above, so need a seperate loop
!$OMP PARALLEL DO PRIVATE(myPhase,myInstance,myStructure,neighboring_e,neighboring_i, & !$OMP PARALLEL DO PRIVATE(myPhase,myInstance,myStructure,neighboring_e,neighboring_i,neighboringPhase,neighboringInstance,neighboringStructure)
!$OMP & neighboringPhase,neighboringInstance,neighboringStructure)
do e = FEsolving_execElem(1),FEsolving_execElem(2) do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
myPhase = material_phase(1,i,e) ! get my phase myPhase = material_phase(1,i,e) ! get my phase
@ -3162,8 +3165,6 @@ function crystallite_postResults(&
dt, g, i, e) dt, g, i, e)
c = c + constitutive_sizePostResults(g,i,e) c = c + constitutive_sizePostResults(g,i,e)
return
endfunction endfunction

40
code/homogenization_RGC.f90
View File

@ -219,8 +219,6 @@ subroutine homogenization_RGC_init(&
! (6) Volume discrepancy, (7) Avg relaxation rate component, (8) Max relaxation rate component ! (6) Volume discrepancy, (7) Avg relaxation rate component, (8) Max relaxation rate component
enddo enddo
return
endsubroutine endsubroutine
@ -239,8 +237,6 @@ function homogenization_RGC_stateInit(myInstance)
! open(1978,file='homogenization_RGC_debugging.out',status='unknown') ! open(1978,file='homogenization_RGC_debugging.out',status='unknown')
homogenization_RGC_stateInit = 0.0_pReal homogenization_RGC_stateInit = 0.0_pReal
return
endfunction endfunction
@ -320,8 +316,6 @@ subroutine homogenization_RGC_partitionDeformation(&
enddo enddo
return
endsubroutine endsubroutine
!******************************************************************** !********************************************************************
@ -774,8 +768,6 @@ function homogenization_RGC_updateState(&
deallocate(tract,resid,jmatrix,jnverse,relax,drelax,pmatrix,smatrix,p_relax,p_resid) deallocate(tract,resid,jmatrix,jnverse,relax,drelax,pmatrix,smatrix,p_relax,p_resid)
!*** End of calculation of state update !*** End of calculation of state update
return
endfunction endfunction
!******************************************************************** !********************************************************************
@ -815,7 +807,7 @@ subroutine homogenization_RGC_averageStressAndItsTangent(&
if (debug_verbosity == 4) then if (debug_verbosity == 4) then
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
do iGrain = 1,Ngrains do iGrain = 1,Ngrains
dPdF99 = math_Plain3333to99(dPdF(:,:,:,:,iGrain)) dPdF99 = math_Plain3333to99(dPdF(1:3,1:3,1:3,1:3,iGrain))
write(6,'(x,a30,x,i3)')'Stress tangent of grain: ',iGrain write(6,'(x,a30,x,i3)')'Stress tangent of grain: ',iGrain
do i = 1,9 do i = 1,9
write(6,'(x,(e14.8,x))') (dPdF99(i,j), j = 1,9) write(6,'(x,(e14.8,x))') (dPdF99(i,j), j = 1,9)
@ -830,8 +822,6 @@ subroutine homogenization_RGC_averageStressAndItsTangent(&
avgP = sum(P,3)/dble(Ngrains) avgP = sum(P,3)/dble(Ngrains)
dAvgPdAvgF = sum(dPdF,5)/dble(Ngrains) dAvgPdAvgF = sum(dPdF,5)/dble(Ngrains)
return
endsubroutine endsubroutine
!******************************************************************** !********************************************************************
@ -858,8 +848,6 @@ function homogenization_RGC_averageTemperature(&
Ngrains = homogenization_Ngrains(mesh_element(3,el)) Ngrains = homogenization_Ngrains(mesh_element(3,el))
homogenization_RGC_averageTemperature = sum(Temperature(1:Ngrains))/dble(Ngrains) homogenization_RGC_averageTemperature = sum(Temperature(1:Ngrains))/dble(Ngrains)
return
endfunction endfunction
!******************************************************************** !********************************************************************
@ -916,8 +904,6 @@ pure function homogenization_RGC_postResults(&
end select end select
enddo enddo
return
endfunction endfunction
!******************************************************************** !********************************************************************
@ -1051,8 +1037,6 @@ subroutine homogenization_RGC_stressPenalty(&
enddo enddo
!*** End of mismatch and penalty stress tensor calculation !*** End of mismatch and penalty stress tensor calculation
return
endsubroutine endsubroutine
!******************************************************************** !********************************************************************
@ -1112,8 +1096,6 @@ subroutine homogenization_RGC_volumePenalty(&
enddo enddo
return
endsubroutine endsubroutine
!******************************************************************** !********************************************************************
@ -1159,8 +1141,6 @@ function homogenization_RGC_surfaceCorrection(&
sqrt(homogenization_RGC_surfaceCorrection(iBase))*detF sqrt(homogenization_RGC_surfaceCorrection(iBase))*detF
enddo enddo
return
endfunction endfunction
!******************************************************************** !********************************************************************
@ -1196,8 +1176,6 @@ function homogenization_RGC_equivalentModuli(&
!* Obtain the length of Burgers vector !* Obtain the length of Burgers vector
homogenization_RGC_equivalentModuli(2) = constitutive_averageBurgers(grainID,ip,el) homogenization_RGC_equivalentModuli(2) = constitutive_averageBurgers(grainID,ip,el)
return
endfunction endfunction
!******************************************************************** !********************************************************************
@ -1228,8 +1206,6 @@ function homogenization_RGC_relaxationVector(&
if (iNum .gt. 0_pInt) homogenization_RGC_relaxationVector = state%p((3*iNum-2):(3*iNum)) if (iNum .gt. 0_pInt) homogenization_RGC_relaxationVector = state%p((3*iNum-2):(3*iNum))
! get the corresponding entries ! get the corresponding entries
return
endfunction endfunction
!******************************************************************** !********************************************************************
@ -1268,8 +1244,6 @@ function homogenization_RGC_interfaceNormal(&
! call flush(6) ! call flush(6)
! endif ! endif
return
endfunction endfunction
!******************************************************************** !********************************************************************
@ -1297,8 +1271,6 @@ function homogenization_RGC_getInterface(&
if (iDir < 0_pInt) & ! to have a correlation with coordinate/position in real space if (iDir < 0_pInt) & ! to have a correlation with coordinate/position in real space
homogenization_RGC_getInterface(1_pInt-iDir) = homogenization_RGC_getInterface(1_pInt-iDir)-1_pInt homogenization_RGC_getInterface(1_pInt-iDir) = homogenization_RGC_getInterface(1_pInt-iDir)-1_pInt
return
endfunction endfunction
!******************************************************************** !********************************************************************
@ -1324,8 +1296,6 @@ function homogenization_RGC_grain1to3(&
homogenization_RGC_grain1to3(2) = 1+mod((grain1-1)/nGDim(1),nGDim(2)) homogenization_RGC_grain1to3(2) = 1+mod((grain1-1)/nGDim(1),nGDim(2))
homogenization_RGC_grain1to3(1) = 1+mod((grain1-1),nGDim(1)) homogenization_RGC_grain1to3(1) = 1+mod((grain1-1),nGDim(1))
return
endfunction endfunction
!******************************************************************** !********************************************************************
@ -1350,8 +1320,6 @@ function homogenization_RGC_grain3to1(&
nGDim = homogenization_RGC_Ngrains(:,homID) nGDim = homogenization_RGC_Ngrains(:,homID)
homogenization_RGC_grain3to1 = grain3(1) + nGDim(1)*(grain3(2)-1) + nGDim(1)*nGDim(2)*(grain3(3)-1) homogenization_RGC_grain3to1 = grain3(1) + nGDim(1)*(grain3(2)-1) + nGDim(1)*nGDim(2)*(grain3(3)-1)
return
endfunction endfunction
!******************************************************************** !********************************************************************
@ -1393,8 +1361,6 @@ function homogenization_RGC_interface4to1(&
if ((iFace4D(4) == 0_pInt) .or. (iFace4D(4) == nGDim(3))) homogenization_RGC_interface4to1 = 0_pInt if ((iFace4D(4) == 0_pInt) .or. (iFace4D(4) == nGDim(3))) homogenization_RGC_interface4to1 = 0_pInt
endif endif
return
endfunction endfunction
!******************************************************************** !********************************************************************
@ -1441,8 +1407,6 @@ function homogenization_RGC_interface1to4(&
homogenization_RGC_interface1to4(4) = int(dble(iFace1D-nIntFace(2)-nIntFace(1)-1)/dble(nGDim(1))/dble(nGDim(2)))+1 homogenization_RGC_interface1to4(4) = int(dble(iFace1D-nIntFace(2)-nIntFace(1)-1)/dble(nGDim(1))/dble(nGDim(2)))+1
endif endif
return
endfunction endfunction
!******************************************************************** !********************************************************************
@ -1493,8 +1457,6 @@ subroutine homogenization_RGC_grainDeformation(&
F(:,:,iGrain) = F(:,:,iGrain) + avgF(:,:) ! relaxed deformation gradient F(:,:,iGrain) = F(:,:,iGrain) + avgF(:,:) ! relaxed deformation gradient
enddo enddo
return
endsubroutine endsubroutine
END MODULE END MODULE

59
code/lattice.f90
View File

@ -811,7 +811,6 @@ function lattice_initializeStructure(struct,CoverA)
logical :: processMe logical :: processMe
integer(pInt) lattice_initializeStructure integer(pInt) lattice_initializeStructure
processMe = .false. processMe = .false.
select case(struct(1:3)) ! check first three chars of structure name select case(struct(1:3)) ! check first three chars of structure name
@ -825,14 +824,14 @@ function lattice_initializeStructure(struct,CoverA)
if (lattice_fcc_Nstructure == 1_pInt) then ! me is first fcc structure if (lattice_fcc_Nstructure == 1_pInt) then ! me is first fcc structure
processMe = .true. processMe = .true.
do i = 1,myNslip ! calculate slip system vectors do i = 1,myNslip ! calculate slip system vectors
sd(:,i) = lattice_fcc_systemSlip(1:3,i)/sqrt(math_mul3x3(lattice_fcc_systemSlip(1:3,i),lattice_fcc_systemSlip(1:3,i))) 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(:,i) = lattice_fcc_systemSlip(4:6,i)/sqrt(math_mul3x3(lattice_fcc_systemSlip(4:6,i),lattice_fcc_systemSlip(4:6,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(:,i) = math_vectorproduct(sd(:,i),sn(:,i)) st(1:3,i) = math_vectorproduct(sd(1:3,i),sn(1:3,i))
enddo enddo
do i = 1,myNtwin ! calculate twin system vectors and (assign) shears do i = 1,myNtwin ! calculate twin system vectors and (assign) shears
td(:,i) = lattice_fcc_systemTwin(1:3,i)/sqrt(math_mul3x3(lattice_fcc_systemTwin(1:3,i),lattice_fcc_systemTwin(1:3,i))) 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(:,i) = lattice_fcc_systemTwin(4:6,i)/sqrt(math_mul3x3(lattice_fcc_systemTwin(4:6,i),lattice_fcc_systemTwin(4:6,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(:,i) = math_vectorproduct(td(:,i),tn(:,i)) tt(1:3,i) = math_vectorproduct(td(1:3,i),tn(1:3,i))
ts(i) = lattice_fcc_shearTwin(i) ts(i) = lattice_fcc_shearTwin(i)
enddo enddo
interactionSlipSlip => lattice_fcc_interactionSlipSlip interactionSlipSlip => lattice_fcc_interactionSlipSlip
@ -851,14 +850,14 @@ function lattice_initializeStructure(struct,CoverA)
if (lattice_bcc_Nstructure == 1_pInt) then ! me is first bcc structure if (lattice_bcc_Nstructure == 1_pInt) then ! me is first bcc structure
processMe = .true. processMe = .true.
do i = 1,myNslip ! calculate slip system vectors do i = 1,myNslip ! calculate slip system vectors
sd(:,i) = lattice_bcc_systemSlip(1:3,i)/sqrt(math_mul3x3(lattice_bcc_systemSlip(1:3,i),lattice_bcc_systemSlip(1:3,i))) 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(:,i) = lattice_bcc_systemSlip(4:6,i)/sqrt(math_mul3x3(lattice_bcc_systemSlip(4:6,i),lattice_bcc_systemSlip(4:6,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(:,i) = math_vectorproduct(sd(:,i),sn(:,i)) st(1:3,i) = math_vectorproduct(sd(1:3,i),sn(1:3,i))
enddo enddo
do i = 1,myNtwin ! calculate twin system vectors and (assign) shears do i = 1,myNtwin ! calculate twin system vectors and (assign) shears
td(:,i) = lattice_bcc_systemTwin(1:3,i)/sqrt(math_mul3x3(lattice_bcc_systemTwin(1:3,i),lattice_bcc_systemTwin(1:3,i))) 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(:,i) = lattice_bcc_systemTwin(4:6,i)/sqrt(math_mul3x3(lattice_bcc_systemTwin(4:6,i),lattice_bcc_systemTwin(4:6,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(:,i) = math_vectorproduct(td(:,i),tn(:,i)) tt(1:3,i) = math_vectorproduct(td(1:3,i),tn(1:3,i))
ts(i) = lattice_bcc_shearTwin(i) ts(i) = lattice_bcc_shearTwin(i)
enddo enddo
interactionSlipSlip => lattice_bcc_interactionSlipSlip interactionSlipSlip => lattice_bcc_interactionSlipSlip
@ -885,9 +884,9 @@ function lattice_initializeStructure(struct,CoverA)
hex_n(2) = (lattice_hex_systemSlip(5,i)+2.0_pReal*lattice_hex_systemSlip(6,i))/sqrt(3.0_pReal) hex_n(2) = (lattice_hex_systemSlip(5,i)+2.0_pReal*lattice_hex_systemSlip(6,i))/sqrt(3.0_pReal)
hex_n(3) = lattice_hex_systemSlip(8,i)/CoverA hex_n(3) = lattice_hex_systemSlip(8,i)/CoverA
sd(:,i) = hex_d/sqrt(math_mul3x3(hex_d,hex_d)) sd(1:3,i) = hex_d/sqrt(math_mul3x3(hex_d,hex_d))
sn(:,i) = hex_n/sqrt(math_mul3x3(hex_n,hex_n)) sn(1:3,i) = hex_n/sqrt(math_mul3x3(hex_n,hex_n))
st(:,i) = math_vectorproduct(sd(:,i),sn(:,i)) st(1:3,i) = math_vectorproduct(sd(1:3,i),sn(1:3,i))
enddo enddo
do i = 1,myNtwin do i = 1,myNtwin
hex_d(1) = lattice_hex_systemTwin(1,i)*1.5_pReal hex_d(1) = lattice_hex_systemTwin(1,i)*1.5_pReal
@ -897,9 +896,9 @@ function lattice_initializeStructure(struct,CoverA)
hex_n(2) = (lattice_hex_systemTwin(5,i)+2.0_pReal*lattice_hex_systemTwin(6,i))/sqrt(3.0_pReal) hex_n(2) = (lattice_hex_systemTwin(5,i)+2.0_pReal*lattice_hex_systemTwin(6,i))/sqrt(3.0_pReal)
hex_n(3) = lattice_hex_systemTwin(8,i)/CoverA hex_n(3) = lattice_hex_systemTwin(8,i)/CoverA
td(:,i) = hex_d/sqrt(math_mul3x3(hex_d,hex_d)) td(1:3,i) = hex_d/sqrt(math_mul3x3(hex_d,hex_d))
tn(:,i) = hex_n/sqrt(math_mul3x3(hex_n,hex_n)) tn(1:3,i) = hex_n/sqrt(math_mul3x3(hex_n,hex_n))
tt(:,i) = math_vectorproduct(td(:,i),tn(:,i)) tt(1:3,i) = math_vectorproduct(td(1:3,i),tn(1:3,i))
select case(lattice_hex_shearTwin(i)) ! from Christian & Mahajan 1995 p.29 select case(lattice_hex_shearTwin(i)) ! from Christian & Mahajan 1995 p.29
case (1) ! {10.2}<-10.1> case (1) ! {10.2}<-10.1>
@ -924,19 +923,19 @@ function lattice_initializeStructure(struct,CoverA)
if (myStructure > lattice_Nstructure) & if (myStructure > lattice_Nstructure) &
call IO_error(666,0,0,0,'structure index too large') ! check for memory leakage call IO_error(666,0,0,0,'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,myNslip ! store slip system vectors and Schmid matrix for my structure
lattice_sd(:,i,myStructure) = sd(:,i) lattice_sd(1:3,i,myStructure) = sd(1:3,i)
lattice_st(:,i,myStructure) = st(:,i) lattice_st(1:3,i,myStructure) = st(1:3,i)
lattice_sn(:,i,myStructure) = sn(:,i) lattice_sn(1:3,i,myStructure) = sn(1:3,i)
lattice_Sslip(:,:,i,myStructure) = math_tensorproduct(sd(:,i),sn(:,i)) lattice_Sslip(1:3,1:3,i,myStructure) = math_tensorproduct(sd(1:3,i),sn(1:3,i))
lattice_Sslip_v(:,i,myStructure) = math_Mandel33to6(math_symmetric3x3(lattice_Sslip(:,:,i,myStructure))) lattice_Sslip_v(1:6,i,myStructure) = math_Mandel33to6(math_symmetric3x3(lattice_Sslip(1:3,1:3,i,myStructure)))
enddo enddo
do i = 1,myNtwin ! store twin system vectors and Schmid plus rotation matrix for my structure do i = 1,myNtwin ! store twin system vectors and Schmid plus rotation matrix for my structure
lattice_td(:,i,myStructure) = td(:,i) lattice_td(1:3,i,myStructure) = td(1:3,i)
lattice_tt(:,i,myStructure) = tt(:,i) lattice_tt(1:3,i,myStructure) = tt(1:3,i)
lattice_tn(:,i,myStructure) = tn(:,i) lattice_tn(1:3,i,myStructure) = tn(1:3,i)
lattice_Stwin(:,:,i,myStructure) = math_tensorproduct(td(:,i),tn(:,i)) lattice_Stwin(1:3,1:3,i,myStructure) = math_tensorproduct(td(1:3,i),tn(1:3,i))
lattice_Stwin_v(:,i,myStructure) = math_Mandel33to6(math_symmetric3x3(lattice_Stwin(:,:,i,myStructure))) lattice_Stwin_v(1:6,i,myStructure) = math_Mandel33to6(math_symmetric3x3(lattice_Stwin(1:3,1:3,i,myStructure)))
lattice_Qtwin(:,:,i,myStructure) = math_AxisAngleToR(tn(:,i),180.0_pReal*inRad) lattice_Qtwin(1:3,1:3,i,myStructure) = math_AxisAngleToR(tn(1:3,i),180.0_pReal*inRad)
lattice_shearTwin(i,myStructure) = ts(i) lattice_shearTwin(i,myStructure) = ts(i)
enddo enddo
lattice_NslipSystem(1:lattice_maxNslipFamily,myStructure) = myNslipSystem ! number of slip systems in each family lattice_NslipSystem(1:lattice_maxNslipFamily,myStructure) = myNslipSystem ! number of slip systems in each family

96
code/makefile
View File

@ -10,7 +10,7 @@
# Install fftw3 (v3.2.2 is tested) with "./configure --enable-threads --enable-float" and "make", "make install" is not needed # Install fftw3 (v3.2.2 is tested) with "./configure --enable-threads --enable-float" and "make", "make install" is not needed
# as long as the two library files are copied to the source code directory. # as long as the two library files are copied to the source code directory.
# OPTIONS = standart (alternative): meaning # OPTIONS = standard (alternative): meaning
#------------------------------------------------------------- #-------------------------------------------------------------
# PRECISION = double (single): floating point precision # PRECISION = double (single): floating point precision
# F90 = ifort (gfortran): compiler, choose Intel or GNU # F90 = ifort (gfortran): compiler, choose Intel or GNU
@ -19,6 +19,25 @@
# OPTIMIZATION = DEFENSIVE (OFF,AGGRESSIVE): Optimization mode, O0, O2, O3 # OPTIMIZATION = DEFENSIVE (OFF,AGGRESSIVE): Optimization mode, O0, O2, O3
# OPENMP = TRUE (FALSE): OpenMP multiprocessor support # OPENMP = TRUE (FALSE): OpenMP multiprocessor support
# ACML = OFF (ON): link with AMD math core library (v. 4.4 need to be installed) # ACML = OFF (ON): link with AMD math core library (v. 4.4 need to be installed)
# PREFIX: specifie an arbitrary prefix
# COMPILERNAME = overwrite name of Compiler, e.g. using mpich-g90 instead of ifort
# Here are some usefull debugging switches. Switch on by uncommenting last line:
#--------------------------------------------------------------------------------
# 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. Porblably it helps to also unlimit other limits
DEBUG4 =-heap-arrays
#SUFFIX =$(DEBUG1) $(DEBUG2) $(DEBUG3)
#SUFFIX =$(DEBUG1) $(DEBUG3)
ifeq ($(F90), ) ifeq ($(F90), )
F90 =ifort F90 =ifort
@ -43,7 +62,7 @@ endif
endif endif
ifneq ($(OPENMP), OFF) ifneq ($(OPENMP), OFF)
OPENMP_FLAG_ifort =-openmp OPENMP_FLAG_ifort =-openmp -openmp-report0 -parallel
OPENMP_FLAG_gfortran =-fopenmp OPENMP_FLAG_gfortran =-fopenmp
OPENMP =ON OPENMP =ON
endif endif
@ -59,49 +78,46 @@ OPTIMIZATION_OFF_ifort =-O0
OPTIMIZATION_OFF_gfortran =-O0 OPTIMIZATION_OFF_gfortran =-O0
OPTIMIZATION_DEFENSIVE_ifort =-O2 OPTIMIZATION_DEFENSIVE_ifort =-O2
OPTIMIZATION_DEFENSIVE_gfortran =-O2 OPTIMIZATION_DEFENSIVE_gfortran =-O2
OPTIMIZATION_AGGRESSIVE_ifort =-O3 -static $(PORTABLE_SWITCH) OPTIMIZATION_AGGRESSIVE_ifort =-O3 $(PORTABLE_SWITCH) -ip
OPTIMIZATION_AGGRESSIVE_gfortran =-O3 OPTIMIZATION_AGGRESSIVE_gfortran =-O3
COMPILE_OPTIONS_ifort =-fpp -diag-disable 8291,8290 COMPILE_OPTIONS_ifort =-fpp -diag-disable 8291,8290
COMPILE_OPTIONS_gfortran =-xf95-cpp-input COMPILE_OPTIONS_gfortran =-xf95-cpp-input -ffree-line-length-none
HEAP_ARRAYS_ifort =-heap-arrays 500000000
HEAP_ARRAYS_gfortran =
COMPILE =${OPENMP_FLAG_${F90}} ${COMPILE_OPTIONS_${F90}} ${OPTIMIZATION_${OPTIMIZATION}_${F90}} -c COMPILE =${OPENMP_FLAG_${F90}} ${COMPILE_OPTIONS_${F90}} ${OPTIMIZATION_${OPTIMIZATION}_${F90}} -c
COMPILE_HEAP =$(COMPILE) ${HEAP_ARRAYS_${F90}} COMPILE_MAXOPTI =${OPENMP_FLAG_${F90}} ${COMPILE_OPTIONS_${F90}} ${OPTIMIZATION_${MAXOPTI}_${F90}} -c
COMPILE_HEAP_MAXOPTI =${OPENMP_FLAG_${F90}} ${COMPILE_OPTIONS_${F90}} ${OPTIMIZATION_${MAXOPTI}_${F90}} ${HEAP_ARRAYS_${F90}} -c
ifndef COMPILERNAME
COMPILERNAME=$(F90)
endif
ifeq ($(PRECISION),single) ifeq ($(PRECISION),single)
DAMASK_spectral_single.exe: DAMASK_spectral_single.o CPFEM.a DAMASK_spectral_single.exe: DAMASK_spectral_single.o CPFEM.a
$(F90) ${OPENMP_FLAG_${F90}} -o DAMASK_spectral_single.exe DAMASK_spectral_single.o CPFEM.a include/libfftw3f_threads.a include/libfftw3f.a constitutive.a advanced.a basics.a\ $(PREFIX) $(COMPILERNAME) ${OPENMP_FLAG_${F90}} -o DAMASK_spectral_single.exe DAMASK_spectral_single.o CPFEM.a include/libfftw3f_threads.a include/libfftw3f.a\
-lpthread ${BLAS_${OPENMP}_${F90}} constitutive.a advanced.a basics.a -lpthread ${BLAS_${OPENMP}_${F90}}
DAMASK_spectral_single.o: DAMASK_spectral_single.f90 CPFEM.o DAMASK_spectral_single.o: DAMASK_spectral_single.f90 CPFEM.o
$(F90) $(COMPILE_HEAP_MAXOPTI) DAMASK_spectral_single.f90 $(PREFIX) $(COMPILERNAME) $(COMPILE_MAXOPTI) DAMASK_spectral_single.f90 $(SUFFIX)
else else
DAMASK_spectral.exe: DAMASK_spectral.o CPFEM.a DAMASK_spectral.exe: DAMASK_spectral.o CPFEM.a
$(F90) ${OPENMP_FLAG_${F90}} -o DAMASK_spectral.exe DAMASK_spectral.o CPFEM.a include/libfftw3_threads.a include/libfftw3.a constitutive.a advanced.a basics.a\ $(PREFIX) $(COMPILERNAME) ${OPENMP_FLAG_${F90}} -o DAMASK_spectral.exe DAMASK_spectral.o CPFEM.a include/libfftw3_threads.a include/libfftw3.a\
-lpthread ${BLAS_${OPENMP}_${F90}} constitutive.a advanced.a basics.a -lpthread ${BLAS_${OPENMP}_${F90}}
DAMASK_spectral.o: DAMASK_spectral.f90 CPFEM.o DAMASK_spectral.o: DAMASK_spectral.f90 CPFEM.o
$(F90) $(COMPILE_HEAP_MAXOPTI) DAMASK_spectral.f90 $(PREFIX) $(COMPILERNAME) $(COMPILE_MAXOPTI) DAMASK_spectral.f90 $(SUFFIX)
endif endif
CPFEM.a: CPFEM.o CPFEM.a: CPFEM.o
ar rc CPFEM.a homogenization.o homogenization_RGC.o homogenization_isostrain.o crystallite.o CPFEM.o constitutive.o ar rc CPFEM.a homogenization.o homogenization_RGC.o homogenization_isostrain.o crystallite.o CPFEM.o constitutive.o
CPFEM.o: CPFEM.f90 homogenization.o CPFEM.o: CPFEM.f90 homogenization.o
$(F90) $(COMPILE_HEAP) CPFEM.f90 $(PREFIX) $(COMPILERNAME) $(COMPILE) CPFEM.f90 $(SUFFIX)
homogenization.o: homogenization.f90 homogenization_isostrain.o homogenization_RGC.o crystallite.o homogenization.o: homogenization.f90 homogenization_isostrain.o homogenization_RGC.o crystallite.o
$(F90) $(COMPILE_HEAP) homogenization.f90 $(PREFIX) $(COMPILERNAME) $(COMPILE) homogenization.f90 $(SUFFIX)
homogenization_RGC.o: homogenization_RGC.f90 constitutive.a homogenization_RGC.o: homogenization_RGC.f90 constitutive.a
$(F90) $(COMPILE_HEAP) homogenization_RGC.f90 $(PREFIX) $(COMPILERNAME) $(COMPILE) homogenization_RGC.f90 $(SUFFIX)
homogenization_isostrain.o: homogenization_isostrain.f90 basics.a advanced.a homogenization_isostrain.o: homogenization_isostrain.f90 basics.a advanced.a
$(F90) $(COMPILE_HEAP) homogenization_isostrain.f90 $(PREFIX) $(COMPILERNAME) $(COMPILE) homogenization_isostrain.f90 $(SUFFIX)
crystallite.o: crystallite.f90 constitutive.a crystallite.o: crystallite.f90 constitutive.a
$(F90) $(COMPILE_HEAP) crystallite.f90 $(PREFIX) $(COMPILERNAME) $(COMPILE) crystallite.f90 $(SUFFIX)
@ -109,22 +125,22 @@ constitutive.a: constitutive.o
ar rc constitutive.a constitutive.o constitutive_titanmod.o constitutive_nonlocal.o constitutive_dislotwin.o constitutive_j2.o constitutive_phenopowerlaw.o basics.a advanced.a ar rc constitutive.a constitutive.o constitutive_titanmod.o constitutive_nonlocal.o constitutive_dislotwin.o constitutive_j2.o constitutive_phenopowerlaw.o basics.a advanced.a
constitutive.o: constitutive.f90 constitutive_titanmod.o constitutive_nonlocal.o constitutive_dislotwin.o constitutive_j2.o constitutive_phenopowerlaw.o constitutive.o: constitutive.f90 constitutive_titanmod.o constitutive_nonlocal.o constitutive_dislotwin.o constitutive_j2.o constitutive_phenopowerlaw.o
$(F90) $(COMPILE_HEAP) constitutive.f90 $(PREFIX) $(COMPILERNAME) $(COMPILE) constitutive.f90 $(SUFFIX)
constitutive_titanmod.o: constitutive_titanmod.f90 basics.a advanced.a constitutive_titanmod.o: constitutive_titanmod.f90 basics.a advanced.a
$(F90) $(COMPILE_HEAP) constitutive_titanmod.f90 $(PREFIX) $(COMPILERNAME) $(COMPILE) constitutive_titanmod.f90 $(SUFFIX)
constitutive_nonlocal.o: constitutive_nonlocal.f90 basics.a advanced.a constitutive_nonlocal.o: constitutive_nonlocal.f90 basics.a advanced.a
$(F90) $(COMPILE_HEAP) constitutive_nonlocal.f90 $(PREFIX) $(COMPILERNAME) $(COMPILE) constitutive_nonlocal.f90 $(SUFFIX)
constitutive_dislotwin.o: constitutive_dislotwin.f90 basics.a advanced.a constitutive_dislotwin.o: constitutive_dislotwin.f90 basics.a advanced.a
$(F90) $(COMPILE_HEAP) constitutive_dislotwin.f90 $(PREFIX) $(COMPILERNAME) $(COMPILE) constitutive_dislotwin.f90 $(SUFFIX)
constitutive_j2.o: constitutive_j2.f90 basics.a advanced.a constitutive_j2.o: constitutive_j2.f90 basics.a advanced.a
$(F90) $(COMPILE_HEAP) constitutive_j2.f90 $(PREFIX) $(COMPILERNAME) $(COMPILE) constitutive_j2.f90 $(SUFFIX)
constitutive_phenopowerlaw.o: constitutive_phenopowerlaw.f90 basics.a advanced.a constitutive_phenopowerlaw.o: constitutive_phenopowerlaw.f90 basics.a advanced.a
$(F90) $(COMPILE_HEAP) constitutive_phenopowerlaw.f90 $(PREFIX) $(COMPILERNAME) $(COMPILE) constitutive_phenopowerlaw.f90 $(SUFFIX)
@ -133,13 +149,13 @@ advanced.a: lattice.o
lattice.o: lattice.f90 material.o lattice.o: lattice.f90 material.o
$(F90) $(COMPILE_HEAP) lattice.f90 $(PREFIX) $(COMPILERNAME) $(COMPILE) lattice.f90 $(SUFFIX)
material.o: material.f90 mesh.o material.o: material.f90 mesh.o
$(F90) $(COMPILE_HEAP) material.f90 $(PREFIX) $(COMPILERNAME) $(COMPILE) material.f90 $(SUFFIX)
mesh.o: mesh.f90 FEsolving.o mesh.o: mesh.f90 FEsolving.o
$(F90) $(COMPILE_HEAP) mesh.f90 $(PREFIX) $(COMPILERNAME) $(COMPILE) mesh.f90 $(SUFFIX)
FEsolving.o: FEsolving.f90 basics.a FEsolving.o: FEsolving.f90 basics.a
$(F90) $(COMPILE_HEAP) FEsolving.f90 $(PREFIX) $(COMPILERNAME) $(COMPILE) FEsolving.f90 $(SUFFIX)
ifeq ($(PRECISION),single) ifeq ($(PRECISION),single)
basics.a: debug.o math.o basics.a: debug.o math.o
@ -150,25 +166,25 @@ basics.a: debug.o math.o
endif endif
debug.o: debug.f90 numerics.o debug.o: debug.f90 numerics.o
$(F90) $(COMPILE) debug.f90 $(PREFIX) $(COMPILERNAME) $(COMPILE) debug.f90 $(SUFFIX)
math.o: math.f90 numerics.o math.o: math.f90 numerics.o
$(F90) $(COMPILE) math.f90 $(PREFIX) $(COMPILERNAME) $(COMPILE) math.f90 $(SUFFIX)
numerics.o: numerics.f90 IO.o numerics.o: numerics.f90 IO.o
$(F90) $(COMPILE) numerics.f90 $(PREFIX) $(COMPILERNAME) $(COMPILE) numerics.f90 $(SUFFIX)
IO.o: IO.f90 DAMASK_spectral_interface.o IO.o: IO.f90 DAMASK_spectral_interface.o
$(F90) $(COMPILE) IO.f90 $(PREFIX) $(COMPILERNAME) $(COMPILE) IO.f90 $(SUFFIX)
ifeq ($(PRECISION),single) ifeq ($(PRECISION),single)
DAMASK_spectral_interface.o: DAMASK_spectral_interface.f90 prec_single.o DAMASK_spectral_interface.o: DAMASK_spectral_interface.f90 prec_single.o
$(F90) $(COMPILE) DAMASK_spectral_interface.f90 $(PREFIX) $(COMPILERNAME) $(COMPILE) DAMASK_spectral_interface.f90 $(SUFFIX)
prec_single.o: prec_single.f90 prec_single.o: prec_single.f90
$(F90) $(COMPILE) prec_single.f90 $(PREFIX) $(COMPILERNAME) $(COMPILE) prec_single.f90 $(SUFFIX)
else else
DAMASK_spectral_interface.o: DAMASK_spectral_interface.f90 prec.o DAMASK_spectral_interface.o: DAMASK_spectral_interface.f90 prec.o
$(F90) $(COMPILE) DAMASK_spectral_interface.f90 $(PREFIX) $(COMPILERNAME) $(COMPILE) DAMASK_spectral_interface.f90 $(SUFFIX)
prec.o: prec.f90 prec.o: prec.f90
$(F90) $(COMPILE) prec.f90 $(PREFIX) $(COMPILERNAME) $(COMPILE) prec.f90 $(SUFFIX)
endif endif

3
code/material.f90
View File

@ -245,7 +245,6 @@ subroutine material_parseHomogenization(file,myPart)
enddo enddo
100 homogenization_maxNgrains = maxval(homogenization_Ngrains,homogenization_active) 100 homogenization_maxNgrains = maxval(homogenization_Ngrains,homogenization_active)
return
endsubroutine endsubroutine
@ -800,8 +799,6 @@ subroutine material_populateGrains()
deallocate(textureOfGrain) deallocate(textureOfGrain)
deallocate(orientationOfGrain) deallocate(orientationOfGrain)
return
endsubroutine endsubroutine

87
code/mesh.f90
View File

@ -449,7 +449,6 @@
exit exit
endif endif
enddo enddo
return
endfunction endfunction
@ -1353,8 +1352,6 @@ FE_ipNeighbor(:FE_NipNeighbors(8),:FE_Nips(8),8) = & ! element 117
6, 7, 7, 6 & 6, 7, 7, 6 &
/),(/FE_NipFaceNodes,FE_NipNeighbors(10),FE_Nips(10)/)) /),(/FE_NipFaceNodes,FE_NipNeighbors(10),FE_Nips(10)/))
return
endsubroutine endsubroutine
@ -1392,9 +1389,7 @@ FE_ipNeighbor(:FE_NipNeighbors(8),:FE_Nips(8),8) = & ! element 117
endif endif
enddo enddo
620 return 620 endsubroutine
endsubroutine
!******************************************************************** !********************************************************************
@ -1441,9 +1436,7 @@ do
endif endif
enddo enddo
620 return 620 endsubroutine
endsubroutine
!******************************************************************** !********************************************************************
@ -1490,9 +1483,7 @@ endsubroutine
endif endif
enddo enddo
100 return 100 endsubroutine
endsubroutine
!******************************************************************** !********************************************************************
! count overall number of nodes and elements in mesh ! count overall number of nodes and elements in mesh
@ -1528,9 +1519,7 @@ endsubroutine
endif endif
enddo enddo
620 return 620 endsubroutine
endsubroutine
!******************************************************************** !********************************************************************
! count overall number of nodes and elements in mesh ! count overall number of nodes and elements in mesh
@ -1589,8 +1578,6 @@ endsubroutine
620 if (mesh_Nnodes < 2) call IO_error(900) 620 if (mesh_Nnodes < 2) call IO_error(900)
if (mesh_Nelems == 0) call IO_error(901) if (mesh_Nelems == 0) call IO_error(901)
return
endsubroutine endsubroutine
@ -1629,9 +1616,7 @@ endsubroutine
endif endif
enddo enddo
620 return 620 endsubroutine
endsubroutine
!******************************************************************** !********************************************************************
@ -1673,7 +1658,6 @@ endsubroutine
620 continue 620 continue
if (mesh_NelemSets == 0) call IO_error(902) if (mesh_NelemSets == 0) call IO_error(902)
return
endsubroutine endsubroutine
@ -1716,7 +1700,6 @@ endsubroutine
620 if (mesh_Nmaterials == 0) call IO_error(903) 620 if (mesh_Nmaterials == 0) call IO_error(903)
return
endsubroutine endsubroutine
@ -1731,7 +1714,6 @@ endsubroutine
implicit none implicit none
mesh_NcpElems = mesh_Nelems mesh_NcpElems = mesh_Nelems
return
endsubroutine endsubroutine
@ -1771,9 +1753,7 @@ endsubroutine
endif endif
enddo enddo
620 return 620 endsubroutine
endsubroutine
!******************************************************************** !********************************************************************
@ -1825,7 +1805,6 @@ endsubroutine
620 if (mesh_NcpElems == 0) call IO_error(906) 620 if (mesh_NcpElems == 0) call IO_error(906)
return
endsubroutine endsubroutine
@ -1865,9 +1844,7 @@ endsubroutine
endif endif
enddo enddo
640 return 640 endsubroutine
endsubroutine
!******************************************************************** !********************************************************************
@ -1917,7 +1894,6 @@ endsubroutine
if (mesh_mapElemSet(1,i) == 0) call IO_error(ID=904,ext_msg=mesh_nameElemSet(i)) if (mesh_mapElemSet(1,i) == 0) call IO_error(ID=904,ext_msg=mesh_nameElemSet(i))
enddo enddo
return
endsubroutine endsubroutine
@ -1984,7 +1960,6 @@ endsubroutine
if (mesh_nameMaterial(i)=='' .or. mesh_mapMaterial(i)=='') call IO_error(905) if (mesh_nameMaterial(i)=='' .or. mesh_mapMaterial(i)=='') call IO_error(905)
enddo enddo
return
endsubroutine endsubroutine
@ -2004,9 +1979,7 @@ endsubroutine
allocate (mesh_mapFEtoCPnode(2,mesh_Nnodes)) ; mesh_mapFEtoCPnode = 0_pInt allocate (mesh_mapFEtoCPnode(2,mesh_Nnodes)) ; mesh_mapFEtoCPnode = 0_pInt
forall (i = 1:mesh_Nnodes) & forall (i = 1:mesh_Nnodes) &
mesh_mapFEtoCPnode(:,i) = i mesh_mapFEtoCPnode(1:2,i) = i
return
endsubroutine endsubroutine
@ -2055,8 +2028,6 @@ endsubroutine
650 call qsort(mesh_mapFEtoCPnode,1,size(mesh_mapFEtoCPnode,2)) 650 call qsort(mesh_mapFEtoCPnode,1,size(mesh_mapFEtoCPnode,2))
return
endsubroutine endsubroutine
@ -2119,7 +2090,6 @@ endsubroutine
650 call qsort(mesh_mapFEtoCPnode,1,size(mesh_mapFEtoCPnode,2)) 650 call qsort(mesh_mapFEtoCPnode,1,size(mesh_mapFEtoCPnode,2))
if (size(mesh_mapFEtoCPnode) == 0) call IO_error(908) if (size(mesh_mapFEtoCPnode) == 0) call IO_error(908)
return
endsubroutine endsubroutine
@ -2139,9 +2109,7 @@ endsubroutine
allocate (mesh_mapFEtoCPelem(2,mesh_NcpElems)) ; mesh_mapFEtoCPelem = 0_pInt allocate (mesh_mapFEtoCPelem(2,mesh_NcpElems)) ; mesh_mapFEtoCPelem = 0_pInt
forall (i = 1:mesh_NcpElems) & forall (i = 1:mesh_NcpElems) &
mesh_mapFEtoCPelem(:,i) = i mesh_mapFEtoCPelem(1:2,i) = i
return
endsubroutine endsubroutine
@ -2191,7 +2159,6 @@ endsubroutine
660 call qsort(mesh_mapFEtoCPelem,1,size(mesh_mapFEtoCPelem,2)) ! should be mesh_NcpElems 660 call qsort(mesh_mapFEtoCPelem,1,size(mesh_mapFEtoCPelem,2)) ! should be mesh_NcpElems
return
endsubroutine endsubroutine
@ -2254,7 +2221,6 @@ endsubroutine
if (size(mesh_mapFEtoCPelem) < 2) call IO_error(907) if (size(mesh_mapFEtoCPelem) < 2) call IO_error(907)
return
endsubroutine endsubroutine
@ -2328,8 +2294,7 @@ subroutine mesh_marc_count_cpSizes (unit)
endif endif
enddo enddo
630 return 630 endsubroutine
endsubroutine
!******************************************************************** !********************************************************************
@ -2392,9 +2357,7 @@ subroutine mesh_marc_count_cpSizes (unit)
endif endif
enddo enddo
620 return 620 endsubroutine
endsubroutine
!******************************************************************** !********************************************************************
@ -2483,9 +2446,7 @@ subroutine mesh_marc_count_cpSizes (unit)
mesh_node = mesh_node0 mesh_node = mesh_node0
100 return 100 endsubroutine
endsubroutine
!******************************************************************** !********************************************************************
@ -2528,7 +2489,6 @@ subroutine mesh_marc_count_cpSizes (unit)
enddo enddo
670 mesh_node = mesh_node0 670 mesh_node = mesh_node0
return
endsubroutine endsubroutine
@ -2588,7 +2548,6 @@ return
670 if (size(mesh_node0,2) /= mesh_Nnodes) call IO_error(909) 670 if (size(mesh_node0,2) /= mesh_Nnodes) call IO_error(909)
mesh_node = mesh_node0 mesh_node = mesh_node0
return
endsubroutine endsubroutine
@ -2676,9 +2635,7 @@ return
mesh_maxValStateVar(2) = max(mesh_maxValStateVar(2),mesh_element(4,e)) mesh_maxValStateVar(2) = max(mesh_maxValStateVar(2),mesh_element(4,e))
enddo enddo
110 return 110 endsubroutine
endsubroutine
@ -2762,9 +2719,7 @@ return
endif endif
enddo enddo
630 return 630 endsubroutine
endsubroutine
@ -2867,9 +2822,7 @@ return
endselect endselect
enddo enddo
630 return 630 endsubroutine
endsubroutine
!******************************************************************** !********************************************************************
@ -3247,7 +3200,6 @@ endsubroutine
allocate(mesh_ipArea(mesh_maxNipNeighbors,mesh_maxNips,mesh_NcpElems)) ; mesh_ipArea = 0.0_pReal allocate(mesh_ipArea(mesh_maxNipNeighbors,mesh_maxNips,mesh_NcpElems)) ; mesh_ipArea = 0.0_pReal
allocate(mesh_ipAreaNormal(3,mesh_maxNipNeighbors,mesh_maxNips,mesh_NcpElems)) ; mesh_ipAreaNormal = 0.0_pReal allocate(mesh_ipAreaNormal(3,mesh_maxNipNeighbors,mesh_maxNips,mesh_NcpElems)) ; mesh_ipAreaNormal = 0.0_pReal
do e = 1,mesh_NcpElems ! loop over cpElems do e = 1,mesh_NcpElems ! loop over cpElems
t = mesh_element(2,e) ! get elemType t = mesh_element(2,e) ! get elemType
do i = 1,FE_Nips(t) ! loop over IPs of elem do i = 1,FE_Nips(t) ! loop over IPs of elem
@ -3259,14 +3211,13 @@ endsubroutine
area(j,n) = sqrt(sum(normal(:,j,n)*normal(:,j,n))) ! and area area(j,n) = sqrt(sum(normal(:,j,n)*normal(:,j,n))) ! and area
end forall end forall
forall (n = 1:FE_NipFaceNodes, j = 1:Ntriangles, area(j,n) > 0.0_pReal) & forall (n = 1:FE_NipFaceNodes, j = 1:Ntriangles, area(j,n) > 0.0_pReal) &
normal(:,j,n) = normal(:,j,n) / area(j,n) ! make unit normal normal(1:3,j,n) = normal(1:3,j,n) / area(j,n) ! make unit normal
mesh_ipArea(f,i,e) = sum(area) / (FE_NipFaceNodes*2.0_pReal) ! area of parallelograms instead of triangles 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) / count(area > 0.0_pReal) ! average of all valid normals mesh_ipAreaNormal(:,f,i,e) = sum(sum(normal,3),2) / count(area > 0.0_pReal) ! average of all valid normals
enddo enddo
enddo enddo
enddo enddo
return
endsubroutine endsubroutine
@ -3399,11 +3350,11 @@ if (debug_verbosity > 0) then
write (6,*) mesh_maxValStateVar(1), " : maximum homogenization index" write (6,*) mesh_maxValStateVar(1), " : maximum homogenization index"
write (6,*) mesh_maxValStateVar(2), " : maximum microstructure index" write (6,*) mesh_maxValStateVar(2), " : maximum microstructure index"
write (6,*) write (6,*)
write (fmt,"(a,i5,a)") "(9(x),a2,x,",mesh_maxValStateVar(2),"(i8))" write (fmt,"(a,i32.32,a)") "(9(x),a2,x,",mesh_maxValStateVar(2),"(i8))"
write (6,fmt) "+-",math_range(mesh_maxValStateVar(2)) write (6,fmt) "+-",math_range(mesh_maxValStateVar(2))
write (fmt,"(a,i5,a)") "(i8,x,a2,x,",mesh_maxValStateVar(2),"(i8))" write (fmt,"(a,i32.32,a)") "(i8,x,a2,x,",mesh_maxValStateVar(2),"(i8))"
do i=1,mesh_maxValStateVar(1) ! loop over all (possibly assigned) homogenizations do i=1,mesh_maxValStateVar(1) ! loop over all (possibly assigned) homogenizations
write (6,fmt) i,"| ",mesh_HomogMicro(i,:) ! loop over all (possibly assigned) microstrcutures write (6,fmt) i,"| ",mesh_HomogMicro(i,:) ! loop over all (possibly assigned) microstructures
enddo enddo
write(6,*) write(6,*)
write(6,*) "Input Parser: ADDITIONAL MPIE OPTIONS" write(6,*) "Input Parser: ADDITIONAL MPIE OPTIONS"