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

Merge branch 'no-pInt-implicitNone' into 'development' 

No pInt and implicit none

See merge request damask/DAMASK!81
This commit is contained in:
Martin Diehl 2019-05-15 21:40:13 +02:00
parent 880bc330f0 346c7c4a7f
commit 789420c9d6
27 changed files with 1709 additions and 1941 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
src/CPFEM.f90
View File

@ -4,9 +4,7 @@
!> @brief CPFEM engine
!--------------------------------------------------------------------------------------------------
module CPFEM
use prec, only: &
pReal, &
pInt
use prec
implicit none
private
@ -57,8 +55,6 @@ contains
!> @brief call (thread safe) all module initializations
!--------------------------------------------------------------------------------------------------
subroutine CPFEM_initAll(el,ip)
use prec, only: &
prec_init
use numerics, only: &
numerics_init
use debug, only: &
@ -91,7 +87,6 @@ subroutine CPFEM_initAll(el,ip)
IO_init
use DAMASK_interface
implicit none
integer(pInt), intent(in) :: el, & !< FE el number
ip !< FE integration point number
@ -155,7 +150,6 @@ subroutine CPFEM_init
crystallite_Li0, &
crystallite_S0
implicit none
integer :: k,l,m,ph,homog
write(6,'(/,a)') ' <<<+- CPFEM init -+>>>'
@ -325,7 +319,6 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt
IO_warning
use DAMASK_interface
implicit none
integer(pInt), intent(in) :: elFE, & !< FE element number
ip !< integration point number
real(pReal), intent(in) :: dt !< time increment
@ -639,8 +632,6 @@ end subroutine CPFEM_general
!> @brief triggers writing of the results
!--------------------------------------------------------------------------------------------------
subroutine CPFEM_results(inc,time)
use prec, only: &
pInt
#ifdef DAMASK_HDF5
use results
use HDF5_utilities
@ -650,7 +641,6 @@ subroutine CPFEM_results(inc,time)
use crystallite, only: &
crystallite_results
implicit none
integer(pInt), intent(in) :: inc
real(pReal), intent(in) :: time

39
src/CPFEM2.f90
View File

@ -12,6 +12,7 @@ module CPFEM2
CPFEM_age, &
CPFEM_initAll, &
CPFEM_results
contains
@ -20,7 +21,6 @@ contains
!--------------------------------------------------------------------------------------------------
subroutine CPFEM_initAll()
use prec, only: &
pInt, &
prec_init
use numerics, only: &
numerics_init
@ -57,8 +57,6 @@ subroutine CPFEM_initAll()
FEM_Zoo_init
#endif
implicit none
call DAMASK_interface_init ! Spectral and FEM interface to commandline
call prec_init
call IO_init
@ -87,8 +85,6 @@ end subroutine CPFEM_initAll
!> @brief allocate the arrays defined in module CPFEM and initialize them
!--------------------------------------------------------------------------------------------------
subroutine CPFEM_init
use prec, only: &
pInt, pReal
use IO, only: &
IO_error
use numerics, only: &
@ -124,8 +120,8 @@ subroutine CPFEM_init
use DAMASK_interface, only: &
getSolverJobName
implicit none
integer(pInt) :: ph,homog
integer :: ph,homog
character(len=1024) :: rankStr, PlasticItem, HomogItem
integer(HID_T) :: fileHandle, groupPlasticID, groupHomogID
@ -134,7 +130,7 @@ subroutine CPFEM_init
! *** restore the last converged values of each essential variable from the binary file
if (restartRead) then
if (iand(debug_level(debug_CPFEM), debug_levelExtensive) /= 0_pInt) then
if (iand(debug_level(debug_CPFEM), debug_levelExtensive) /= 0) then
write(6,'(a)') '<< CPFEM >> restored state variables of last converged step from hdf5 file'
flush(6)
endif
@ -152,14 +148,14 @@ subroutine CPFEM_init
call HDF5_read(fileHandle,crystallite_S0, 'convergedS')
groupPlasticID = HDF5_openGroup(fileHandle,'PlasticPhases')
do ph = 1_pInt,size(phase_plasticity)
do ph = 1,size(phase_plasticity)
write(PlasticItem,*) ph,'_'
call HDF5_read(groupPlasticID,plasticState(ph)%state0,trim(PlasticItem)//'convergedStateConst')
enddo
call HDF5_closeGroup(groupPlasticID)
groupHomogID = HDF5_openGroup(fileHandle,'HomogStates')
do homog = 1_pInt, material_Nhomogenization
do homog = 1, material_Nhomogenization
write(HomogItem,*) homog,'_'
call HDF5_read(groupHomogID,homogState(homog)%state0, trim(HomogItem)//'convergedStateHomog')
enddo
@ -178,8 +174,7 @@ end subroutine CPFEM_init
!--------------------------------------------------------------------------------------------------
subroutine CPFEM_age()
use prec, only: &
pReal, &
pInt
pReal
use numerics, only: &
worldrank
use debug, only: &
@ -224,12 +219,11 @@ subroutine CPFEM_age()
use DAMASK_interface, only: &
getSolverJobName
implicit none
integer(pInt) :: i, ph, homog, mySource
integer :: i, ph, homog, mySource
character(len=32) :: rankStr, PlasticItem, HomogItem
integer(HID_T) :: fileHandle, groupPlastic, groupHomog
if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0_pInt) &
if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0) &
write(6,'(a)') '<< CPFEM >> aging states'
crystallite_F0 = crystallite_partionedF
@ -246,14 +240,14 @@ subroutine CPFEM_age()
do mySource = 1,phase_Nsources(i)
sourceState(i)%p(mySource)%state0 = sourceState(i)%p(mySource)%state
enddo; enddo
do homog = 1_pInt, material_Nhomogenization
do homog = 1, material_Nhomogenization
homogState (homog)%state0 = homogState (homog)%state
thermalState (homog)%state0 = thermalState (homog)%state
damageState (homog)%state0 = damageState (homog)%state
enddo
if (restartWrite) then
if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0_pInt) &
if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0) &
write(6,'(a)') '<< CPFEM >> writing restart variables of last converged step to hdf5 file'
write(rankStr,'(a1,i0)')'_',worldrank
@ -268,14 +262,14 @@ subroutine CPFEM_age()
call HDF5_write(fileHandle,crystallite_S0, 'convergedS')
groupPlastic = HDF5_addGroup(fileHandle,'PlasticPhases')
do ph = 1_pInt,size(phase_plasticity)
do ph = 1,size(phase_plasticity)
write(PlasticItem,*) ph,'_'
call HDF5_write(groupPlastic,plasticState(ph)%state0,trim(PlasticItem)//'convergedStateConst')
enddo
call HDF5_closeGroup(groupPlastic)
groupHomog = HDF5_addGroup(fileHandle,'HomogStates')
do homog = 1_pInt, material_Nhomogenization
do homog = 1, material_Nhomogenization
write(HomogItem,*) homog,'_'
call HDF5_write(groupHomog,homogState(homog)%state0,trim(HomogItem)//'convergedStateHomog')
enddo
@ -285,7 +279,7 @@ subroutine CPFEM_age()
restartWrite = .false.
endif
if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0_pInt) &
if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0) &
write(6,'(a)') '<< CPFEM >> done aging states'
end subroutine CPFEM_age
@ -295,8 +289,6 @@ end subroutine CPFEM_age
!> @brief triggers writing of the results
!--------------------------------------------------------------------------------------------------
subroutine CPFEM_results(inc,time)
use prec, only: &
pInt
use results
use HDF5_utilities
use homogenization, only: &
@ -306,8 +298,7 @@ subroutine CPFEM_results(inc,time)
use crystallite, only: &
crystallite_results
implicit none
integer(pInt), intent(in) :: inc
integer, intent(in) :: inc
real(pReal), intent(in) :: time
call results_openJobFile

44
src/DAMASK_marc.f90
View File

@ -96,14 +96,12 @@ end subroutine DAMASK_interface_init
!> @brief solver job name (no extension) as combination of geometry and load case name
!--------------------------------------------------------------------------------------------------
function getSolverJobName()
use prec, only: &
pReal, &
pInt
use prec
implicit none
character(1024) :: getSolverJobName, inputName
character(len=*), parameter :: pathSep = achar(47)//achar(92) ! forward and backward slash
integer(pInt) :: extPos
integer :: extPos
getSolverJobName=''
inputName=''
@ -133,9 +131,7 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, &
dispt,coord,ffn,frotn,strechn,eigvn,ffn1,frotn1, &
strechn1,eigvn1,ncrd,itel,ndeg,ndm,nnode, &
jtype,lclass,ifr,ifu)
use prec, only: &
pReal, &
pInt
use prec
use numerics, only: &
!$ DAMASK_NumThreadsInt, &
numerics_unitlength, &
@ -180,7 +176,7 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, &
implicit none
!$ include "omp_lib.h" ! the openMP function library
integer(pInt), intent(in) :: & ! according to MSC.Marc 2012 Manual D
integer, intent(in) :: & ! according to MSC.Marc 2012 Manual D
ngens, & !< size of stress-strain law
nn, & !< integration point number
ndi, & !< number of direct components
@ -193,7 +189,7 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, &
jtype, & !< element type
ifr, & !< set to 1 if R has been calculated
ifu !< set to 1 if stretch has been calculated
integer(pInt), dimension(2), intent(in) :: & ! according to MSC.Marc 2012 Manual D
integer, dimension(2), intent(in) :: & ! according to MSC.Marc 2012 Manual D
m, & !< (1) user element number, (2) internal element number
matus, & !< (1) user material identification number, (2) internal material identification number
kcus, & !< (1) layer number, (2) internal layer number
@ -236,10 +232,10 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, &
logical :: cutBack
real(pReal), dimension(6) :: stress
real(pReal), dimension(6,6) :: ddsdde
integer(pInt) :: computationMode, i, cp_en, node, CPnodeID
integer :: computationMode, i, cp_en, node, CPnodeID
!$ integer(4) :: defaultNumThreadsInt !< default value set by Marc
if (iand(debug_level(debug_MARC),debug_LEVELBASIC) /= 0_pInt) then
if (iand(debug_level(debug_MARC),debug_LEVELBASIC) /= 0) then
write(6,'(a,/,i8,i8,i2)') ' MSC.MARC information on shape of element(2), IP:', m, nn
write(6,'(a,2(i1))') ' Jacobian: ', ngens,ngens
write(6,'(a,i1)') ' Direct stress: ', ndi
@ -260,7 +256,7 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, &
!$ call omp_set_num_threads(DAMASK_NumThreadsInt) ! set number of threads for parallel execution set by DAMASK_NUM_THREADS
computationMode = 0_pInt ! save initialization value, since it does not result in any calculation
computationMode = 0 ! save initialization value, since it does not result in any calculation
if (lovl == 4 ) then ! jacobian requested by marc
if (timinc < theDelta .and. theInc == inc .and. lastLovl /= lovl) & ! first after cutback
computationMode = CPFEM_RESTOREJACOBIAN
@ -307,7 +303,7 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, &
if (lastLovl /= lovl) then ! first after ping pong
call debug_reset() ! resets debugging
outdatedFFN1 = .false.
cycleCounter = cycleCounter + 1_pInt
cycleCounter = cycleCounter + 1
mesh_cellnode = mesh_build_cellnodes(mesh_node,mesh_Ncellnodes) ! update cell node coordinates
call mesh_build_ipCoordinates() ! update ip coordinates
endif
@ -324,7 +320,7 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, &
lastIncConverged = .false. ! reset flag
endif
do node = 1,theMesh%elem%nNodes
CPnodeID = mesh_element(4_pInt+node,cp_en)
CPnodeID = mesh_element(4+node,cp_en)
mesh_node(1:ndeg,CPnodeID) = mesh_node0(1:ndeg,CPnodeID) + numerics_unitlength * dispt(1:ndeg,node)
enddo
endif
@ -336,7 +332,7 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, &
call debug_info() ! first reports (meaningful) debugging
call debug_reset() ! and resets debugging
outdatedFFN1 = .false.
cycleCounter = cycleCounter + 1_pInt
cycleCounter = cycleCounter + 1
mesh_cellnode = mesh_build_cellnodes(mesh_node,mesh_Ncellnodes) ! update cell node coordinates
call mesh_build_ipCoordinates() ! update ip coordinates
endif
@ -376,9 +372,7 @@ end subroutine hypela2
!> @brief calculate internal heat generated due to inelastic energy dissipation
!--------------------------------------------------------------------------------------------------
subroutine flux(f,ts,n,time)
use prec, only: &
pReal, &
pInt
use prec
use thermal_conduction, only: &
thermal_conduction_getSourceAndItsTangent
use mesh, only: &
@ -387,7 +381,7 @@ subroutine flux(f,ts,n,time)
implicit none
real(pReal), dimension(6), intent(in) :: &
ts
integer(pInt), dimension(10), intent(in) :: &
integer, dimension(10), intent(in) :: &
n
real(pReal), intent(in) :: &
time
@ -404,9 +398,7 @@ subroutine flux(f,ts,n,time)
!> @details select a variable contour plotting (user subroutine).
!--------------------------------------------------------------------------------------------------
subroutine uedinc(inc,incsub)
use prec, only: &
pReal, &
pInt
use prec
use CPFEM, only: &
CPFEM_results
@ -424,9 +416,7 @@ end subroutine uedinc
!> @details select a variable contour plotting (user subroutine).
!--------------------------------------------------------------------------------------------------
subroutine plotv(v,s,sp,etot,eplas,ecreep,t,m,nn,layer,ndi,nshear,jpltcd)
use prec, only: &
pReal, &
pInt
use prec
use mesh, only: &
mesh_FEasCP
use IO, only: &
@ -436,7 +426,7 @@ subroutine plotv(v,s,sp,etot,eplas,ecreep,t,m,nn,layer,ndi,nshear,jpltcd)
materialpoint_sizeResults
implicit none
integer(pInt), intent(in) :: &
integer, intent(in) :: &
m, & !< element number
nn, & !< integration point number
layer, & !< layer number
@ -453,7 +443,7 @@ subroutine plotv(v,s,sp,etot,eplas,ecreep,t,m,nn,layer,ndi,nshear,jpltcd)
real(pReal), intent(out) :: &
v !< variable
if (jpltcd > materialpoint_sizeResults) call IO_error(700_pInt,jpltcd) ! complain about out of bounds error
if (jpltcd > materialpoint_sizeResults) call IO_error(700,jpltcd) ! complain about out of bounds error
v = materialpoint_results(jpltcd,nn,mesh_FEasCP('elem', m))
end subroutine plotv

57
src/FEsolving.f90
View File

@ -5,14 +5,12 @@
!> @todo Descriptions for public variables needed
!--------------------------------------------------------------------------------------------------
module FEsolving
use prec, only: &
pInt, &
pReal
use prec
implicit none
private
integer(pInt), public :: &
restartInc = 1_pInt !< needs description
integer, public :: &
restartInc = 1 !< needs description
logical, public :: &
symmetricSolver = .false., & !< use a symmetric FEM solver
@ -20,10 +18,10 @@ module FEsolving
restartRead = .false., & !< restart information to continue calculation from saved state
terminallyIll = .false. !< at least one material point is terminally ill
integer(pInt), dimension(:,:), allocatable, public :: &
integer, dimension(:,:), allocatable, public :: &
FEsolving_execIP !< for ping-pong scheme always range to max IP, otherwise one specific IP
integer(pInt), dimension(2), public :: &
integer, dimension(2), public :: &
FEsolving_execElem !< for ping-pong scheme always whole range, otherwise one specific element
character(len=1024), public :: &
@ -59,13 +57,12 @@ subroutine FE_init
IO_warning
use DAMASK_interface
implicit none
#if defined(Marc4DAMASK) || defined(Abaqus)
integer(pInt), parameter :: &
FILEUNIT = 222_pInt
integer(pInt) :: j
integer, parameter :: &
FILEUNIT = 222
integer :: j
character(len=65536) :: tag, line
integer(pInt), allocatable, dimension(:) :: chunkPos
integer, allocatable, dimension(:) :: chunkPos
#endif
write(6,'(/,a)') ' <<<+- FEsolving init -+>>>'
@ -75,35 +72,35 @@ subroutine FE_init
#if defined(Grid) || defined(FEM)
restartInc = interface_RestartInc
if(restartInc < 0_pInt) then
call IO_warning(warning_ID=34_pInt)
restartInc = 0_pInt
if(restartInc < 0) then
call IO_warning(warning_ID=34)
restartInc = 0
endif
restartRead = restartInc > 0_pInt ! only read in if "true" restart requested
restartRead = restartInc > 0 ! only read in if "true" restart requested
#else
call IO_open_inputFile(FILEUNIT,modelName)
rewind(FILEUNIT)
do
read (FILEUNIT,'(a1024)',END=100) line
chunkPos = IO_stringPos(line)
tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) ! extract key
tag = IO_lc(IO_stringValue(line,chunkPos,1)) ! extract key
select case(tag)
case ('solver')
read (FILEUNIT,'(a1024)',END=100) line ! next line
chunkPos = IO_stringPos(line)
symmetricSolver = (IO_intValue(line,chunkPos,2_pInt) /= 1_pInt)
symmetricSolver = (IO_intValue(line,chunkPos,2) /= 1)
case ('restart')
read (FILEUNIT,'(a1024)',END=100) line ! next line
chunkPos = IO_stringPos(line)
restartWrite = iand(IO_intValue(line,chunkPos,1_pInt),1_pInt) > 0_pInt
restartRead = iand(IO_intValue(line,chunkPos,1_pInt),2_pInt) > 0_pInt
restartWrite = iand(IO_intValue(line,chunkPos,1),1) > 0
restartRead = iand(IO_intValue(line,chunkPos,1),2) > 0
case ('*restart')
do j=2_pInt,chunkPos(1)
do j=2,chunkPos(1)
restartWrite = (IO_lc(IO_StringValue(line,chunkPos,j)) == 'write') .or. restartWrite
restartRead = (IO_lc(IO_StringValue(line,chunkPos,j)) == 'read') .or. restartRead
enddo
if(restartWrite) then
do j=2_pInt,chunkPos(1)
do j=2,chunkPos(1)
restartWrite = (IO_lc(IO_StringValue(line,chunkPos,j)) /= 'frequency=0') .and. restartWrite
enddo
endif
@ -118,11 +115,11 @@ subroutine FE_init
do
read (FILEUNIT,'(a1024)',END=200) line
chunkPos = IO_stringPos(line)
if ( IO_lc(IO_stringValue(line,chunkPos,1_pInt)) == 'restart' &
.and. IO_lc(IO_stringValue(line,chunkPos,2_pInt)) == 'file' &
.and. IO_lc(IO_stringValue(line,chunkPos,3_pInt)) == 'job' &
.and. IO_lc(IO_stringValue(line,chunkPos,4_pInt)) == 'id' ) &
modelName = IO_StringValue(line,chunkPos,6_pInt)
if ( IO_lc(IO_stringValue(line,chunkPos,1)) == 'restart' &
.and. IO_lc(IO_stringValue(line,chunkPos,2)) == 'file' &
.and. IO_lc(IO_stringValue(line,chunkPos,3)) == 'job' &
.and. IO_lc(IO_stringValue(line,chunkPos,4)) == 'id' ) &
modelName = IO_StringValue(line,chunkPos,6)
enddo
#else ! QUESTION: is this meaningful for the spectral/FEM case?
call IO_open_inputFile(FILEUNIT,modelName)
@ -130,10 +127,10 @@ subroutine FE_init
do
read (FILEUNIT,'(a1024)',END=200) line
chunkPos = IO_stringPos(line)
if (IO_lc(IO_stringValue(line,chunkPos,1_pInt))=='*heading') then
if (IO_lc(IO_stringValue(line,chunkPos,1))=='*heading') then
read (FILEUNIT,'(a1024)',END=200) line
chunkPos = IO_stringPos(line)
modelName = IO_StringValue(line,chunkPos,1_pInt)
modelName = IO_StringValue(line,chunkPos,1)
endif
enddo
#endif
@ -141,7 +138,7 @@ subroutine FE_init
endif
#endif
if (iand(debug_level(debug_FEsolving),debug_levelBasic) /= 0_pInt) then
if (iand(debug_level(debug_FEsolving),debug_levelBasic) /= 0) then
write(6,'(a21,l1)') ' restart writing: ', restartWrite
write(6,'(a21,l1)') ' restart reading: ', restartRead
if (restartRead) write(6,'(a,/)') ' restart Job: '//trim(modelName)

511
src/IO.f90
View File

File diff suppressed because it is too large Load Diff

19
src/config.f90
View File

@ -6,13 +6,11 @@
!! parts 'homogenization', 'crystallite', 'phase', 'texture', and 'microstucture'
!--------------------------------------------------------------------------------------------------
module config
use prec, only: &
pReal
use list, only: &
tPartitionedStringList
use prec
use list
implicit none
private
type(tPartitionedStringList), public, protected, allocatable, dimension(:) :: &
config_phase, &
config_microstructure, &
@ -47,8 +45,6 @@ contains
!> @brief reads material.config and stores its content per part
!--------------------------------------------------------------------------------------------------
subroutine config_init
use prec, only: &
pStringLen
use DAMASK_interface, only: &
getSolverJobName
use IO, only: &
@ -61,7 +57,6 @@ subroutine config_init
debug_material, &
debug_levelBasic
implicit none
integer :: myDebug,i
character(len=pStringLen) :: &
@ -149,7 +144,6 @@ recursive function read_materialConfig(fileName,cnt) result(fileContent)
use IO, only: &
IO_warning
implicit none
character(len=*), intent(in) :: fileName
integer, intent(in), optional :: cnt !< recursion counter
character(len=pStringLen), dimension(:), allocatable :: fileContent !< file content, separated per lines
@ -231,12 +225,10 @@ end function read_materialConfig
!--------------------------------------------------------------------------------------------------
subroutine parse_materialConfig(sectionNames,part,line, &
fileContent)
use prec, only: &
pStringLen
use IO, only: &
IO_intOut
implicit none
character(len=64), allocatable, dimension(:), intent(out) :: sectionNames
type(tPartitionedStringList), allocatable, dimension(:), intent(inout) :: part
character(len=pStringLen), intent(inout) :: line
@ -288,7 +280,7 @@ end subroutine parse_materialConfig
!--------------------------------------------------------------------------------------------------
subroutine parse_debugAndNumericsConfig(config_list, &
fileContent)
implicit none
type(tPartitionedStringList), intent(out) :: config_list
character(len=pStringLen), dimension(:), intent(in) :: fileContent
integer :: i
@ -309,7 +301,6 @@ subroutine config_deallocate(what)
use IO, only: &
IO_error
implicit none
character(len=*), intent(in) :: what
select case(trim(what))

44
src/damage_local.f90
View File

@ -3,19 +3,17 @@
!> @brief material subroutine for locally evolving damage field
!--------------------------------------------------------------------------------------------------
module damage_local
use prec, only: &
pReal, &
pInt
use prec
implicit none
private
integer(pInt), dimension(:,:), allocatable, target, public :: &
integer, dimension(:,:), allocatable, target, public :: &
damage_local_sizePostResult !< size of each post result output
character(len=64), dimension(:,:), allocatable, target, public :: &
damage_local_output !< name of each post result output
integer(pInt), dimension(:), allocatable, target, public :: &
integer, dimension(:), allocatable, target, public :: &
damage_local_Noutput !< number of outputs per instance of this damage
enum, bind(c)
@ -62,11 +60,10 @@ subroutine damage_local_init
use config, only: &
config_homogenization
implicit none
integer(pInt) :: maxNinstance,homog,instance,o,i
integer(pInt) :: sizeState
integer(pInt) :: NofMyHomog, h
integer :: maxNinstance,homog,instance,o,i
integer :: sizeState
integer :: NofMyHomog, h
integer(kind(undefined_ID)) :: &
outputID
character(len=65536), dimension(0), parameter :: emptyStringArray = [character(len=65536)::]
@ -76,13 +73,13 @@ subroutine damage_local_init
write(6,'(/,a)') ' <<<+- damage_'//DAMAGE_local_label//' init -+>>>'
maxNinstance = int(count(damage_type == DAMAGE_local_ID),pInt)
if (maxNinstance == 0_pInt) return
if (maxNinstance == 0) return
allocate(damage_local_sizePostResult (maxval(homogenization_Noutput),maxNinstance),source=0_pInt)
allocate(damage_local_sizePostResult (maxval(homogenization_Noutput),maxNinstance),source=0)
allocate(damage_local_output (maxval(homogenization_Noutput),maxNinstance))
damage_local_output = ''
allocate(damage_local_outputID (maxval(homogenization_Noutput),maxNinstance),source=undefined_ID)
allocate(damage_local_Noutput (maxNinstance), source=0_pInt)
allocate(damage_local_Noutput (maxNinstance), source=0)
allocate(param(maxNinstance))
@ -116,7 +113,7 @@ subroutine damage_local_init
! allocate state arrays
sizeState = 1_pInt
sizeState = 1
damageState(homog)%sizeState = sizeState
damageState(homog)%sizePostResults = sum(damage_local_sizePostResult(:,instance))
allocate(damageState(homog)%state0 (sizeState,NofMyHomog), source=damage_initialPhi(homog))
@ -147,15 +144,14 @@ function damage_local_updateState(subdt, ip, el)
mappingHomogenization, &
damageState
implicit none
integer(pInt), intent(in) :: &
integer, intent(in) :: &
ip, & !< integration point number
el !< element number
real(pReal), intent(in) :: &
subdt
logical, dimension(2) :: &
damage_local_updateState
integer(pInt) :: &
integer :: &
homog, &
offset
real(pReal) :: &
@ -201,13 +197,12 @@ subroutine damage_local_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, el
use source_damage_anisoDuctile, only: &
source_damage_anisoductile_getRateAndItsTangent
implicit none
integer(pInt), intent(in) :: &
integer, intent(in) :: &
ip, & !< integration point number
el !< element number
real(pReal), intent(in) :: &
phi
integer(pInt) :: &
integer :: &
phase, &
grain, &
source, &
@ -259,27 +254,26 @@ function damage_local_postResults(ip,el)
damageMapping, &
damage
implicit none
integer(pInt), intent(in) :: &
integer, intent(in) :: &
ip, & !< integration point
el !< element
real(pReal), dimension(sum(damage_local_sizePostResult(:,damage_typeInstance(material_homogenizationAt(el))))) :: &
damage_local_postResults
integer(pInt) :: &
integer :: &
instance, homog, offset, o, c
homog = material_homogenizationAt(el)
offset = damageMapping(homog)%p(ip,el)
instance = damage_typeInstance(homog)
associate(prm => param(instance))
c = 0_pInt
c = 0
outputsLoop: do o = 1_pInt,size(prm%outputID)
outputsLoop: do o = 1,size(prm%outputID)
select case(prm%outputID(o))
case (damage_ID)
damage_local_postResults(c+1_pInt) = damage(homog)%p(offset)
damage_local_postResults(c+1) = damage(homog)%p(offset)
c = c + 1
end select
enddo outputsLoop

1
src/damage_none.f90
View File

@ -27,7 +27,6 @@ subroutine damage_none_init()
DAMAGE_NONE_LABEL, &
DAMAGE_NONE_ID
implicit none
integer :: &
homog, &
NofMyHomog

58
src/damage_nonlocal.f90
View File

@ -4,19 +4,17 @@
!> @details to be done
!--------------------------------------------------------------------------------------------------
module damage_nonlocal
use prec, only: &
pReal, &
pInt
use prec
implicit none
private
integer(pInt), dimension(:,:), allocatable, target, public :: &
integer, dimension(:,:), allocatable, target, public :: &
damage_nonlocal_sizePostResult !< size of each post result output
character(len=64), dimension(:,:), allocatable, target, public :: &
damage_nonlocal_output !< name of each post result output
integer(pInt), dimension(:), allocatable, target, public :: &
integer, dimension(:), allocatable, target, public :: &
damage_nonlocal_Noutput !< number of outputs per instance of this damage
enum, bind(c)
@ -62,11 +60,10 @@ subroutine damage_nonlocal_init
use config, only: &
config_homogenization
implicit none
integer(pInt) :: maxNinstance,homog,instance,o,i
integer(pInt) :: sizeState
integer(pInt) :: NofMyHomog, h
integer :: maxNinstance,homog,instance,o,i
integer :: sizeState
integer :: NofMyHomog, h
integer(kind(undefined_ID)) :: &
outputID
character(len=65536), dimension(0), parameter :: emptyStringArray = [character(len=65536)::]
@ -75,13 +72,13 @@ subroutine damage_nonlocal_init
write(6,'(/,a)') ' <<<+- damage_'//DAMAGE_nonlocal_label//' init -+>>>'
maxNinstance = int(count(damage_type == DAMAGE_nonlocal_ID),pInt)
if (maxNinstance == 0_pInt) return
maxNinstance = int(count(damage_type == DAMAGE_nonlocal_ID))
if (maxNinstance == 0) return
allocate(damage_nonlocal_sizePostResult (maxval(homogenization_Noutput),maxNinstance),source=0_pInt)
allocate(damage_nonlocal_sizePostResult (maxval(homogenization_Noutput),maxNinstance),source=0)
allocate(damage_nonlocal_output (maxval(homogenization_Noutput),maxNinstance))
damage_nonlocal_output = ''
allocate(damage_nonlocal_Noutput (maxNinstance), source=0_pInt)
allocate(damage_nonlocal_Noutput (maxNinstance), source=0)
allocate(param(maxNinstance))
@ -114,7 +111,7 @@ subroutine damage_nonlocal_init
! allocate state arrays
sizeState = 1_pInt
sizeState = 1
damageState(homog)%sizeState = sizeState
damageState(homog)%sizePostResults = sum(damage_nonlocal_sizePostResult(:,instance))
allocate(damageState(homog)%state0 (sizeState,NofMyHomog), source=damage_initialPhi(homog))
@ -154,13 +151,12 @@ subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip,
use source_damage_anisoDuctile, only: &
source_damage_anisoductile_getRateAndItsTangent
implicit none
integer(pInt), intent(in) :: &
integer, intent(in) :: &
ip, & !< integration point number
el !< element number
real(pReal), intent(in) :: &
phi
integer(pInt) :: &
integer :: &
phase, &
grain, &
source, &
@ -217,13 +213,12 @@ function damage_nonlocal_getDiffusion33(ip,el)
use crystallite, only: &
crystallite_push33ToRef
implicit none
integer(pInt), intent(in) :: &
integer, intent(in) :: &
ip, & !< integration point number
el !< element number
real(pReal), dimension(3,3) :: &
damage_nonlocal_getDiffusion33
integer(pInt) :: &
integer :: &
homog, &
grain
@ -235,7 +230,7 @@ function damage_nonlocal_getDiffusion33(ip,el)
enddo
damage_nonlocal_getDiffusion33 = &
charLength**2_pInt*damage_nonlocal_getDiffusion33/real(homogenization_Ngrains(homog),pReal)
charLength**2*damage_nonlocal_getDiffusion33/real(homogenization_Ngrains(homog),pReal)
end function damage_nonlocal_getDiffusion33
@ -251,11 +246,10 @@ real(pReal) function damage_nonlocal_getMobility(ip,el)
material_phase, &
homogenization_Ngrains
implicit none
integer(pInt), intent(in) :: &
integer, intent(in) :: &
ip, & !< integration point number
el !< element number
integer(pInt) :: &
integer :: &
ipc
damage_nonlocal_getMobility = 0.0_pReal
@ -278,13 +272,12 @@ subroutine damage_nonlocal_putNonLocalDamage(phi,ip,el)
damageMapping, &
damage
implicit none
integer(pInt), intent(in) :: &
integer, intent(in) :: &
ip, & !< integration point number
el !< element number
real(pReal), intent(in) :: &
phi
integer(pInt) :: &
integer :: &
homog, &
offset
@ -304,27 +297,26 @@ function damage_nonlocal_postResults(ip,el)
damageMapping, &
damage
implicit none
integer(pInt), intent(in) :: &
integer, intent(in) :: &
ip, & !< integration point
el !< element
real(pReal), dimension(sum(damage_nonlocal_sizePostResult(:,damage_typeInstance(material_homogenizationAt(el))))) :: &
damage_nonlocal_postResults
integer(pInt) :: &
integer :: &
instance, homog, offset, o, c
homog = material_homogenizationAt(el)
offset = damageMapping(homog)%p(ip,el)
instance = damage_typeInstance(homog)
associate(prm => param(instance))
c = 0_pInt
c = 0
outputsLoop: do o = 1_pInt,size(prm%outputID)
outputsLoop: do o = 1,size(prm%outputID)
select case(prm%outputID(o))
case (damage_ID)
damage_nonlocal_postResults(c+1_pInt) = damage(homog)%p(offset)
damage_nonlocal_postResults(c+1) = damage(homog)%p(offset)
c = c + 1
end select
enddo outputsLoop

4
src/element.f90
View File

@ -3,8 +3,7 @@
!> @author Christoph Koords, Max-Planck-Institut für Eisenforschung GmbH
!--------------------------------------------------------------------------------------------------
module element
use prec, only: &
pReal
use prec
implicit none
private
@ -802,7 +801,6 @@ module element
use IO, only: &
IO_error
implicit none
class(tElement) :: self
integer, intent(in) :: elemType
self%elemType = elemType

274
src/grid/DAMASK_grid.f90
View File

@ -9,69 +9,19 @@
program DAMASK_spectral
#include <petsc/finclude/petscsys.h>
use PETScsys
use prec, only: &
pInt, &
pLongInt, &
pReal, &
tol_math_check, &
dNeq
use DAMASK_interface, only: &
DAMASK_interface_init, &
loadCaseFile, &
geometryFile, &
getSolverJobName, &
interface_restartInc
use IO, only: &
IO_isBlank, &
IO_stringPos, &
IO_stringValue, &
IO_floatValue, &
IO_intValue, &
IO_error, &
IO_lc, &
IO_intOut, &
IO_warning
use config, only: &
config_numerics
use debug, only: &
debug_level, &
debug_spectral, &
debug_levelBasic
use math ! need to include the whole module for FFTW
use mesh, only: &
grid, &
geomSize
use CPFEM2, only: &
CPFEM_initAll, &
CPFEM_results
use FEsolving, only: &
restartWrite, &
restartInc
use numerics, only: &
worldrank, &
worldsize, &
stagItMax, &
maxCutBack, &
continueCalculation
use homogenization, only: &
materialpoint_sizeResults, &
materialpoint_results, &
materialpoint_postResults
use material, only: &
thermal_type, &
damage_type, &
THERMAL_conduction_ID, &
DAMAGE_nonlocal_ID
use spectral_utilities, only: &
utilities_init, &
tSolutionState, &
tLoadCase, &
cutBack, &
nActiveFields, &
FIELD_UNDEFINED_ID, &
FIELD_MECH_ID, &
FIELD_THERMAL_ID, &
FIELD_DAMAGE_ID
use prec
use DAMASK_interface
use IO
use config
use debug
use math
use mesh
use CPFEM2
use FEsolving
use numerics
use homogenization
use material
use spectral_utilities
use grid_mech_spectral_basic
use grid_mech_spectral_polarisation
use grid_mech_FEM
@ -86,11 +36,11 @@ program DAMASK_spectral
! variables related to information from load case and geom file
real(pReal), dimension(9) :: temp_valueVector = 0.0_pReal !< temporarily from loadcase file when reading in tensors (initialize to 0.0)
logical, dimension(9) :: temp_maskVector = .false. !< temporarily from loadcase file when reading in tensors
integer(pInt), allocatable, dimension(:) :: chunkPos
integer(pInt) :: &
N_t = 0_pInt, & !< # of time indicators found in load case file
N_n = 0_pInt, & !< # of increment specifiers found in load case file
N_def = 0_pInt !< # of rate of deformation specifiers found in load case file
integer, allocatable, dimension(:) :: chunkPos
integer :: &
N_t = 0, & !< # of time indicators found in load case file
N_n = 0, & !< # of increment specifiers found in load case file
N_def = 0 !< # of rate of deformation specifiers found in load case file
character(len=65536) :: &
line
@ -99,8 +49,8 @@ program DAMASK_spectral
real(pReal), dimension(3,3), parameter :: &
ones = 1.0_pReal, &
zeros = 0.0_pReal
integer(pInt), parameter :: &
subStepFactor = 2_pInt !< for each substep, divide the last time increment by 2.0
integer, parameter :: &
subStepFactor = 2 !< for each substep, divide the last time increment by 2.0
real(pReal) :: &
time = 0.0_pReal, & !< elapsed time
time0 = 0.0_pReal, & !< begin of interval
@ -110,21 +60,21 @@ program DAMASK_spectral
logical :: &
guess, & !< guess along former trajectory
stagIterate
integer(pInt) :: &
integer :: &
i, j, k, l, field, &
errorID = 0_pInt, &
cutBackLevel = 0_pInt, & !< cut back level \f$ t = \frac{t_{inc}}{2^l} \f$
stepFraction = 0_pInt !< fraction of current time interval
integer(pInt) :: &
currentLoadcase = 0_pInt, & !< current load case
errorID = 0, &
cutBackLevel = 0, & !< cut back level \f$ t = \frac{t_{inc}}{2^l} \f$
stepFraction = 0 !< fraction of current time interval
integer :: &
currentLoadcase = 0, & !< current load case
inc, & !< current increment in current load case
totalIncsCounter = 0_pInt, & !< total # of increments
convergedCounter = 0_pInt, & !< # of converged increments
notConvergedCounter = 0_pInt, & !< # of non-converged increments
fileUnit = 0_pInt, & !< file unit for reading load case and writing results
totalIncsCounter = 0, & !< total # of increments
convergedCounter = 0, & !< # of converged increments
notConvergedCounter = 0, & !< # of non-converged increments
fileUnit = 0, & !< file unit for reading load case and writing results
myStat, &
statUnit = 0_pInt, & !< file unit for statistics output
lastRestartWritten = 0_pInt, & !< total increment # at which last restart information was written
statUnit = 0, & !< file unit for statistics output
lastRestartWritten = 0, & !< total increment # at which last restart information was written
stagIter
character(len=6) :: loadcase_string
character(len=1024) :: &
@ -134,8 +84,8 @@ program DAMASK_spectral
type(tSolutionState), allocatable, dimension(:) :: solres
integer(MPI_OFFSET_KIND) :: fileOffset
integer(MPI_OFFSET_KIND), dimension(:), allocatable :: outputSize
integer(pInt), parameter :: maxByteOut = 2147483647-4096 !< limit of one file output write https://trac.mpich.org/projects/mpich/ticket/1742
integer(pInt), parameter :: maxRealOut = maxByteOut/pReal
integer, parameter :: maxByteOut = 2147483647-4096 !< limit of one file output write https://trac.mpich.org/projects/mpich/ticket/1742
integer, parameter :: maxRealOut = maxByteOut/pReal
integer(pLongInt), dimension(2) :: outputIndex
PetscErrorCode :: ierr
procedure(grid_mech_spectral_basic_init), pointer :: &
@ -174,20 +124,20 @@ program DAMASK_spectral
case ('polarisation')
if(iand(debug_level(debug_spectral),debug_levelBasic)/= 0) &
call IO_warning(42_pInt, ext_msg='debug Divergence')
call IO_warning(42, ext_msg='debug Divergence')
mech_init => grid_mech_spectral_polarisation_init
mech_forward => grid_mech_spectral_polarisation_forward
mech_solution => grid_mech_spectral_polarisation_solution
case ('fem')
if(iand(debug_level(debug_spectral),debug_levelBasic)/= 0) &
call IO_warning(42_pInt, ext_msg='debug Divergence')
call IO_warning(42, ext_msg='debug Divergence')
mech_init => grid_mech_FEM_init
mech_forward => grid_mech_FEM_forward
mech_solution => grid_mech_FEM_solution
case default
call IO_error(error_ID = 891_pInt, ext_msg = config_numerics%getString('spectral_solver'))
call IO_error(error_ID = 891, ext_msg = config_numerics%getString('spectral_solver'))
end select
@ -195,27 +145,27 @@ program DAMASK_spectral
! reading information from load case file and to sanity checks
allocate (loadCases(0)) ! array of load cases
open(newunit=fileunit,iostat=myStat,file=trim(loadCaseFile),action='read')
if (myStat /= 0_pInt) call IO_error(100_pInt,el=myStat,ext_msg=trim(loadCaseFile))
if (myStat /= 0) call IO_error(100,el=myStat,ext_msg=trim(loadCaseFile))
do
read(fileUnit, '(A)', iostat=myStat) line
if ( myStat /= 0_pInt) exit
if ( myStat /= 0) exit
if (IO_isBlank(line)) cycle ! skip empty lines
currentLoadCase = currentLoadCase + 1_pInt
currentLoadCase = currentLoadCase + 1
chunkPos = IO_stringPos(line)
do i = 1_pInt, chunkPos(1) ! reading compulsory parameters for loadcase
do i = 1, chunkPos(1) ! reading compulsory parameters for loadcase
select case (IO_lc(IO_stringValue(line,chunkPos,i)))
case('l','velocitygrad','velgrad','velocitygradient','fdot','dotf','f')
N_def = N_def + 1_pInt
N_def = N_def + 1
case('t','time','delta')
N_t = N_t + 1_pInt
N_t = N_t + 1
case('n','incs','increments','steps','logincs','logincrements','logsteps')
N_n = N_n + 1_pInt
N_n = N_n + 1
end select
enddo
if ((N_def /= N_n) .or. (N_n /= N_t) .or. N_n < 1_pInt) & ! sanity check
call IO_error(error_ID=837_pInt,el=currentLoadCase,ext_msg = trim(loadCaseFile)) ! error message for incomplete loadcase
if ((N_def /= N_n) .or. (N_n /= N_t) .or. N_n < 1) & ! sanity check
call IO_error(error_ID=837,el=currentLoadCase,ext_msg = trim(loadCaseFile)) ! error message for incomplete loadcase
newLoadCase%stress%myType='stress'
field = 1
@ -229,7 +179,7 @@ program DAMASK_spectral
newLoadCase%ID(field) = FIELD_DAMAGE_ID
endif damageActive
readIn: do i = 1_pInt, chunkPos(1)
readIn: do i = 1, chunkPos(1)
select case (IO_lc(IO_stringValue(line,chunkPos,i)))
case('fdot','dotf','l','velocitygrad','velgrad','velocitygradient','f') ! assign values for the deformation BC matrix
temp_valueVector = 0.0_pReal
@ -241,7 +191,7 @@ program DAMASK_spectral
else
newLoadCase%deformation%myType = 'l'
endif
do j = 1_pInt, 9_pInt
do j = 1, 9
temp_maskVector(j) = IO_stringValue(line,chunkPos,i+j) /= '*' ! true if not a *
if (temp_maskVector(j)) temp_valueVector(j) = IO_floatValue(line,chunkPos,i+j) ! read value where applicable
enddo
@ -250,7 +200,7 @@ program DAMASK_spectral
newLoadCase%deformation%values = math_9to33(temp_valueVector) ! values in 3x3 notation
case('p','pk1','piolakirchhoff','stress', 's')
temp_valueVector = 0.0_pReal
do j = 1_pInt, 9_pInt
do j = 1, 9
temp_maskVector(j) = IO_stringValue(line,chunkPos,i+j) /= '*' ! true if not an asterisk
if (temp_maskVector(j)) temp_valueVector(j) = IO_floatValue(line,chunkPos,i+j) ! read value where applicable
enddo
@ -258,54 +208,54 @@ program DAMASK_spectral
newLoadCase%stress%maskFloat = merge(ones,zeros,newLoadCase%stress%maskLogical)
newLoadCase%stress%values = math_9to33(temp_valueVector)
case('t','time','delta') ! increment time
newLoadCase%time = IO_floatValue(line,chunkPos,i+1_pInt)
newLoadCase%time = IO_floatValue(line,chunkPos,i+1)
case('n','incs','increments','steps') ! number of increments
newLoadCase%incs = IO_intValue(line,chunkPos,i+1_pInt)
newLoadCase%incs = IO_intValue(line,chunkPos,i+1)
case('logincs','logincrements','logsteps') ! number of increments (switch to log time scaling)
newLoadCase%incs = IO_intValue(line,chunkPos,i+1_pInt)
newLoadCase%logscale = 1_pInt
newLoadCase%incs = IO_intValue(line,chunkPos,i+1)
newLoadCase%logscale = 1
case('freq','frequency','outputfreq') ! frequency of result writings
newLoadCase%outputfrequency = IO_intValue(line,chunkPos,i+1_pInt)
newLoadCase%outputfrequency = IO_intValue(line,chunkPos,i+1)
case('r','restart','restartwrite') ! frequency of writing restart information
newLoadCase%restartfrequency = &
max(0_pInt,IO_intValue(line,chunkPos,i+1_pInt))
max(0,IO_intValue(line,chunkPos,i+1))
case('guessreset','dropguessing')
newLoadCase%followFormerTrajectory = .false. ! do not continue to predict deformation along former trajectory
case('euler') ! rotation of load case given in euler angles
temp_valueVector = 0.0_pReal
l = 1_pInt ! assuming values given in degrees
k = 1_pInt ! assuming keyword indicating degree/radians present
select case (IO_lc(IO_stringValue(line,chunkPos,i+1_pInt)))
l = 1 ! assuming values given in degrees
k = 1 ! assuming keyword indicating degree/radians present
select case (IO_lc(IO_stringValue(line,chunkPos,i+1)))
case('deg','degree')
case('rad','radian') ! don't convert from degree to radian
l = 0_pInt
l = 0
case default
k = 0_pInt
k = 0
end select
do j = 1_pInt, 3_pInt
do j = 1, 3
temp_valueVector(j) = IO_floatValue(line,chunkPos,i+k+j)
enddo
if (l == 1_pInt) temp_valueVector(1:3) = temp_valueVector(1:3) * inRad ! convert to rad
if (l == 1) temp_valueVector(1:3) = temp_valueVector(1:3) * INRAD ! convert to rad
newLoadCase%rotation = math_EulerToR(temp_valueVector(1:3)) ! convert rad Eulers to rotation matrix
case('rotation','rot') ! assign values for the rotation matrix
temp_valueVector = 0.0_pReal
do j = 1_pInt, 9_pInt
do j = 1, 9
temp_valueVector(j) = IO_floatValue(line,chunkPos,i+j)
enddo
newLoadCase%rotation = math_9to33(temp_valueVector)
end select
enddo readIn
newLoadCase%followFormerTrajectory = merge(.true.,.false.,currentLoadCase > 1_pInt) ! by default, guess from previous load case
newLoadCase%followFormerTrajectory = merge(.true.,.false.,currentLoadCase > 1) ! by default, guess from previous load case
reportAndCheck: if (worldrank == 0) then
write (loadcase_string, '(i6)' ) currentLoadCase
write(6,'(/,1x,a,i6)') 'load case: ', currentLoadCase
if (.not. newLoadCase%followFormerTrajectory) write(6,'(2x,a)') 'drop guessing along trajectory'
if (newLoadCase%deformation%myType == 'l') then
do j = 1_pInt, 3_pInt
do j = 1, 3
if (any(newLoadCase%deformation%maskLogical(j,1:3) .eqv. .true.) .and. &
any(newLoadCase%deformation%maskLogical(j,1:3) .eqv. .false.)) errorID = 832_pInt ! each row should be either fully or not at all defined
any(newLoadCase%deformation%maskLogical(j,1:3) .eqv. .false.)) errorID = 832 ! each row should be either fully or not at all defined
enddo
write(6,'(2x,a)') 'velocity gradient:'
else if (newLoadCase%deformation%myType == 'f') then
@ -313,7 +263,7 @@ program DAMASK_spectral
else
write(6,'(2x,a)') 'deformation gradient rate:'
endif
do i = 1_pInt, 3_pInt; do j = 1_pInt, 3_pInt
do i = 1, 3; do j = 1, 3
if(newLoadCase%deformation%maskLogical(i,j)) then
write(6,'(2x,f12.7)',advance='no') newLoadCase%deformation%values(i,j)
else
@ -322,13 +272,13 @@ program DAMASK_spectral
enddo; write(6,'(/)',advance='no')
enddo
if (any(newLoadCase%stress%maskLogical .eqv. &
newLoadCase%deformation%maskLogical)) errorID = 831_pInt ! exclusive or masking only
newLoadCase%deformation%maskLogical)) errorID = 831 ! exclusive or masking only
if (any(newLoadCase%stress%maskLogical .and. &
transpose(newLoadCase%stress%maskLogical) .and. &
reshape([ .false.,.true.,.true.,.true.,.false.,.true.,.true.,.true.,.false.],[ 3,3]))) &
errorID = 838_pInt ! no rotation is allowed by stress BC
errorID = 838 ! no rotation is allowed by stress BC
write(6,'(2x,a)') 'stress / GPa:'
do i = 1_pInt, 3_pInt; do j = 1_pInt, 3_pInt
do i = 1, 3; do j = 1, 3
if(newLoadCase%stress%maskLogical(i,j)) then
write(6,'(2x,f12.7)',advance='no') newLoadCase%stress%values(i,j)*1e-9_pReal
else
@ -340,18 +290,18 @@ program DAMASK_spectral
transpose(newLoadCase%rotation))-math_I3) > &
reshape(spread(tol_math_check,1,9),[ 3,3]))&
.or. abs(math_det33(newLoadCase%rotation)) > &
1.0_pReal + tol_math_check) errorID = 846_pInt ! given rotation matrix contains strain
1.0_pReal + tol_math_check) errorID = 846 ! given rotation matrix contains strain
if (any(dNeq(newLoadCase%rotation, math_I3))) &
write(6,'(2x,a,/,3(3(3x,f12.7,1x)/))',advance='no') 'rotation of loadframe:',&
transpose(newLoadCase%rotation)
if (newLoadCase%time < 0.0_pReal) errorID = 834_pInt ! negative time increment
if (newLoadCase%time < 0.0_pReal) errorID = 834 ! negative time increment
write(6,'(2x,a,f12.6)') 'time: ', newLoadCase%time
if (newLoadCase%incs < 1_pInt) errorID = 835_pInt ! non-positive incs count
if (newLoadCase%incs < 1) errorID = 835 ! non-positive incs count
write(6,'(2x,a,i5)') 'increments: ', newLoadCase%incs
if (newLoadCase%outputfrequency < 1_pInt) errorID = 836_pInt ! non-positive result frequency
if (newLoadCase%outputfrequency < 1) errorID = 836 ! non-positive result frequency
write(6,'(2x,a,i5)') 'output frequency: ', newLoadCase%outputfrequency
write(6,'(2x,a,i5)') 'restart frequency: ', newLoadCase%restartfrequency
if (errorID > 0_pInt) call IO_error(error_ID = errorID, ext_msg = loadcase_string) ! exit with error message
if (errorID > 0) call IO_error(error_ID = errorID, ext_msg = loadcase_string) ! exit with error message
endif reportAndCheck
loadCases = [loadCases,newLoadCase] ! load case is ok, append it
enddo
@ -383,7 +333,7 @@ program DAMASK_spectral
!--------------------------------------------------------------------------------------------------
! write header of output file
if (worldrank == 0) then
writeHeader: if (interface_restartInc < 1_pInt) then
writeHeader: if (interface_restartInc < 1) then
open(newunit=fileUnit,file=trim(getSolverJobName())//&
'.spectralOut',form='UNFORMATTED',status='REPLACE')
write(fileUnit) 'load:', trim(loadCaseFile) ! ... and write header
@ -417,59 +367,59 @@ program DAMASK_spectral
allocate(outputSize(worldsize), source = 0_MPI_OFFSET_KIND)
outputSize(worldrank+1) = size(materialpoint_results,kind=MPI_OFFSET_KIND)*int(pReal,MPI_OFFSET_KIND)
call MPI_allreduce(MPI_IN_PLACE,outputSize,worldsize,MPI_LONG,MPI_SUM,PETSC_COMM_WORLD,ierr) ! get total output size over each process
if (ierr /= 0_pInt) call IO_error(error_ID=894_pInt, ext_msg='MPI_allreduce')
if (ierr /= 0) call IO_error(error_ID=894, ext_msg='MPI_allreduce')
call MPI_file_open(PETSC_COMM_WORLD, trim(getSolverJobName())//'.spectralOut', &
MPI_MODE_WRONLY + MPI_MODE_APPEND, &
MPI_INFO_NULL, &
fileUnit, &
ierr)
if (ierr /= 0_pInt) call IO_error(error_ID=894_pInt, ext_msg='MPI_file_open')
if (ierr /= 0) call IO_error(error_ID=894, ext_msg='MPI_file_open')
call MPI_file_get_position(fileUnit,fileOffset,ierr) ! get offset from header
if (ierr /= 0_pInt) call IO_error(error_ID=894_pInt, ext_msg='MPI_file_get_position')
if (ierr /= 0) call IO_error(error_ID=894, ext_msg='MPI_file_get_position')
fileOffset = fileOffset + sum(outputSize(1:worldrank)) ! offset of my process in file (header + processes before me)
call MPI_file_seek (fileUnit,fileOffset,MPI_SEEK_SET,ierr)
if (ierr /= 0_pInt) call IO_error(error_ID=894_pInt, ext_msg='MPI_file_seek')
if (ierr /= 0) call IO_error(error_ID=894, ext_msg='MPI_file_seek')
writeUndeformed: if (interface_restartInc < 1_pInt) then
writeUndeformed: if (interface_restartInc < 1) then
write(6,'(1/,a)') ' ... writing initial configuration to file ........................'
call CPFEM_results(0_pInt,0.0_pReal)
call CPFEM_results(0,0.0_pReal)
do i = 1, size(materialpoint_results,3)/(maxByteOut/(materialpoint_sizeResults*pReal))+1 ! slice the output of my process in chunks not exceeding the limit for one output
outputIndex = int([(i-1_pInt)*((maxRealOut)/materialpoint_sizeResults)+1_pInt, & ! QUESTION: why not starting i at 0 instead of murky 1?
outputIndex = int([(i-1)*((maxRealOut)/materialpoint_sizeResults)+1, &
min(i*((maxRealOut)/materialpoint_sizeResults),size(materialpoint_results,3))],pLongInt)
call MPI_file_write(fileUnit,reshape(materialpoint_results(:,:,outputIndex(1):outputIndex(2)), &
[(outputIndex(2)-outputIndex(1)+1)*int(materialpoint_sizeResults,pLongInt)]), &
int((outputIndex(2)-outputIndex(1)+1)*int(materialpoint_sizeResults,pLongInt)), &
MPI_DOUBLE, MPI_STATUS_IGNORE, ierr)
if (ierr /= 0_pInt) call IO_error(error_ID=894_pInt, ext_msg='MPI_file_write')
if (ierr /= 0) call IO_error(error_ID=894, ext_msg='MPI_file_write')
enddo
fileOffset = fileOffset + sum(outputSize) ! forward to current file position
endif writeUndeformed
loadCaseLooping: do currentLoadCase = 1_pInt, size(loadCases)
loadCaseLooping: do currentLoadCase = 1, size(loadCases)
time0 = time ! load case start time
guess = loadCases(currentLoadCase)%followFormerTrajectory ! change of load case? homogeneous guess for the first inc
incLooping: do inc = 1_pInt, loadCases(currentLoadCase)%incs
totalIncsCounter = totalIncsCounter + 1_pInt
incLooping: do inc = 1, loadCases(currentLoadCase)%incs
totalIncsCounter = totalIncsCounter + 1
!--------------------------------------------------------------------------------------------------
! forwarding time
timeIncOld = timeinc ! last timeinc that brought former inc to an end
if (loadCases(currentLoadCase)%logscale == 0_pInt) then ! linear scale
if (loadCases(currentLoadCase)%logscale == 0) then ! linear scale
timeinc = loadCases(currentLoadCase)%time/real(loadCases(currentLoadCase)%incs,pReal)
else
if (currentLoadCase == 1_pInt) then ! 1st load case of logarithmic scale
if (inc == 1_pInt) then ! 1st inc of 1st load case of logarithmic scale
timeinc = loadCases(1)%time*(2.0_pReal**real( 1_pInt-loadCases(1)%incs ,pReal)) ! assume 1st inc is equal to 2nd
if (currentLoadCase == 1) then ! 1st load case of logarithmic scale
if (inc == 1) then ! 1st inc of 1st load case of logarithmic scale
timeinc = loadCases(1)%time*(2.0_pReal**real( 1-loadCases(1)%incs ,pReal)) ! assume 1st inc is equal to 2nd
else ! not-1st inc of 1st load case of logarithmic scale
timeinc = loadCases(1)%time*(2.0_pReal**real(inc-1_pInt-loadCases(1)%incs ,pReal))
timeinc = loadCases(1)%time*(2.0_pReal**real(inc-1-loadCases(1)%incs ,pReal))
endif
else ! not-1st load case of logarithmic scale
timeinc = time0 * &
( (1.0_pReal + loadCases(currentLoadCase)%time/time0 )**(real( inc ,pReal)/&
real(loadCases(currentLoadCase)%incs ,pReal))&
-(1.0_pReal + loadCases(currentLoadCase)%time/time0 )**(real( inc-1_pInt ,pReal)/&
-(1.0_pReal + loadCases(currentLoadCase)%time/time0 )**(real( inc-1 ,pReal)/&
real(loadCases(currentLoadCase)%incs ,pReal)))
endif
endif
@ -479,12 +429,12 @@ program DAMASK_spectral
time = time + timeinc ! just advance time, skip already performed calculation
guess = .true. ! QUESTION:why forced guessing instead of inheriting loadcase preference
else skipping
stepFraction = 0_pInt ! fraction scaled by stepFactor**cutLevel
stepFraction = 0 ! fraction scaled by stepFactor**cutLevel
subStepLooping: do while (stepFraction < subStepFactor**cutBackLevel)
remainingLoadCaseTime = loadCases(currentLoadCase)%time+time0 - time
time = time + timeinc ! forward target time
stepFraction = stepFraction + 1_pInt ! count step
stepFraction = stepFraction + 1 ! count step
!--------------------------------------------------------------------------------------------------
! report begin of new step
@ -524,7 +474,7 @@ program DAMASK_spectral
!--------------------------------------------------------------------------------------------------
! solve fields
stagIter = 0_pInt
stagIter = 0
stagIterate = .true.
do while (stagIterate)
do field = 1, nActiveFields
@ -546,7 +496,7 @@ program DAMASK_spectral
if (.not. solres(field)%converged) exit ! no solution found
enddo
stagIter = stagIter + 1_pInt
stagIter = stagIter + 1
stagIterate = stagIter < stagItMax &
.and. all(solres(:)%converged) &
.and. .not. all(solres(:)%stagConverged) ! stationary with respect to staggered iteration
@ -567,52 +517,52 @@ program DAMASK_spectral
endif
elseif (cutBackLevel < maxCutBack) then ! further cutbacking tolerated?
cutBack = .true.
stepFraction = (stepFraction - 1_pInt) * subStepFactor ! adjust to new denominator
cutBackLevel = cutBackLevel + 1_pInt
stepFraction = (stepFraction - 1) * subStepFactor ! adjust to new denominator
cutBackLevel = cutBackLevel + 1
time = time - timeinc ! rewind time
timeinc = timeinc/real(subStepFactor,pReal) ! cut timestep
write(6,'(/,a)') ' cutting back '
else ! no more options to continue
call IO_warning(850_pInt)
call IO_warning(850)
call MPI_file_close(fileUnit,ierr)
close(statUnit)
call quit(-1_pInt*(lastRestartWritten+1_pInt)) ! quit and provide information about last restart inc written
call quit(-1*(lastRestartWritten+1)) ! quit and provide information about last restart inc written
endif
enddo subStepLooping
cutBackLevel = max(0_pInt, cutBackLevel - 1_pInt) ! try half number of subincs next inc
cutBackLevel = max(0, cutBackLevel - 1) ! try half number of subincs next inc
if (all(solres(:)%converged)) then
convergedCounter = convergedCounter + 1_pInt
convergedCounter = convergedCounter + 1
write(6,'(/,a,'//IO_intOut(totalIncsCounter)//',a)') & ! report converged inc
' increment ', totalIncsCounter, ' converged'
else
notConvergedCounter = notConvergedCounter + 1_pInt
notConvergedCounter = notConvergedCounter + 1
write(6,'(/,a,'//IO_intOut(totalIncsCounter)//',a)') & ! report non-converged inc
' increment ', totalIncsCounter, ' NOT converged'
endif; flush(6)
if (mod(inc,loadCases(currentLoadCase)%outputFrequency) == 0_pInt) then ! at output frequency
if (mod(inc,loadCases(currentLoadCase)%outputFrequency) == 0) then ! at output frequency
write(6,'(1/,a)') ' ... writing results to file ......................................'
flush(6)
call materialpoint_postResults()
call MPI_file_seek (fileUnit,fileOffset,MPI_SEEK_SET,ierr)
if (ierr /= 0_pInt) call IO_error(894_pInt, ext_msg='MPI_file_seek')
if (ierr /= 0) call IO_error(894, ext_msg='MPI_file_seek')
do i=1, size(materialpoint_results,3)/(maxByteOut/(materialpoint_sizeResults*pReal))+1 ! slice the output of my process in chunks not exceeding the limit for one output
outputIndex=int([(i-1_pInt)*((maxRealOut)/materialpoint_sizeResults)+1_pInt, &
outputIndex=int([(i-1)*((maxRealOut)/materialpoint_sizeResults)+1, &
min(i*((maxRealOut)/materialpoint_sizeResults),size(materialpoint_results,3))],pLongInt)
call MPI_file_write(fileUnit,reshape(materialpoint_results(:,:,outputIndex(1):outputIndex(2)),&
[(outputIndex(2)-outputIndex(1)+1)*int(materialpoint_sizeResults,pLongInt)]), &
int((outputIndex(2)-outputIndex(1)+1)*int(materialpoint_sizeResults,pLongInt)),&
MPI_DOUBLE, MPI_STATUS_IGNORE, ierr)
if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='MPI_file_write')
if(ierr /=0) call IO_error(894, ext_msg='MPI_file_write')
enddo
fileOffset = fileOffset + sum(outputSize) ! forward to current file position
call CPFEM_results(totalIncsCounter,time)
endif
if ( loadCases(currentLoadCase)%restartFrequency > 0_pInt & ! writing of restart info requested ...
.and. mod(inc,loadCases(currentLoadCase)%restartFrequency) == 0_pInt) then ! ... and at frequency of writing restart information
if ( loadCases(currentLoadCase)%restartFrequency > 0 & ! writing of restart info requested ...
.and. mod(inc,loadCases(currentLoadCase)%restartFrequency) == 0) then ! ... and at frequency of writing restart information
restartWrite = .true. ! set restart parameter for FEsolving
lastRestartWritten = inc ! QUESTION: first call to CPFEM_general will write?
endif
@ -636,7 +586,7 @@ program DAMASK_spectral
call MPI_file_close(fileUnit,ierr)
close(statUnit)
if (notConvergedCounter > 0_pInt) call quit(2_pInt) ! error if some are not converged
call quit(0_pInt) ! no complains ;)
if (notConvergedCounter > 0) call quit(2) ! error if some are not converged
call quit(0) ! no complains ;)
end program DAMASK_spectral

10
src/homogenization.f90
View File

@ -5,8 +5,7 @@
!> @brief homogenization manager, organizing deformation partitioning and stress homogenization
!--------------------------------------------------------------------------------------------------
module homogenization
use prec, only: &
pReal
use prec
use material
!--------------------------------------------------------------------------------------------------
@ -113,7 +112,6 @@ subroutine homogenization_init
use numerics, only: &
worldrank
implicit none
integer, parameter :: FILEUNIT = 200
integer :: e,i,p
integer, dimension(:,:), pointer :: thisSize
@ -351,7 +349,6 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
debug_i
#endif
implicit none
real(pReal), intent(in) :: dt !< time increment
logical, intent(in) :: updateJaco !< initiating Jacobian update
integer :: &
@ -653,7 +650,6 @@ subroutine materialpoint_postResults
crystallite_sizePostResults, &
crystallite_postResults
implicit none
integer :: &
thePos, &
theSize, &
@ -709,7 +705,6 @@ subroutine partitionDeformation(ip,el)
use homogenization_mech_RGC, only: &
homogenization_RGC_partitionDeformation
implicit none
integer, intent(in) :: &
ip, & !< integration point
el !< element number
@ -754,7 +749,6 @@ function updateState(ip,el)
use damage_local, only: &
damage_local_updateState
implicit none
integer, intent(in) :: &
ip, & !< integration point
el !< element number
@ -807,7 +801,6 @@ subroutine averageStressAndItsTangent(ip,el)
use homogenization_mech_RGC, only: &
homogenization_RGC_averageStressAndItsTangent
implicit none
integer, intent(in) :: &
ip, & !< integration point
el !< element number
@ -855,7 +848,6 @@ function postResults(ip,el)
use damage_nonlocal, only: &
damage_nonlocal_postResults
implicit none
integer, intent(in) :: &
ip, & !< integration point
el !< element number

23
src/homogenization_mech_RGC.f90
View File

@ -7,8 +7,7 @@
!> Nconstituents is defined as p x q x r (cluster)
!--------------------------------------------------------------------------------------------------
module homogenization_mech_RGC
use prec, only: &
pReal
use prec
use material
implicit none
@ -109,7 +108,6 @@ subroutine homogenization_RGC_init()
use config, only: &
config_homogenization
implicit none
integer :: &
Ninstance, &
h, i, &
@ -251,7 +249,6 @@ subroutine homogenization_RGC_partitionDeformation(F,avgF,instance,of)
debug_levelExtensive
#endif
implicit none
real(pReal), dimension (:,:,:), intent(out) :: F !< partioned F per grain
real(pReal), dimension (:,:), intent(in) :: avgF !< averaged F
@ -302,8 +299,6 @@ end subroutine homogenization_RGC_partitionDeformation
! "happy" with result
!--------------------------------------------------------------------------------------------------
function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el)
use prec, only: &
dEq0
#ifdef DEBUG
use debug, only: &
debug_level, &
@ -323,8 +318,6 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el)
viscModus_RGC, &
refRelaxRate_RGC
implicit none
real(pReal), dimension(:,:,:), intent(in) :: &
P,& !< array of P
F,& !< array of F
@ -748,7 +741,6 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el)
use numerics, only: &
xSmoo_RGC
implicit none
real(pReal), dimension (:,:,:), intent(out) :: rPen !< stress-like penalty
real(pReal), dimension (:,:), intent(out) :: nMis !< total amount of mismatch
@ -868,7 +860,6 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el)
volDiscrMod_RGC,&
volDiscrPow_RGC
implicit none
real(pReal), dimension (:,:,:), intent(out) :: vPen ! stress-like penalty due to volume
real(pReal), intent(out) :: vDiscrep ! total volume discrepancy
@ -919,7 +910,6 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el)
use math, only: &
math_invert33
implicit none
real(pReal), dimension(3) :: surfaceCorrection
real(pReal), dimension(3,3), intent(in) :: avgF !< average F
@ -953,7 +943,6 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el)
use constitutive, only: &
constitutive_homogenizedC
implicit none
real(pReal), dimension(2) :: equivalentModuli
integer, intent(in) :: &
@ -989,7 +978,6 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el)
!--------------------------------------------------------------------------------------------------
subroutine grainDeformation(F, avgF, instance, of)
implicit none
real(pReal), dimension(:,:,:), intent(out) :: F !< partioned F per grain
real(pReal), dimension(:,:), intent(in) :: avgF !< averaged F
@ -1032,7 +1020,6 @@ end function homogenization_RGC_updateState
!--------------------------------------------------------------------------------------------------
subroutine homogenization_RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,instance)
implicit none
real(pReal), dimension (3,3), intent(out) :: avgP !< average stress at material point
real(pReal), dimension (3,3,3,3), intent(out) :: dAvgPdAvgF !< average stiffness at material point
@ -1051,7 +1038,6 @@ end subroutine homogenization_RGC_averageStressAndItsTangent
!--------------------------------------------------------------------------------------------------
pure function homogenization_RGC_postResults(instance,of) result(postResults)
implicit none
integer, intent(in) :: &
instance, &
of
@ -1148,7 +1134,6 @@ end subroutine mech_RGC_results
!--------------------------------------------------------------------------------------------------
pure function relaxationVector(intFace,instance,of)
implicit none
real(pReal), dimension (3) :: relaxationVector
integer, intent(in) :: instance,of
@ -1176,7 +1161,6 @@ end function relaxationVector
!--------------------------------------------------------------------------------------------------
pure function interfaceNormal(intFace,instance,of)
implicit none
real(pReal), dimension(3) :: interfaceNormal
integer, dimension(4), intent(in) :: intFace !< interface ID in 4D array (normal and position)
@ -1202,7 +1186,6 @@ end function interfaceNormal
!--------------------------------------------------------------------------------------------------
pure function getInterface(iFace,iGrain3)
implicit none
integer, dimension(4) :: getInterface
integer, dimension(3), intent(in) :: iGrain3 !< grain ID in 3D array
@ -1227,7 +1210,6 @@ end function getInterface
!--------------------------------------------------------------------------------------------------
pure function grain1to3(grain1,nGDim)
implicit none
integer, dimension(3) :: grain1to3
integer, intent(in) :: grain1 !< grain ID in 1D array
@ -1245,7 +1227,6 @@ end function grain1to3
!--------------------------------------------------------------------------------------------------
integer pure function grain3to1(grain3,nGDim)
implicit none
integer, dimension(3), intent(in) :: grain3 !< grain ID in 3D array (pos.x,pos.y,pos.z)
integer, dimension(3), intent(in) :: nGDim
@ -1261,7 +1242,6 @@ end function grain3to1
!--------------------------------------------------------------------------------------------------
integer pure function interface4to1(iFace4D, nGDim)
implicit none
integer, dimension(4), intent(in) :: iFace4D !< interface ID in 4D array (n.dir,pos.x,pos.y,pos.z)
integer, dimension(3), intent(in) :: nGDim
@ -1308,7 +1288,6 @@ end function interface4to1
!--------------------------------------------------------------------------------------------------
pure function interface1to4(iFace1D, nGDim)
implicit none
integer, dimension(4) :: interface1to4
integer, intent(in) :: iFace1D !< interface ID in 1D array

3
src/homogenization_mech_isostrain.f90
View File

@ -39,7 +39,6 @@ module subroutine mech_isostrain_init
use config, only: &
config_homogenization
implicit none
integer :: &
Ninstance, &
h, &
@ -91,7 +90,6 @@ end subroutine mech_isostrain_init
!--------------------------------------------------------------------------------------------------
module subroutine mech_isostrain_partitionDeformation(F,avgF)
implicit none
real(pReal), dimension (:,:,:), intent(out) :: F !< partitioned deformation gradient
real(pReal), dimension (3,3), intent(in) :: avgF !< average deformation gradient at material point
@ -106,7 +104,6 @@ end subroutine mech_isostrain_partitionDeformation
!--------------------------------------------------------------------------------------------------
module subroutine mech_isostrain_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,instance)
implicit none
real(pReal), dimension (3,3), intent(out) :: avgP !< average stress at material point
real(pReal), dimension (3,3,3,3), intent(out) :: dAvgPdAvgF !< average stiffness at material point

1
src/homogenization_mech_none.f90
View File

@ -21,7 +21,6 @@ module subroutine mech_none_init
use config, only: &
config_homogenization
implicit none
integer :: &
Ninstance, &
h, &

62
src/kinematics_cleavage_opening.f90
View File

@ -5,18 +5,16 @@
!> @details to be done
!--------------------------------------------------------------------------------------------------
module kinematics_cleavage_opening
use prec, only: &
pReal, &
pInt
use prec
implicit none
private
integer(pInt), dimension(:), allocatable, private :: kinematics_cleavage_opening_instance
integer, dimension(:), allocatable, private :: kinematics_cleavage_opening_instance
type, private :: tParameters !< container type for internal constitutive parameters
integer(pInt) :: &
integer :: &
totalNcleavage
integer(pInt), dimension(:), allocatable :: &
integer, dimension(:), allocatable :: &
Ncleavage !< active number of cleavage systems per family
real(pReal) :: &
sdot0, &
@ -27,10 +25,10 @@ module kinematics_cleavage_opening
end type
! Begin Deprecated
integer(pInt), dimension(:), allocatable, private :: &
integer, dimension(:), allocatable, private :: &
kinematics_cleavage_opening_totalNcleavage !< total number of cleavage systems
integer(pInt), dimension(:,:), allocatable, private :: &
integer, dimension(:,:), allocatable, private :: &
kinematics_cleavage_opening_Ncleavage !< number of cleavage systems per family
real(pReal), dimension(:), allocatable, private :: &
@ -70,33 +68,32 @@ subroutine kinematics_cleavage_opening_init()
lattice_maxNcleavageFamily, &
lattice_NcleavageSystem
implicit none
integer(pInt), allocatable, dimension(:) :: tempInt
integer, allocatable, dimension(:) :: tempInt
real(pReal), allocatable, dimension(:) :: tempFloat
integer(pInt) :: maxNinstance,p,instance,kinematics
integer :: maxNinstance,p,instance,kinematics
write(6,'(/,a)') ' <<<+- kinematics_'//KINEMATICS_cleavage_opening_LABEL//' init -+>>>'
maxNinstance = int(count(phase_kinematics == KINEMATICS_cleavage_opening_ID),pInt)
if (maxNinstance == 0_pInt) return
maxNinstance = int(count(phase_kinematics == KINEMATICS_cleavage_opening_ID))
if (maxNinstance == 0) return
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) &
write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
allocate(kinematics_cleavage_opening_instance(size(config_phase)), source=0_pInt)
do p = 1_pInt, size(config_phase)
allocate(kinematics_cleavage_opening_instance(size(config_phase)), source=0)
do p = 1, size(config_phase)
kinematics_cleavage_opening_instance(p) = count(phase_kinematics(:,1:p) == kinematics_cleavage_opening_ID) ! ToDo: count correct?
enddo
allocate(kinematics_cleavage_opening_critDisp(lattice_maxNcleavageFamily,maxNinstance), source=0.0_pReal)
allocate(kinematics_cleavage_opening_critLoad(lattice_maxNcleavageFamily,maxNinstance), source=0.0_pReal)
allocate(kinematics_cleavage_opening_Ncleavage(lattice_maxNcleavageFamily,maxNinstance), source=0_pInt)
allocate(kinematics_cleavage_opening_totalNcleavage(maxNinstance), source=0_pInt)
allocate(kinematics_cleavage_opening_Ncleavage(lattice_maxNcleavageFamily,maxNinstance), source=0)
allocate(kinematics_cleavage_opening_totalNcleavage(maxNinstance), source=0)
allocate(kinematics_cleavage_opening_sdot_0(maxNinstance), source=0.0_pReal)
allocate(kinematics_cleavage_opening_N(maxNinstance), source=0.0_pReal)
do p = 1_pInt, size(config_phase)
do p = 1, size(config_phase)
if (all(phase_kinematics(:,p) /= KINEMATICS_cleavage_opening_ID)) cycle
instance = kinematics_cleavage_opening_instance(p)
kinematics_cleavage_opening_sdot_0(instance) = config_phase(p)%getFloat('anisobrittle_sdot0')
@ -115,13 +112,13 @@ subroutine kinematics_cleavage_opening_init()
kinematics_cleavage_opening_Ncleavage(1:lattice_maxNcleavageFamily,instance))
kinematics_cleavage_opening_totalNcleavage(instance) = sum(kinematics_cleavage_opening_Ncleavage(:,instance)) ! how many cleavage systems altogether
if (kinematics_cleavage_opening_sdot_0(instance) <= 0.0_pReal) &
call IO_error(211_pInt,el=instance,ext_msg='sdot_0 ('//KINEMATICS_cleavage_opening_LABEL//')')
call IO_error(211,el=instance,ext_msg='sdot_0 ('//KINEMATICS_cleavage_opening_LABEL//')')
if (any(kinematics_cleavage_opening_critDisp(1:size(tempInt),instance) < 0.0_pReal)) &
call IO_error(211_pInt,el=instance,ext_msg='critical_displacement ('//KINEMATICS_cleavage_opening_LABEL//')')
call IO_error(211,el=instance,ext_msg='critical_displacement ('//KINEMATICS_cleavage_opening_LABEL//')')
if (any(kinematics_cleavage_opening_critLoad(1:size(tempInt),instance) < 0.0_pReal)) &
call IO_error(211_pInt,el=instance,ext_msg='critical_load ('//KINEMATICS_cleavage_opening_LABEL//')')
call IO_error(211,el=instance,ext_msg='critical_load ('//KINEMATICS_cleavage_opening_LABEL//')')
if (kinematics_cleavage_opening_N(instance) <= 0.0_pReal) &
call IO_error(211_pInt,el=instance,ext_msg='rate_sensitivity ('//KINEMATICS_cleavage_opening_LABEL//')')
call IO_error(211,el=instance,ext_msg='rate_sensitivity ('//KINEMATICS_cleavage_opening_LABEL//')')
enddo
end subroutine kinematics_cleavage_opening_init
@ -130,8 +127,6 @@ end subroutine kinematics_cleavage_opening_init
!> @brief contains the constitutive equation for calculating the velocity gradient
!--------------------------------------------------------------------------------------------------
subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, ip, el)
use prec, only: &
tol_math_check
use math, only: &
math_mul33xx33
use material, only: &
@ -144,8 +139,7 @@ subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, i
lattice_maxNcleavageFamily, &
lattice_NcleavageSystem
implicit none
integer(pInt), intent(in) :: &
integer, intent(in) :: &
ipc, & !< grain number
ip, & !< integration point number
el !< element number
@ -155,7 +149,7 @@ subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, i
Ld !< damage velocity gradient
real(pReal), intent(out), dimension(3,3,3,3) :: &
dLd_dTstar !< derivative of Ld with respect to Tstar (4th-order tensor)
integer(pInt) :: &
integer :: &
instance, phase, &
homog, damageOffset, &
f, i, index_myFamily, k, l, m, n
@ -170,9 +164,9 @@ subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, i
Ld = 0.0_pReal
dLd_dTstar = 0.0_pReal
do f = 1_pInt,lattice_maxNcleavageFamily
index_myFamily = sum(lattice_NcleavageSystem(1:f-1_pInt,phase)) ! at which index starts my family
do i = 1_pInt,kinematics_cleavage_opening_Ncleavage(f,instance) ! process each (active) cleavage system in family
do f = 1,lattice_maxNcleavageFamily
index_myFamily = sum(lattice_NcleavageSystem(1:f-1,phase)) ! at which index starts my family
do i = 1,kinematics_cleavage_opening_Ncleavage(f,instance) ! process each (active) cleavage system in family
traction_d = math_mul33xx33(S,lattice_Scleavage(1:3,1:3,1,index_myFamily+i,phase))
traction_t = math_mul33xx33(S,lattice_Scleavage(1:3,1:3,2,index_myFamily+i,phase))
traction_n = math_mul33xx33(S,lattice_Scleavage(1:3,1:3,3,index_myFamily+i,phase))
@ -186,7 +180,7 @@ subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, i
Ld = Ld + udotd*lattice_Scleavage(1:3,1:3,1,index_myFamily+i,phase)
dudotd_dt = sign(1.0_pReal,traction_d)*udotd*kinematics_cleavage_opening_N(instance)/ &
max(0.0_pReal, abs(traction_d) - traction_crit)
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + &
dudotd_dt*lattice_Scleavage(k,l,1,index_myFamily+i,phase)* &
lattice_Scleavage(m,n,1,index_myFamily+i,phase)
@ -200,7 +194,7 @@ subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, i
Ld = Ld + udott*lattice_Scleavage(1:3,1:3,2,index_myFamily+i,phase)
dudott_dt = sign(1.0_pReal,traction_t)*udott*kinematics_cleavage_opening_N(instance)/ &
max(0.0_pReal, abs(traction_t) - traction_crit)
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + &
dudott_dt*lattice_Scleavage(k,l,2,index_myFamily+i,phase)* &
lattice_Scleavage(m,n,2,index_myFamily+i,phase)
@ -214,7 +208,7 @@ subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, i
Ld = Ld + udotn*lattice_Scleavage(1:3,1:3,3,index_myFamily+i,phase)
dudotn_dt = sign(1.0_pReal,traction_n)*udotn*kinematics_cleavage_opening_N(instance)/ &
max(0.0_pReal, abs(traction_n) - traction_crit)
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + &
dudotn_dt*lattice_Scleavage(k,l,3,index_myFamily+i,phase)* &
lattice_Scleavage(m,n,3,index_myFamily+i,phase)

36
src/kinematics_slipplane_opening.f90
View File

@ -5,18 +5,16 @@
!> @details to be done
!--------------------------------------------------------------------------------------------------
module kinematics_slipplane_opening
use prec, only: &
pReal, &
pInt
use prec
implicit none
private
integer(pInt), dimension(:), allocatable, private :: kinematics_slipplane_opening_instance
integer, dimension(:), allocatable, private :: kinematics_slipplane_opening_instance
type, private :: tParameters !< container type for internal constitutive parameters
integer(pInt) :: &
integer :: &
totalNslip
integer(pInt), dimension(:), allocatable :: &
integer, dimension(:), allocatable :: &
Nslip !< active number of slip systems per family
real(pReal) :: &
sdot0, &
@ -58,26 +56,25 @@ subroutine kinematics_slipplane_opening_init()
KINEMATICS_slipplane_opening_ID
use lattice
implicit none
integer(pInt) :: maxNinstance,p,instance,kinematics
integer :: maxNinstance,p,instance,kinematics
write(6,'(/,a)') ' <<<+- kinematics_'//KINEMATICS_slipplane_opening_LABEL//' init -+>>>'
maxNinstance = count(phase_kinematics == KINEMATICS_slipplane_opening_ID)
if (maxNinstance == 0) return
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) &
write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
allocate(kinematics_slipplane_opening_instance(size(config_phase)), source=0_pInt)
do p = 1_pInt, size(config_phase)
allocate(kinematics_slipplane_opening_instance(size(config_phase)), source=0)
do p = 1, size(config_phase)
kinematics_slipplane_opening_instance(p) = count(phase_kinematics(:,1:p) == kinematics_slipplane_opening_ID) ! ToDo: count correct?
enddo
allocate(param(maxNinstance))
do p = 1_pInt, size(config_phase)
do p = 1, size(config_phase)
if (all(phase_kinematics(:,p) /= KINEMATICS_slipplane_opening_ID)) cycle
associate(prm => param(kinematics_slipplane_opening_instance(p)), &
config => config_phase(p))
@ -99,11 +96,11 @@ prm%slip_direction = lattice_slip_direction (prm%Nslip,config%getString('latti
config%getFloat('c/a',defaultVal=0.0_pReal))
! if (kinematics_slipplane_opening_sdot_0(instance) <= 0.0_pReal) &
! call IO_error(211_pInt,el=instance,ext_msg='sdot_0 ('//KINEMATICS_slipplane_opening_LABEL//')')
! call IO_error(211,el=instance,ext_msg='sdot_0 ('//KINEMATICS_slipplane_opening_LABEL//')')
! if (any(kinematics_slipplane_opening_critPlasticStrain(:,instance) < 0.0_pReal)) &
! call IO_error(211_pInt,el=instance,ext_msg='criticaPlasticStrain ('//KINEMATICS_slipplane_opening_LABEL//')')
! call IO_error(211,el=instance,ext_msg='criticaPlasticStrain ('//KINEMATICS_slipplane_opening_LABEL//')')
! if (kinematics_slipplane_opening_N(instance) <= 0.0_pReal) &
! call IO_error(211_pInt,el=instance,ext_msg='rate_sensitivity ('//KINEMATICS_slipplane_opening_LABEL//')')
! call IO_error(211,el=instance,ext_msg='rate_sensitivity ('//KINEMATICS_slipplane_opening_LABEL//')')
end associate
enddo
@ -114,8 +111,6 @@ end subroutine kinematics_slipplane_opening_init
!> @brief contains the constitutive equation for calculating the velocity gradient
!--------------------------------------------------------------------------------------------------
subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, ip, el)
use prec, only: &
tol_math_check
use math, only: &
math_mul33xx33, &
math_outer
@ -125,7 +120,6 @@ subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc,
damage, &
damageMapping
implicit none
integer, intent(in) :: &
ipc, & !< grain number
ip, & !< integration point number
@ -173,7 +167,7 @@ subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc,
if (abs(udotd) > tol_math_check) then
Ld = Ld + udotd*projection_d
dudotd_dt = udotd*prm%n/traction_d
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + &
dudotd_dt*projection_d(k,l)*projection_d(m,n)
endif
@ -185,7 +179,7 @@ subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc,
if (abs(udott) > tol_math_check) then
Ld = Ld + udott*projection_t
dudott_dt = udott*prm%n/traction_t
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + &
dudott_dt*projection_t(k,l)*projection_t(m,n)
endif
@ -197,7 +191,7 @@ subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc,
if (abs(udotn) > tol_math_check) then
Ld = Ld + udotn*projection_n
dudotn_dt = udotn*prm%n/traction_n
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + &
dudotn_dt*projection_n(k,l)*projection_n(m,n)
endif

7
src/kinematics_thermal_expansion.f90
View File

@ -4,9 +4,7 @@
!> @details to be done
!--------------------------------------------------------------------------------------------------
module kinematics_thermal_expansion
use prec, only: &
pReal, &
pInt
use prec
implicit none
private
@ -42,7 +40,6 @@ subroutine kinematics_thermal_expansion_init()
use config, only: &
config_phase
implicit none
integer(pInt) :: &
Ninstance, &
p, i
@ -87,7 +84,6 @@ pure function kinematics_thermal_expansion_initialStrain(homog,phase,offset)
lattice_thermalExpansion33, &
lattice_referenceTemperature
implicit none
integer(pInt), intent(in) :: &
phase, &
homog, offset
@ -120,7 +116,6 @@ subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar, ipc, ip,
lattice_thermalExpansion33, &
lattice_referenceTemperature
implicit none
integer(pInt), intent(in) :: &
ipc, & !< grain number
ip, & !< integration point number

234
src/material.f90
View File

@ -98,10 +98,10 @@ module material
integer(kind(HOMOGENIZATION_undefined_ID)), dimension(:), allocatable, public, protected :: &
homogenization_type !< type of each homogenization
integer(pInt), public, protected :: &
integer, public, protected :: &
homogenization_maxNgrains !< max number of grains in any USED homogenization
integer(pInt), dimension(:), allocatable, public, protected :: &
integer, dimension(:), allocatable, public, protected :: &
phase_Nsources, & !< number of source mechanisms active in each phase
phase_Nkinematics, & !< number of kinematic mechanisms active in each phase
phase_NstiffnessDegradations, & !< number of stiffness degradation mechanisms active in each phase
@ -132,7 +132,7 @@ module material
! END NEW MAPPINGS
! DEPRECATED: use material_phaseAt
integer(pInt), dimension(:,:,:), allocatable, public :: &
integer, dimension(:,:,:), allocatable, public :: &
material_phase !< phase (index) of each grain,IP,element
type(tPlasticState), allocatable, dimension(:), public :: &
@ -144,7 +144,7 @@ module material
thermalState, &
damageState
integer(pInt), dimension(:,:,:), allocatable, public, protected :: &
integer, dimension(:,:,:), allocatable, public, protected :: &
material_texture !< texture (index) of each grain,IP,element
real(pReal), dimension(:,:,:,:), allocatable, public, protected :: &
@ -155,15 +155,15 @@ module material
microstructure_elemhomo, & !< flag to indicate homogeneous microstructure distribution over element's IPs
phase_localPlasticity !< flags phases with local constitutive law
integer(pInt), private :: &
integer, private :: &
microstructure_maxNconstituents, & !< max number of constituents in any phase
texture_maxNgauss !< max number of Gauss components in any texture
integer(pInt), dimension(:), allocatable, private :: &
integer, dimension(:), allocatable, private :: &
microstructure_Nconstituents, & !< number of constituents in each microstructure
texture_Ngauss !< number of Gauss components per texture
integer(pInt), dimension(:,:), allocatable, private :: &
integer, dimension(:,:), allocatable, private :: &
microstructure_phase, & !< phase IDs of each microstructure
microstructure_texture !< texture IDs of each microstructure
@ -178,11 +178,11 @@ module material
homogenization_active
! BEGIN DEPRECATED
integer(pInt), dimension(:,:,:), allocatable, public :: phaseAt !< phase ID of every material point (ipc,ip,el)
integer(pInt), dimension(:,:,:), allocatable, public :: phasememberAt !< memberID of given phase at every material point (ipc,ip,el)
integer, dimension(:,:,:), allocatable, public :: phaseAt !< phase ID of every material point (ipc,ip,el)
integer, dimension(:,:,:), allocatable, public :: phasememberAt !< memberID of given phase at every material point (ipc,ip,el)
integer(pInt), dimension(:,:,:), allocatable, public, target :: mappingHomogenization !< mapping from material points to offset in heterogenous state/field
integer(pInt), dimension(:,:), allocatable, private, target :: mappingHomogenizationConst !< mapping from material points to offset in constant state/field
integer, dimension(:,:,:), allocatable, public, target :: mappingHomogenization !< mapping from material points to offset in heterogenous state/field
integer, dimension(:,:), allocatable, private, target :: mappingHomogenizationConst !< mapping from material points to offset in constant state/field
! END DEPRECATED
type(tHomogMapping), allocatable, dimension(:), public :: &
@ -256,13 +256,13 @@ subroutine material_init
use mesh, only: &
theMesh
integer(pInt), parameter :: FILEUNIT = 210_pInt
integer(pInt) :: m,c,h, myDebug, myPhase, myHomog
integer(pInt) :: &
integer, parameter :: FILEUNIT = 210
integer :: m,c,h, myDebug, myPhase, myHomog
integer :: &
g, & !< grain number
i, & !< integration point number
e !< element number
integer(pInt), dimension(:), allocatable :: &
integer, dimension(:), allocatable :: &
CounterPhase, &
CounterHomogenization
@ -271,19 +271,19 @@ subroutine material_init
write(6,'(/,a)') ' <<<+- material init -+>>>'
call material_parsePhase()
if (iand(myDebug,debug_levelBasic) /= 0_pInt) write(6,'(a)') ' Phase parsed'; flush(6)
if (iand(myDebug,debug_levelBasic) /= 0) write(6,'(a)') ' Phase parsed'; flush(6)
call material_parseMicrostructure()
if (iand(myDebug,debug_levelBasic) /= 0_pInt) write(6,'(a)') ' Microstructure parsed'; flush(6)
if (iand(myDebug,debug_levelBasic) /= 0) write(6,'(a)') ' Microstructure parsed'; flush(6)
call material_parseCrystallite()
if (iand(myDebug,debug_levelBasic) /= 0_pInt) write(6,'(a)') ' Crystallite parsed'; flush(6)
if (iand(myDebug,debug_levelBasic) /= 0) write(6,'(a)') ' Crystallite parsed'; flush(6)
call material_parseHomogenization()
if (iand(myDebug,debug_levelBasic) /= 0_pInt) write(6,'(a)') ' Homogenization parsed'; flush(6)
if (iand(myDebug,debug_levelBasic) /= 0) write(6,'(a)') ' Homogenization parsed'; flush(6)
call material_parseTexture()
if (iand(myDebug,debug_levelBasic) /= 0_pInt) write(6,'(a)') ' Texture parsed'; flush(6)
if (iand(myDebug,debug_levelBasic) /= 0) write(6,'(a)') ' Texture parsed'; flush(6)
allocate(plasticState (size(config_phase)))
allocate(sourceState (size(config_phase)))
@ -303,34 +303,34 @@ subroutine material_init
allocate(temperatureRate (size(config_homogenization)))
do m = 1_pInt,size(config_microstructure)
if(microstructure_crystallite(m) < 1_pInt .or. &
do m = 1,size(config_microstructure)
if(microstructure_crystallite(m) < 1 .or. &
microstructure_crystallite(m) > size(config_crystallite)) &
call IO_error(150_pInt,m,ext_msg='crystallite')
if(minval(microstructure_phase(1:microstructure_Nconstituents(m),m)) < 1_pInt .or. &
call IO_error(150,m,ext_msg='crystallite')
if(minval(microstructure_phase(1:microstructure_Nconstituents(m),m)) < 1 .or. &
maxval(microstructure_phase(1:microstructure_Nconstituents(m),m)) > size(config_phase)) &
call IO_error(150_pInt,m,ext_msg='phase')
if(minval(microstructure_texture(1:microstructure_Nconstituents(m),m)) < 1_pInt .or. &
call IO_error(150,m,ext_msg='phase')
if(minval(microstructure_texture(1:microstructure_Nconstituents(m),m)) < 1 .or. &
maxval(microstructure_texture(1:microstructure_Nconstituents(m),m)) > size(config_texture)) &
call IO_error(150_pInt,m,ext_msg='texture')
if(microstructure_Nconstituents(m) < 1_pInt) &
call IO_error(151_pInt,m)
call IO_error(150,m,ext_msg='texture')
if(microstructure_Nconstituents(m) < 1) &
call IO_error(151,m)
enddo
debugOut: if (iand(myDebug,debug_levelExtensive) /= 0_pInt) then
debugOut: if (iand(myDebug,debug_levelExtensive) /= 0) then
write(6,'(/,a,/)') ' MATERIAL configuration'
write(6,'(a32,1x,a16,1x,a6)') 'homogenization ','type ','grains'
do h = 1_pInt,size(config_homogenization)
do h = 1,size(config_homogenization)
write(6,'(1x,a32,1x,a16,1x,i6)') homogenization_name(h),homogenization_type(h),homogenization_Ngrains(h)
enddo
write(6,'(/,a14,18x,1x,a11,1x,a12,1x,a13)') 'microstructure','crystallite','constituents','homogeneous'
do m = 1_pInt,size(config_microstructure)
do m = 1,size(config_microstructure)
write(6,'(1x,a32,1x,i11,1x,i12,1x,l13)') microstructure_name(m), &
microstructure_crystallite(m), &
microstructure_Nconstituents(m), &
microstructure_elemhomo(m)
if (microstructure_Nconstituents(m) > 0_pInt) then
do c = 1_pInt,microstructure_Nconstituents(m)
if (microstructure_Nconstituents(m) > 0) then
do c = 1,microstructure_Nconstituents(m)
write(6,'(a1,1x,a32,1x,a32,1x,f7.4)') '>',phase_name(microstructure_phase(c,m)),&
texture_name(microstructure_texture(c,m)),&
microstructure_fraction(c,m)
@ -383,23 +383,23 @@ subroutine material_init
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! BEGIN DEPRECATED
allocate(phaseAt ( homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%Nelems),source=0_pInt)
allocate(phasememberAt ( homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%Nelems),source=0_pInt)
allocate(mappingHomogenization (2, theMesh%elem%nIPs,theMesh%Nelems),source=0_pInt)
allocate(mappingHomogenizationConst( theMesh%elem%nIPs,theMesh%Nelems),source=1_pInt)
allocate(phaseAt ( homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%Nelems),source=0)
allocate(phasememberAt ( homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%Nelems),source=0)
allocate(mappingHomogenization (2, theMesh%elem%nIPs,theMesh%Nelems),source=0)
allocate(mappingHomogenizationConst( theMesh%elem%nIPs,theMesh%Nelems),source=1)
CounterHomogenization=0
CounterPhase =0
do e = 1_pInt,theMesh%Nelems
do e = 1,theMesh%Nelems
myHomog = theMesh%homogenizationAt(e)
do i = 1_pInt, theMesh%elem%nIPs
CounterHomogenization(myHomog) = CounterHomogenization(myHomog) + 1_pInt
do i = 1, theMesh%elem%nIPs
CounterHomogenization(myHomog) = CounterHomogenization(myHomog) + 1
mappingHomogenization(1:2,i,e) = [CounterHomogenization(myHomog),huge(1)]
do g = 1_pInt,homogenization_Ngrains(myHomog)
do g = 1,homogenization_Ngrains(myHomog)
myPhase = material_phase(g,i,e)
CounterPhase(myPhase) = CounterPhase(myPhase)+1_pInt ! not distinguishing between instances of same phase
CounterPhase(myPhase) = CounterPhase(myPhase)+1 ! not distinguishing between instances of same phase
phaseAt(g,i,e) = myPhase
phasememberAt(g,i,e) = CounterPhase(myPhase)
enddo
@ -429,33 +429,33 @@ subroutine material_parseHomogenization
use IO, only: &
IO_error
integer(pInt) :: h
integer :: h
character(len=65536) :: tag
allocate(homogenization_type(size(config_homogenization)), source=HOMOGENIZATION_undefined_ID)
allocate(thermal_type(size(config_homogenization)), source=THERMAL_isothermal_ID)
allocate(damage_type (size(config_homogenization)), source=DAMAGE_none_ID)
allocate(homogenization_typeInstance(size(config_homogenization)), source=0_pInt)
allocate(thermal_typeInstance(size(config_homogenization)), source=0_pInt)
allocate(damage_typeInstance(size(config_homogenization)), source=0_pInt)
allocate(homogenization_Ngrains(size(config_homogenization)), source=0_pInt)
allocate(homogenization_Noutput(size(config_homogenization)), source=0_pInt)
allocate(homogenization_typeInstance(size(config_homogenization)), source=0)
allocate(thermal_typeInstance(size(config_homogenization)), source=0)
allocate(damage_typeInstance(size(config_homogenization)), source=0)
allocate(homogenization_Ngrains(size(config_homogenization)), source=0)
allocate(homogenization_Noutput(size(config_homogenization)), source=0)
allocate(homogenization_active(size(config_homogenization)), source=.false.) !!!!!!!!!!!!!!!
allocate(thermal_initialT(size(config_homogenization)), source=300.0_pReal)
allocate(damage_initialPhi(size(config_homogenization)), source=1.0_pReal)
forall (h = 1_pInt:size(config_homogenization)) &
forall (h = 1:size(config_homogenization)) &
homogenization_active(h) = any(theMesh%homogenizationAt == h)
do h=1_pInt, size(config_homogenization)
do h=1, size(config_homogenization)
homogenization_Noutput(h) = config_homogenization(h)%countKeys('(output)')
tag = config_homogenization(h)%getString('mech')
select case (trim(tag))
case(HOMOGENIZATION_NONE_label)
homogenization_type(h) = HOMOGENIZATION_NONE_ID
homogenization_Ngrains(h) = 1_pInt
homogenization_Ngrains(h) = 1
case(HOMOGENIZATION_ISOSTRAIN_label)
homogenization_type(h) = HOMOGENIZATION_ISOSTRAIN_ID
homogenization_Ngrains(h) = config_homogenization(h)%getInt('nconstituents')
@ -463,7 +463,7 @@ subroutine material_parseHomogenization
homogenization_type(h) = HOMOGENIZATION_RGC_ID
homogenization_Ngrains(h) = config_homogenization(h)%getInt('nconstituents')
case default
call IO_error(500_pInt,ext_msg=trim(tag))
call IO_error(500,ext_msg=trim(tag))
end select
homogenization_typeInstance(h) = count(homogenization_type==homogenization_type(h))
@ -480,7 +480,7 @@ subroutine material_parseHomogenization
case(THERMAL_conduction_label)
thermal_type(h) = THERMAL_conduction_ID
case default
call IO_error(500_pInt,ext_msg=trim(tag))
call IO_error(500,ext_msg=trim(tag))
end select
endif
@ -497,14 +497,14 @@ subroutine material_parseHomogenization
case(DAMAGE_NONLOCAL_label)
damage_type(h) = DAMAGE_nonlocal_ID
case default
call IO_error(500_pInt,ext_msg=trim(tag))
call IO_error(500,ext_msg=trim(tag))
end select
endif
enddo
do h=1_pInt, size(config_homogenization)
do h=1, size(config_homogenization)
homogenization_typeInstance(h) = count(homogenization_type(1:h) == homogenization_type(h))
thermal_typeInstance(h) = count(thermal_type (1:h) == thermal_type (h))
damage_typeInstance(h) = count(damage_type (1:h) == damage_type (h))
@ -530,58 +530,58 @@ subroutine material_parseMicrostructure
character(len=65536), dimension(:), allocatable :: &
strings
integer(pInt), allocatable, dimension(:) :: chunkPos
integer(pInt) :: e, m, c, i
integer, allocatable, dimension(:) :: chunkPos
integer :: e, m, c, i
character(len=65536) :: &
tag
allocate(microstructure_crystallite(size(config_microstructure)), source=0_pInt)
allocate(microstructure_Nconstituents(size(config_microstructure)), source=0_pInt)
allocate(microstructure_crystallite(size(config_microstructure)), source=0)
allocate(microstructure_Nconstituents(size(config_microstructure)), source=0)
allocate(microstructure_active(size(config_microstructure)), source=.false.)
allocate(microstructure_elemhomo(size(config_microstructure)), source=.false.)
if(any(theMesh%microstructureAt > size(config_microstructure))) &
call IO_error(155_pInt,ext_msg='More microstructures in geometry than sections in material.config')
call IO_error(155,ext_msg='More microstructures in geometry than sections in material.config')
forall (e = 1_pInt:theMesh%Nelems) &
forall (e = 1:theMesh%Nelems) &
microstructure_active(theMesh%microstructureAt(e)) = .true. ! current microstructure used in model? Elementwise view, maximum N operations for N elements
do m=1_pInt, size(config_microstructure)
do m=1, size(config_microstructure)
microstructure_Nconstituents(m) = config_microstructure(m)%countKeys('(constituent)')
microstructure_crystallite(m) = config_microstructure(m)%getInt('crystallite')
microstructure_elemhomo(m) = config_microstructure(m)%keyExists('/elementhomogeneous/')
enddo
microstructure_maxNconstituents = maxval(microstructure_Nconstituents)
allocate(microstructure_phase (microstructure_maxNconstituents,size(config_microstructure)),source=0_pInt)
allocate(microstructure_texture (microstructure_maxNconstituents,size(config_microstructure)),source=0_pInt)
allocate(microstructure_phase (microstructure_maxNconstituents,size(config_microstructure)),source=0)
allocate(microstructure_texture (microstructure_maxNconstituents,size(config_microstructure)),source=0)
allocate(microstructure_fraction(microstructure_maxNconstituents,size(config_microstructure)),source=0.0_pReal)
allocate(strings(1)) ! Intel 16.0 Bug
do m=1_pInt, size(config_microstructure)
do m=1, size(config_microstructure)
strings = config_microstructure(m)%getStrings('(constituent)',raw=.true.)
do c = 1_pInt, size(strings)
do c = 1, size(strings)
chunkPos = IO_stringPos(strings(c))
do i = 1_pInt,5_pInt,2_pInt
do i = 1,5,2
tag = IO_stringValue(strings(c),chunkPos,i)
select case (tag)
case('phase')
microstructure_phase(c,m) = IO_intValue(strings(c),chunkPos,i+1_pInt)
microstructure_phase(c,m) = IO_intValue(strings(c),chunkPos,i+1)
case('texture')
microstructure_texture(c,m) = IO_intValue(strings(c),chunkPos,i+1_pInt)
microstructure_texture(c,m) = IO_intValue(strings(c),chunkPos,i+1)
case('fraction')
microstructure_fraction(c,m) = IO_floatValue(strings(c),chunkPos,i+1_pInt)
microstructure_fraction(c,m) = IO_floatValue(strings(c),chunkPos,i+1)
end select
enddo
enddo
enddo
do m = 1_pInt, size(config_microstructure)
do m = 1, size(config_microstructure)
if (dNeq(sum(microstructure_fraction(:,m)),1.0_pReal)) &
call IO_error(153_pInt,ext_msg=microstructure_name(m))
call IO_error(153,ext_msg=microstructure_name(m))
enddo
end subroutine material_parseMicrostructure
@ -592,10 +592,10 @@ end subroutine material_parseMicrostructure
!--------------------------------------------------------------------------------------------------
subroutine material_parseCrystallite
integer(pInt) :: c
integer :: c
allocate(crystallite_Noutput(size(config_crystallite)),source=0_pInt)
do c=1_pInt, size(config_crystallite)
allocate(crystallite_Noutput(size(config_crystallite)),source=0)
do c=1, size(config_crystallite)
crystallite_Noutput(c) = config_crystallite(c)%countKeys('(output)')
enddo
@ -611,19 +611,19 @@ subroutine material_parsePhase
IO_getTag, &
IO_stringValue
integer(pInt) :: sourceCtr, kinematicsCtr, stiffDegradationCtr, p
integer :: sourceCtr, kinematicsCtr, stiffDegradationCtr, p
character(len=65536), dimension(:), allocatable :: str
allocate(phase_elasticity(size(config_phase)),source=ELASTICITY_undefined_ID)
allocate(phase_plasticity(size(config_phase)),source=PLASTICITY_undefined_ID)
allocate(phase_Nsources(size(config_phase)), source=0_pInt)
allocate(phase_Nkinematics(size(config_phase)), source=0_pInt)
allocate(phase_NstiffnessDegradations(size(config_phase)),source=0_pInt)
allocate(phase_Noutput(size(config_phase)), source=0_pInt)
allocate(phase_Nsources(size(config_phase)), source=0)
allocate(phase_Nkinematics(size(config_phase)), source=0)
allocate(phase_NstiffnessDegradations(size(config_phase)),source=0)
allocate(phase_Noutput(size(config_phase)), source=0)
allocate(phase_localPlasticity(size(config_phase)), source=.false.)
do p=1_pInt, size(config_phase)
do p=1, size(config_phase)
phase_Noutput(p) = config_phase(p)%countKeys('(output)')
phase_Nsources(p) = config_phase(p)%countKeys('(source)')
phase_Nkinematics(p) = config_phase(p)%countKeys('(kinematics)')
@ -634,7 +634,7 @@ subroutine material_parsePhase
case (ELASTICITY_HOOKE_label)
phase_elasticity(p) = ELASTICITY_HOOKE_ID
case default
call IO_error(200_pInt,ext_msg=trim(config_phase(p)%getString('elasticity')))
call IO_error(200,ext_msg=trim(config_phase(p)%getString('elasticity')))
end select
select case (config_phase(p)%getString('plasticity'))
@ -653,7 +653,7 @@ subroutine material_parsePhase
case (PLASTICITY_NONLOCAL_label)
phase_plasticity(p) = PLASTICITY_NONLOCAL_ID
case default
call IO_error(201_pInt,ext_msg=trim(config_phase(p)%getString('plasticity')))
call IO_error(201,ext_msg=trim(config_phase(p)%getString('plasticity')))
end select
enddo
@ -662,7 +662,7 @@ subroutine material_parsePhase
allocate(phase_kinematics(maxval(phase_Nkinematics),size(config_phase)), source=KINEMATICS_undefined_ID)
allocate(phase_stiffnessDegradation(maxval(phase_NstiffnessDegradations),size(config_phase)), &
source=STIFFNESS_DEGRADATION_undefined_ID)
do p=1_pInt, size(config_phase)
do p=1, size(config_phase)
#if defined(__GFORTRAN__) || defined(__PGI)
str = ['GfortranBug86277']
str = config_phase(p)%getStrings('(source)',defaultVal=str)
@ -670,7 +670,7 @@ subroutine material_parsePhase
#else
str = config_phase(p)%getStrings('(source)',defaultVal=[character(len=65536)::])
#endif
do sourceCtr = 1_pInt, size(str)
do sourceCtr = 1, size(str)
select case (trim(str(sourceCtr)))
case (SOURCE_thermal_dissipation_label)
phase_source(sourceCtr,p) = SOURCE_thermal_dissipation_ID
@ -694,7 +694,7 @@ subroutine material_parsePhase
#else
str = config_phase(p)%getStrings('(kinematics)',defaultVal=[character(len=65536)::])
#endif
do kinematicsCtr = 1_pInt, size(str)
do kinematicsCtr = 1, size(str)
select case (trim(str(kinematicsCtr)))
case (KINEMATICS_cleavage_opening_label)
phase_kinematics(kinematicsCtr,p) = KINEMATICS_cleavage_opening_ID
@ -711,7 +711,7 @@ subroutine material_parsePhase
#else
str = config_phase(p)%getStrings('(stiffness_degradation)',defaultVal=[character(len=65536)::])
#endif
do stiffDegradationCtr = 1_pInt, size(str)
do stiffDegradationCtr = 1, size(str)
select case (trim(str(stiffDegradationCtr)))
case (STIFFNESS_DEGRADATION_damage_label)
phase_stiffnessDegradation(stiffDegradationCtr,p) = STIFFNESS_DEGRADATION_damage_ID
@ -719,10 +719,10 @@ subroutine material_parsePhase
enddo
enddo
allocate(phase_plasticityInstance(size(config_phase)), source=0_pInt)
allocate(phase_elasticityInstance(size(config_phase)), source=0_pInt)
allocate(phase_plasticityInstance(size(config_phase)), source=0)
allocate(phase_elasticityInstance(size(config_phase)), source=0)
do p=1_pInt, size(config_phase)
do p=1, size(config_phase)
phase_elasticityInstance(p) = count(phase_elasticity(1:p) == phase_elasticity(p))
phase_plasticityInstance(p) = count(phase_plasticity(1:p) == phase_plasticity(p))
enddo
@ -739,13 +739,13 @@ subroutine material_parseTexture
IO_floatValue, &
IO_stringValue
integer(pInt) :: section, gauss, j, t, i
integer :: section, gauss, j, t, i
character(len=65536), dimension(:), allocatable :: strings ! Values for given key in material config
integer(pInt), dimension(:), allocatable :: chunkPos
integer, dimension(:), allocatable :: chunkPos
allocate(texture_Ngauss(size(config_texture)), source=0_pInt)
allocate(texture_Ngauss(size(config_texture)), source=0)
do t=1_pInt, size(config_texture)
do t=1, size(config_texture)
texture_Ngauss(t) = config_texture(t)%countKeys('(gauss)')
if (config_texture(t)%keyExists('symmetry')) call IO_error(147,ext_msg='symmetry')
if (config_texture(t)%keyExists('(random)')) call IO_error(147,ext_msg='(random)')
@ -757,13 +757,13 @@ subroutine material_parseTexture
allocate(texture_transformation(3,3,size(config_texture)), source=0.0_pReal)
texture_transformation = spread(math_I3,3,size(config_texture))
do t=1_pInt, size(config_texture)
do t=1, size(config_texture)
section = t
gauss = 0_pInt
gauss = 0
if (config_texture(t)%keyExists('axes')) then
strings = config_texture(t)%getStrings('axes')
do j = 1_pInt, 3_pInt ! look for "x", "y", and "z" entries
do j = 1, 3 ! look for "x", "y", and "z" entries
select case (strings(j))
case('x', '+x')
texture_transformation(j,1:3,t) = [ 1.0_pReal, 0.0_pReal, 0.0_pReal] ! original axis is now +x-axis
@ -778,25 +778,25 @@ subroutine material_parseTexture
case('-z')
texture_transformation(j,1:3,t) = [ 0.0_pReal, 0.0_pReal,-1.0_pReal] ! original axis is now -z-axis
case default
call IO_error(157_pInt,t)
call IO_error(157,t)
end select
enddo
if(dNeq(math_det33(texture_transformation(1:3,1:3,t)),1.0_pReal)) call IO_error(157_pInt,t)
if(dNeq(math_det33(texture_transformation(1:3,1:3,t)),1.0_pReal)) call IO_error(157,t)
endif
if (config_texture(t)%keyExists('(gauss)')) then
gauss = gauss + 1_pInt
gauss = gauss + 1
strings = config_texture(t)%getStrings('(gauss)',raw= .true.)
do i = 1_pInt , size(strings)
do i = 1 , size(strings)
chunkPos = IO_stringPos(strings(i))
do j = 1_pInt,9_pInt,2_pInt
do j = 1,9,2
select case (IO_stringValue(strings(i),chunkPos,j))
case('phi1')
texture_Gauss(1,gauss,t) = IO_floatValue(strings(i),chunkPos,j+1_pInt)*inRad
texture_Gauss(1,gauss,t) = IO_floatValue(strings(i),chunkPos,j+1)*inRad
case('phi')
texture_Gauss(2,gauss,t) = IO_floatValue(strings(i),chunkPos,j+1_pInt)*inRad
texture_Gauss(2,gauss,t) = IO_floatValue(strings(i),chunkPos,j+1)*inRad
case('phi2')
texture_Gauss(3,gauss,t) = IO_floatValue(strings(i),chunkPos,j+1_pInt)*inRad
texture_Gauss(3,gauss,t) = IO_floatValue(strings(i),chunkPos,j+1)*inRad
end select
enddo
enddo
@ -817,7 +817,7 @@ subroutine material_allocatePlasticState(phase,NofMyPhase,&
use numerics, only: &
numerics_integrator
integer(pInt), intent(in) :: &
integer, intent(in) :: &
phase, &
NofMyPhase, &
sizeState, &
@ -842,13 +842,13 @@ subroutine material_allocatePlasticState(phase,NofMyPhase,&
allocate(plasticState(phase)%state (sizeState,NofMyPhase), source=0.0_pReal)
allocate(plasticState(phase)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal)
if (numerics_integrator == 1_pInt) then
if (numerics_integrator == 1) then
allocate(plasticState(phase)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal)
allocate(plasticState(phase)%previousDotState2 (sizeDotState,NofMyPhase), source=0.0_pReal)
endif
if (numerics_integrator == 4_pInt) &
if (numerics_integrator == 4) &
allocate(plasticState(phase)%RK4dotState (sizeDotState,NofMyPhase), source=0.0_pReal)
if (numerics_integrator == 5_pInt) &
if (numerics_integrator == 5) &
allocate(plasticState(phase)%RKCK45dotState (6,sizeDotState,NofMyPhase), source=0.0_pReal)
allocate(plasticState(phase)%deltaState (sizeDeltaState,NofMyPhase), source=0.0_pReal)
@ -864,7 +864,7 @@ subroutine material_allocateSourceState(phase,of,NofMyPhase,&
use numerics, only: &
numerics_integrator
integer(pInt), intent(in) :: &
integer, intent(in) :: &
phase, &
of, &
NofMyPhase, &
@ -882,13 +882,13 @@ subroutine material_allocateSourceState(phase,of,NofMyPhase,&
allocate(sourceState(phase)%p(of)%state (sizeState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(of)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal)
if (numerics_integrator == 1_pInt) then
if (numerics_integrator == 1) then
allocate(sourceState(phase)%p(of)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(of)%previousDotState2 (sizeDotState,NofMyPhase), source=0.0_pReal)
endif
if (numerics_integrator == 4_pInt) &
if (numerics_integrator == 4) &
allocate(sourceState(phase)%p(of)%RK4dotState (sizeDotState,NofMyPhase), source=0.0_pReal)
if (numerics_integrator == 5_pInt) &
if (numerics_integrator == 5) &
allocate(sourceState(phase)%p(of)%RKCK45dotState (6,sizeDotState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(of)%deltaState (sizeDeltaState,NofMyPhase), source=0.0_pReal)
@ -905,10 +905,10 @@ subroutine material_populateGrains
use mesh, only: &
theMesh
integer(pInt) :: e,i,c,homog,micro
integer :: e,i,c,homog,micro
allocate(material_phase(homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%Nelems), source=0_pInt)
allocate(material_texture(homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%Nelems), source=0_pInt)
allocate(material_phase(homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%Nelems), source=0)
allocate(material_texture(homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%Nelems), source=0)
allocate(material_EulerAngles(3,homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%Nelems),source=0.0_pReal)
do e = 1, theMesh%Nelems

1028
src/mesh_abaqus.f90
View File

File diff suppressed because it is too large Load Diff

10
src/mesh_base.f90
View File

@ -9,12 +9,8 @@
module mesh_base
use, intrinsic :: iso_c_binding
use prec, only: &
pStringLen, &
pReal, &
pInt
use element, only: &
tElement
use prec
use element
use future
implicit none
@ -54,7 +50,6 @@ module mesh_base
contains
subroutine tMesh_base_init(self,meshType,elemType,nodes)
implicit none
class(tMesh) :: self
character(len=*), intent(in) :: meshType
integer(pInt), intent(in) :: elemType
@ -75,7 +70,6 @@ end subroutine tMesh_base_init
subroutine tMesh_base_setNelems(self,Nelems)
implicit none
class(tMesh) :: self
integer(pInt), intent(in) :: Nelems

806
src/mesh_marc.f90
View File

File diff suppressed because it is too large Load Diff

7
src/numerics.f90
View File

@ -4,9 +4,7 @@
!> @brief Managing of parameters related to numerics
!--------------------------------------------------------------------------------------------------
module numerics
use prec, only: &
pInt, &
pReal
use prec
implicit none
private
@ -115,8 +113,6 @@ contains
! a sanity check
!--------------------------------------------------------------------------------------------------
subroutine numerics_init
use prec, only: &
pStringLen
use IO, only: &
IO_read_ASCII, &
IO_error, &
@ -132,7 +128,6 @@ subroutine numerics_init
use petscsys
#endif
!$ use OMP_LIB, only: omp_set_num_threads
implicit none
!$ integer :: gotDAMASK_NUM_THREADS = 1
integer :: i,j, ierr ! no pInt
integer(pInt), allocatable, dimension(:) :: chunkPos

74
src/source_damage_anisoBrittle.f90
View File

@ -5,23 +5,21 @@
!> @details to be done
!--------------------------------------------------------------------------------------------------
module source_damage_anisoBrittle
use prec, only: &
pReal, &
pInt
use prec
implicit none
private
integer(pInt), dimension(:), allocatable, public, protected :: &
integer, dimension(:), allocatable, public, protected :: &
source_damage_anisoBrittle_offset, & !< which source is my current source mechanism?
source_damage_anisoBrittle_instance !< instance of source mechanism
integer(pInt), dimension(:,:), allocatable, target, public :: &
integer, dimension(:,:), allocatable, target, public :: &
source_damage_anisoBrittle_sizePostResult !< size of each post result output
character(len=64), dimension(:,:), allocatable, target, public :: &
source_damage_anisoBrittle_output !< name of each post result output
integer(pInt), dimension(:,:), allocatable, private :: &
integer, dimension(:,:), allocatable, private :: &
source_damage_anisoBrittle_Ncleavage !< number of cleavage systems per family
enum, bind(c)
@ -40,9 +38,9 @@ module source_damage_anisoBrittle
critLoad
real(pReal), dimension(:,:,:,:), allocatable :: &
cleavage_systems
integer(pInt) :: &
integer :: &
totalNcleavage
integer(pInt), dimension(:), allocatable :: &
integer, dimension(:), allocatable :: &
Ncleavage
integer(kind(undefined_ID)), allocatable, dimension(:) :: &
outputID !< ID of each post result output
@ -65,8 +63,6 @@ contains
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
subroutine source_damage_anisoBrittle_init
use prec, only: &
pStringLen
use debug, only: &
debug_level,&
debug_constitutive,&
@ -91,9 +87,9 @@ subroutine source_damage_anisoBrittle_init
lattice_SchmidMatrix_cleavage, &
lattice_maxNcleavageFamily
integer(pInt) :: Ninstance,phase,instance,source,sourceOffset
integer(pInt) :: NofMyPhase,p ,i
integer(pInt), dimension(0), parameter :: emptyIntArray = [integer(pInt)::]
integer :: Ninstance,phase,instance,source,sourceOffset
integer :: NofMyPhase,p ,i
integer, dimension(0), parameter :: emptyIntArray = [integer::]
character(len=65536), dimension(0), parameter :: emptyStringArray = [character(len=65536)::]
integer(kind(undefined_ID)) :: &
outputID
@ -105,14 +101,14 @@ subroutine source_damage_anisoBrittle_init
write(6,'(/,a)') ' <<<+- source_'//SOURCE_DAMAGE_ANISOBRITTLE_LABEL//' init -+>>>'
Ninstance = int(count(phase_source == SOURCE_damage_anisoBrittle_ID),pInt)
if (Ninstance == 0_pInt) return
Ninstance = int(count(phase_source == SOURCE_damage_anisoBrittle_ID))
if (Ninstance == 0) return
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) &
write(6,'(a16,1x,i5,/)') '# instances:',Ninstance
allocate(source_damage_anisoBrittle_offset(material_Nphase), source=0_pInt)
allocate(source_damage_anisoBrittle_instance(material_Nphase), source=0_pInt)
allocate(source_damage_anisoBrittle_offset(material_Nphase), source=0)
allocate(source_damage_anisoBrittle_instance(material_Nphase), source=0)
do phase = 1, material_Nphase
source_damage_anisoBrittle_instance(phase) = count(phase_source(:,1:phase) == source_damage_anisoBrittle_ID)
do source = 1, phase_Nsources(phase)
@ -121,11 +117,11 @@ subroutine source_damage_anisoBrittle_init
enddo
enddo
allocate(source_damage_anisoBrittle_sizePostResult(maxval(phase_Noutput),Ninstance), source=0_pInt)
allocate(source_damage_anisoBrittle_sizePostResult(maxval(phase_Noutput),Ninstance), source=0)
allocate(source_damage_anisoBrittle_output(maxval(phase_Noutput),Ninstance))
source_damage_anisoBrittle_output = ''
allocate(source_damage_anisoBrittle_Ncleavage(lattice_maxNcleavageFamily,Ninstance), source=0_pInt)
allocate(source_damage_anisoBrittle_Ncleavage(lattice_maxNcleavageFamily,Ninstance), source=0)
allocate(param(Ninstance))
@ -162,18 +158,18 @@ subroutine source_damage_anisoBrittle_init
!--------------------------------------------------------------------------------------------------
! exit if any parameter is out of range
if (extmsg /= '') &
call IO_error(211_pInt,ext_msg=trim(extmsg)//'('//SOURCE_DAMAGE_ANISOBRITTLE_LABEL//')')
call IO_error(211,ext_msg=trim(extmsg)//'('//SOURCE_DAMAGE_ANISOBRITTLE_LABEL//')')
!--------------------------------------------------------------------------------------------------
! output pararameters
outputs = config%getStrings('(output)',defaultVal=emptyStringArray)
allocate(prm%outputID(0))
do i=1_pInt, size(outputs)
do i=1, size(outputs)
outputID = undefined_ID
select case(outputs(i))
case ('anisobrittle_drivingforce')
source_damage_anisoBrittle_sizePostResult(i,source_damage_anisoBrittle_instance(p)) = 1_pInt
source_damage_anisoBrittle_sizePostResult(i,source_damage_anisoBrittle_instance(p)) = 1
source_damage_anisoBrittle_output(i,source_damage_anisoBrittle_instance(p)) = outputs(i)
prm%outputID = [prm%outputID, damage_drivingforce_ID]
@ -189,7 +185,7 @@ subroutine source_damage_anisoBrittle_init
sourceOffset = source_damage_anisoBrittle_offset(phase)
call material_allocateSourceState(phase,sourceOffset,NofMyPhase,1_pInt,1_pInt,0_pInt)
call material_allocateSourceState(phase,sourceOffset,NofMyPhase,1,1,0)
sourceState(phase)%p(sourceOffset)%sizePostResults = sum(source_damage_anisoBrittle_sizePostResult(:,instance))
sourceState(phase)%p(sourceOffset)%aTolState=param(instance)%aTol
@ -217,13 +213,13 @@ subroutine source_damage_anisoBrittle_dotState(S, ipc, ip, el)
lattice_maxNcleavageFamily, &
lattice_NcleavageSystem
integer(pInt), intent(in) :: &
integer, intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
real(pReal), intent(in), dimension(3,3) :: &
S
integer(pInt) :: &
integer :: &
phase, &
constituent, &
instance, &
@ -243,10 +239,10 @@ subroutine source_damage_anisoBrittle_dotState(S, ipc, ip, el)
sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = 0.0_pReal
index = 1_pInt
do f = 1_pInt,lattice_maxNcleavageFamily
index_myFamily = sum(lattice_NcleavageSystem(1:f-1_pInt,phase)) ! at which index starts my family
do i = 1_pInt,source_damage_anisoBrittle_Ncleavage(f,instance) ! process each (active) cleavage system in family
index = 1
do f = 1,lattice_maxNcleavageFamily
index_myFamily = sum(lattice_NcleavageSystem(1:f-1,phase)) ! at which index starts my family
do i = 1,source_damage_anisoBrittle_Ncleavage(f,instance) ! process each (active) cleavage system in family
traction_d = math_mul33xx33(S,lattice_Scleavage(1:3,1:3,1,index_myFamily+i,phase))
traction_t = math_mul33xx33(S,lattice_Scleavage(1:3,1:3,2,index_myFamily+i,phase))
@ -263,7 +259,7 @@ subroutine source_damage_anisoBrittle_dotState(S, ipc, ip, el)
(max(0.0_pReal, abs(traction_n) - traction_crit)/traction_crit)**param(instance)%N)/ &
param(instance)%critDisp(index)
index = index + 1_pInt
index = index + 1
enddo
enddo
@ -276,7 +272,7 @@ subroutine source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalph
use material, only: &
sourceState
integer(pInt), intent(in) :: &
integer, intent(in) :: &
phase, &
constituent
real(pReal), intent(in) :: &
@ -284,7 +280,7 @@ subroutine source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalph
real(pReal), intent(out) :: &
localphiDot, &
dLocalphiDot_dPhi
integer(pInt) :: &
integer :: &
sourceOffset
sourceOffset = source_damage_anisoBrittle_offset(phase)
@ -303,27 +299,27 @@ function source_damage_anisoBrittle_postResults(phase, constituent)
use material, only: &
sourceState
integer(pInt), intent(in) :: &
integer, intent(in) :: &
phase, &
constituent
real(pReal), dimension(sum(source_damage_anisoBrittle_sizePostResult(:, &
source_damage_anisoBrittle_instance(phase)))) :: &
source_damage_anisoBrittle_postResults
integer(pInt) :: &
integer :: &
instance, sourceOffset, o, c
instance = source_damage_anisoBrittle_instance(phase)
sourceOffset = source_damage_anisoBrittle_offset(phase)
c = 0_pInt
c = 0
do o = 1_pInt,size(param(instance)%outputID)
do o = 1,size(param(instance)%outputID)
select case(param(instance)%outputID(o))
case (damage_drivingforce_ID)
source_damage_anisoBrittle_postResults(c+1_pInt) = &
source_damage_anisoBrittle_postResults(c+1) = &
sourceState(phase)%p(sourceOffset)%state(1,constituent)
c = c + 1_pInt
c = c + 1
end select
enddo

54
src/source_damage_anisoDuctile.f90
View File

@ -5,17 +5,15 @@
!> @details to be done
!--------------------------------------------------------------------------------------------------
module source_damage_anisoDuctile
use prec, only: &
pReal, &
pInt
use prec
implicit none
private
integer(pInt), dimension(:), allocatable, public, protected :: &
integer, dimension(:), allocatable, public, protected :: &
source_damage_anisoDuctile_offset, & !< which source is my current damage mechanism?
source_damage_anisoDuctile_instance !< instance of damage source mechanism
integer(pInt), dimension(:,:), allocatable, target, public :: &
integer, dimension(:,:), allocatable, target, public :: &
source_damage_anisoDuctile_sizePostResult !< size of each post result output
character(len=64), dimension(:,:), allocatable, target, public :: &
@ -59,8 +57,6 @@ contains
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
subroutine source_damage_anisoDuctile_init
use prec, only: &
pStringLen
use debug, only: &
debug_level,&
debug_constitutive,&
@ -82,10 +78,10 @@ subroutine source_damage_anisoDuctile_init
config_phase
integer(pInt) :: Ninstance,phase,instance,source,sourceOffset
integer(pInt) :: NofMyPhase,p ,i
integer :: Ninstance,phase,instance,source,sourceOffset
integer :: NofMyPhase,p ,i
integer(pInt), dimension(0), parameter :: emptyIntArray = [integer(pInt)::]
integer, dimension(0), parameter :: emptyIntArray = [integer::]
character(len=65536), dimension(0), parameter :: emptyStringArray = [character(len=65536)::]
integer(kind(undefined_ID)) :: &
outputID
@ -98,13 +94,13 @@ subroutine source_damage_anisoDuctile_init
write(6,'(/,a)') ' <<<+- source_'//SOURCE_DAMAGE_ANISODUCTILE_LABEL//' init -+>>>'
Ninstance = count(phase_source == SOURCE_damage_anisoDuctile_ID)
if (Ninstance == 0_pInt) return
if (Ninstance == 0) return
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) &
write(6,'(a16,1x,i5,/)') '# instances:',Ninstance
allocate(source_damage_anisoDuctile_offset(size(config_phase)), source=0_pInt)
allocate(source_damage_anisoDuctile_instance(size(config_phase)), source=0_pInt)
allocate(source_damage_anisoDuctile_offset(size(config_phase)), source=0)
allocate(source_damage_anisoDuctile_instance(size(config_phase)), source=0)
do phase = 1, size(config_phase)
source_damage_anisoDuctile_instance(phase) = count(phase_source(:,1:phase) == source_damage_anisoDuctile_ID)
do source = 1, phase_Nsources(phase)
@ -113,7 +109,7 @@ subroutine source_damage_anisoDuctile_init
enddo
enddo
allocate(source_damage_anisoDuctile_sizePostResult(maxval(phase_Noutput),Ninstance),source=0_pInt)
allocate(source_damage_anisoDuctile_sizePostResult(maxval(phase_Noutput),Ninstance),source=0)
allocate(source_damage_anisoDuctile_output(maxval(phase_Noutput),Ninstance))
source_damage_anisoDuctile_output = ''
@ -146,18 +142,18 @@ subroutine source_damage_anisoDuctile_init
!--------------------------------------------------------------------------------------------------
! exit if any parameter is out of range
if (extmsg /= '') &
call IO_error(211_pInt,ext_msg=trim(extmsg)//'('//SOURCE_DAMAGE_ANISODUCTILE_LABEL//')')
call IO_error(211,ext_msg=trim(extmsg)//'('//SOURCE_DAMAGE_ANISODUCTILE_LABEL//')')
!--------------------------------------------------------------------------------------------------
! output pararameters
outputs = config%getStrings('(output)',defaultVal=emptyStringArray)
allocate(prm%outputID(0))
do i=1_pInt, size(outputs)
do i=1, size(outputs)
outputID = undefined_ID
select case(outputs(i))
case ('anisoductile_drivingforce')
source_damage_anisoDuctile_sizePostResult(i,source_damage_anisoDuctile_instance(p)) = 1_pInt
source_damage_anisoDuctile_sizePostResult(i,source_damage_anisoDuctile_instance(p)) = 1
source_damage_anisoDuctile_output(i,source_damage_anisoDuctile_instance(p)) = outputs(i)
prm%outputID = [prm%outputID, damage_drivingforce_ID]
@ -173,7 +169,7 @@ subroutine source_damage_anisoDuctile_init
instance = source_damage_anisoDuctile_instance(phase)
sourceOffset = source_damage_anisoDuctile_offset(phase)
call material_allocateSourceState(phase,sourceOffset,NofMyPhase,1_pInt,1_pInt,0_pInt)
call material_allocateSourceState(phase,sourceOffset,NofMyPhase,1,1,0)
sourceState(phase)%p(sourceOffset)%sizePostResults = sum(source_damage_anisoDuctile_sizePostResult(:,instance))
sourceState(phase)%p(sourceOffset)%aTolState=param(instance)%aTol
@ -193,11 +189,11 @@ subroutine source_damage_anisoDuctile_dotState(ipc, ip, el)
damage, &
damageMapping
integer(pInt), intent(in) :: &
integer, intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
integer(pInt) :: &
integer :: &
phase, &
constituent, &
sourceOffset, &
@ -229,7 +225,7 @@ subroutine source_damage_anisoDuctile_getRateAndItsTangent(localphiDot, dLocalph
use material, only: &
sourceState
integer(pInt), intent(in) :: &
integer, intent(in) :: &
phase, &
constituent
real(pReal), intent(in) :: &
@ -237,7 +233,7 @@ subroutine source_damage_anisoDuctile_getRateAndItsTangent(localphiDot, dLocalph
real(pReal), intent(out) :: &
localphiDot, &
dLocalphiDot_dPhi
integer(pInt) :: &
integer :: &
sourceOffset
sourceOffset = source_damage_anisoDuctile_offset(phase)
@ -256,27 +252,27 @@ function source_damage_anisoDuctile_postResults(phase, constituent)
use material, only: &
sourceState
integer(pInt), intent(in) :: &
integer, intent(in) :: &
phase, &
constituent
real(pReal), dimension(sum(source_damage_anisoDuctile_sizePostResult(:, &
source_damage_anisoDuctile_instance(phase)))) :: &
source_damage_anisoDuctile_postResults
integer(pInt) :: &
integer :: &
instance, sourceOffset, o, c
instance = source_damage_anisoDuctile_instance(phase)
sourceOffset = source_damage_anisoDuctile_offset(phase)
c = 0_pInt
c = 0
do o = 1_pInt,size(param(instance)%outputID)
do o = 1,size(param(instance)%outputID)
select case(param(instance)%outputID(o))
case (damage_drivingforce_ID)
source_damage_anisoDuctile_postResults(c+1_pInt) = &
source_damage_anisoDuctile_postResults(c+1) = &
sourceState(phase)%p(sourceOffset)%state(1,constituent)
c = c + 1_pInt
c = c + 1
end select
enddo

52
src/source_damage_isoBrittle.f90
View File

@ -5,17 +5,15 @@
!> @details to be done
!--------------------------------------------------------------------------------------------------
module source_damage_isoBrittle
use prec, only: &
pReal, &
pInt
use prec
implicit none
private
integer(pInt), dimension(:), allocatable, public, protected :: &
integer, dimension(:), allocatable, public, protected :: &
source_damage_isoBrittle_offset, & !< which source is my current damage mechanism?
source_damage_isoBrittle_instance !< instance of damage source mechanism
integer(pInt), dimension(:,:), allocatable, target, public :: &
integer, dimension(:,:), allocatable, target, public :: &
source_damage_isoBrittle_sizePostResult !< size of each post result output
character(len=64), dimension(:,:), allocatable, target, public :: &
@ -53,8 +51,6 @@ contains
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
subroutine source_damage_isoBrittle_init
use prec, only: &
pStringLen
use debug, only: &
debug_level,&
debug_constitutive,&
@ -75,8 +71,8 @@ subroutine source_damage_isoBrittle_init
material_Nphase
integer(pInt) :: Ninstance,phase,instance,source,sourceOffset
integer(pInt) :: NofMyPhase,p,i
integer :: Ninstance,phase,instance,source,sourceOffset
integer :: NofMyPhase,p,i
character(len=65536), dimension(0), parameter :: emptyStringArray = [character(len=65536)::]
integer(kind(undefined_ID)) :: &
outputID
@ -88,14 +84,14 @@ subroutine source_damage_isoBrittle_init
write(6,'(/,a)') ' <<<+- source_'//SOURCE_DAMAGE_ISOBRITTLE_LABEL//' init -+>>>'
Ninstance = int(count(phase_source == SOURCE_damage_isoBrittle_ID),pInt)
if (Ninstance == 0_pInt) return
Ninstance = int(count(phase_source == SOURCE_damage_isoBrittle_ID))
if (Ninstance == 0) return
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) &
write(6,'(a16,1x,i5,/)') '# instances:',Ninstance
allocate(source_damage_isoBrittle_offset(material_Nphase), source=0_pInt)
allocate(source_damage_isoBrittle_instance(material_Nphase), source=0_pInt)
allocate(source_damage_isoBrittle_offset(material_Nphase), source=0)
allocate(source_damage_isoBrittle_instance(material_Nphase), source=0)
do phase = 1, material_Nphase
source_damage_isoBrittle_instance(phase) = count(phase_source(:,1:phase) == source_damage_isoBrittle_ID)
do source = 1, phase_Nsources(phase)
@ -104,7 +100,7 @@ subroutine source_damage_isoBrittle_init
enddo
enddo
allocate(source_damage_isoBrittle_sizePostResult(maxval(phase_Noutput),Ninstance),source=0_pInt)
allocate(source_damage_isoBrittle_sizePostResult(maxval(phase_Noutput),Ninstance),source=0)
allocate(source_damage_isoBrittle_output(maxval(phase_Noutput),Ninstance))
source_damage_isoBrittle_output = ''
@ -129,18 +125,18 @@ subroutine source_damage_isoBrittle_init
!--------------------------------------------------------------------------------------------------
! exit if any parameter is out of range
if (extmsg /= '') &
call IO_error(211_pInt,ext_msg=trim(extmsg)//'('//SOURCE_DAMAGE_ISOBRITTLE_LABEL//')')
call IO_error(211,ext_msg=trim(extmsg)//'('//SOURCE_DAMAGE_ISOBRITTLE_LABEL//')')
!--------------------------------------------------------------------------------------------------
! output pararameters
outputs = config%getStrings('(output)',defaultVal=emptyStringArray)
allocate(prm%outputID(0))
do i=1_pInt, size(outputs)
do i=1, size(outputs)
outputID = undefined_ID
select case(outputs(i))
case ('isobrittle_drivingforce')
source_damage_isoBrittle_sizePostResult(i,source_damage_isoBrittle_instance(p)) = 1_pInt
source_damage_isoBrittle_sizePostResult(i,source_damage_isoBrittle_instance(p)) = 1
source_damage_isoBrittle_output(i,source_damage_isoBrittle_instance(p)) = outputs(i)
prm%outputID = [prm%outputID, damage_drivingforce_ID]
@ -156,7 +152,7 @@ subroutine source_damage_isoBrittle_init
instance = source_damage_isoBrittle_instance(phase)
sourceOffset = source_damage_isoBrittle_offset(phase)
call material_allocateSourceState(phase,sourceOffset,NofMyPhase,1_pInt,1_pInt,1_pInt)
call material_allocateSourceState(phase,sourceOffset,NofMyPhase,1,1,1)
sourceState(phase)%p(sourceOffset)%sizePostResults = sum(source_damage_isoBrittle_sizePostResult(:,instance))
sourceState(phase)%p(sourceOffset)%aTolState=param(instance)%aTol
@ -175,7 +171,7 @@ subroutine source_damage_isoBrittle_deltaState(C, Fe, ipc, ip, el)
math_sym33to6, &
math_I3
integer(pInt), intent(in) :: &
integer, intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
@ -183,7 +179,7 @@ subroutine source_damage_isoBrittle_deltaState(C, Fe, ipc, ip, el)
Fe
real(pReal), intent(in), dimension(6,6) :: &
C
integer(pInt) :: &
integer :: &
phase, constituent, instance, sourceOffset
real(pReal) :: &
strain(6), &
@ -219,7 +215,7 @@ subroutine source_damage_isoBrittle_getRateAndItsTangent(localphiDot, dLocalphiD
use material, only: &
sourceState
integer(pInt), intent(in) :: &
integer, intent(in) :: &
phase, &
constituent
real(pReal), intent(in) :: &
@ -227,7 +223,7 @@ subroutine source_damage_isoBrittle_getRateAndItsTangent(localphiDot, dLocalphiD
real(pReal), intent(out) :: &
localphiDot, &
dLocalphiDot_dPhi
integer(pInt) :: &
integer :: &
instance, sourceOffset
instance = source_damage_isoBrittle_instance(phase)
@ -248,25 +244,25 @@ function source_damage_isoBrittle_postResults(phase, constituent)
use material, only: &
sourceState
integer(pInt), intent(in) :: &
integer, intent(in) :: &
phase, &
constituent
real(pReal), dimension(sum(source_damage_isoBrittle_sizePostResult(:, &
source_damage_isoBrittle_instance(phase)))) :: &
source_damage_isoBrittle_postResults
integer(pInt) :: &
integer :: &
instance, sourceOffset, o, c
instance = source_damage_isoBrittle_instance(phase)
sourceOffset = source_damage_isoBrittle_offset(phase)
c = 0_pInt
c = 0
do o = 1_pInt,size(param(instance)%outputID)
do o = 1,size(param(instance)%outputID)
select case(param(instance)%outputID(o))
case (damage_drivingforce_ID)
source_damage_isoBrittle_postResults(c+1_pInt) = sourceState(phase)%p(sourceOffset)%state(1,constituent)
source_damage_isoBrittle_postResults(c+1) = sourceState(phase)%p(sourceOffset)%state(1,constituent)
c = c + 1
end select

48
src/source_damage_isoDuctile.f90
View File

@ -5,17 +5,15 @@
!> @details to be done
!--------------------------------------------------------------------------------------------------
module source_damage_isoDuctile
use prec, only: &
pReal, &
pInt
use prec
implicit none
private
integer(pInt), dimension(:), allocatable, public, protected :: &
integer, dimension(:), allocatable, public, protected :: &
source_damage_isoDuctile_offset, & !< which source is my current damage mechanism?
source_damage_isoDuctile_instance !< instance of damage source mechanism
integer(pInt), dimension(:,:), allocatable, target, public :: &
integer, dimension(:,:), allocatable, target, public :: &
source_damage_isoDuctile_sizePostResult !< size of each post result output
character(len=64), dimension(:,:), allocatable, target, public :: &
@ -53,8 +51,6 @@ contains
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
subroutine source_damage_isoDuctile_init
use prec, only: &
pStringLen
use debug, only: &
debug_level,&
debug_constitutive,&
@ -75,8 +71,8 @@ subroutine source_damage_isoDuctile_init
material_Nphase
integer(pInt) :: Ninstance,phase,instance,source,sourceOffset
integer(pInt) :: NofMyPhase,p,i
integer :: Ninstance,phase,instance,source,sourceOffset
integer :: NofMyPhase,p,i
character(len=65536), dimension(0), parameter :: emptyStringArray = [character(len=65536)::]
integer(kind(undefined_ID)) :: &
outputID
@ -89,13 +85,13 @@ subroutine source_damage_isoDuctile_init
write(6,'(/,a)') ' <<<+- source_'//SOURCE_DAMAGE_ISODUCTILE_LABEL//' init -+>>>'
Ninstance = count(phase_source == SOURCE_damage_isoDuctile_ID)
if (Ninstance == 0_pInt) return
if (Ninstance == 0) return
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) &
write(6,'(a16,1x,i5,/)') '# instances:',Ninstance
allocate(source_damage_isoDuctile_offset(material_Nphase), source=0_pInt)
allocate(source_damage_isoDuctile_instance(material_Nphase), source=0_pInt)
allocate(source_damage_isoDuctile_offset(material_Nphase), source=0)
allocate(source_damage_isoDuctile_instance(material_Nphase), source=0)
do phase = 1, material_Nphase
source_damage_isoDuctile_instance(phase) = count(phase_source(:,1:phase) == source_damage_isoDuctile_ID)
do source = 1, phase_Nsources(phase)
@ -104,7 +100,7 @@ subroutine source_damage_isoDuctile_init
enddo
enddo
allocate(source_damage_isoDuctile_sizePostResult(maxval(phase_Noutput),Ninstance),source=0_pInt)
allocate(source_damage_isoDuctile_sizePostResult(maxval(phase_Noutput),Ninstance),source=0)
allocate(source_damage_isoDuctile_output(maxval(phase_Noutput),Ninstance))
source_damage_isoDuctile_output = ''
@ -129,18 +125,18 @@ subroutine source_damage_isoDuctile_init
!--------------------------------------------------------------------------------------------------
! exit if any parameter is out of range
if (extmsg /= '') &
call IO_error(211_pInt,ext_msg=trim(extmsg)//'('//SOURCE_DAMAGE_ISODUCTILE_LABEL//')')
call IO_error(211,ext_msg=trim(extmsg)//'('//SOURCE_DAMAGE_ISODUCTILE_LABEL//')')
!--------------------------------------------------------------------------------------------------
! output pararameters
outputs = config%getStrings('(output)',defaultVal=emptyStringArray)
allocate(prm%outputID(0))
do i=1_pInt, size(outputs)
do i=1, size(outputs)
outputID = undefined_ID
select case(outputs(i))
case ('isoductile_drivingforce')
source_damage_isoDuctile_sizePostResult(i,source_damage_isoDuctile_instance(p)) = 1_pInt
source_damage_isoDuctile_sizePostResult(i,source_damage_isoDuctile_instance(p)) = 1
source_damage_isoDuctile_output(i,source_damage_isoDuctile_instance(p)) = outputs(i)
prm%outputID = [prm%outputID, damage_drivingforce_ID]
@ -155,7 +151,7 @@ subroutine source_damage_isoDuctile_init
instance = source_damage_isoDuctile_instance(phase)
sourceOffset = source_damage_isoDuctile_offset(phase)
call material_allocateSourceState(phase,sourceOffset,NofMyPhase,1_pInt,1_pInt,0_pInt)
call material_allocateSourceState(phase,sourceOffset,NofMyPhase,1,1,0)
sourceState(phase)%p(sourceOffset)%sizePostResults = sum(source_damage_isoDuctile_sizePostResult(:,instance))
sourceState(phase)%p(sourceOffset)%aTolState=param(instance)%aTol
@ -176,11 +172,11 @@ subroutine source_damage_isoDuctile_dotState(ipc, ip, el)
damage, &
damageMapping
integer(pInt), intent(in) :: &
integer, intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
integer(pInt) :: &
integer :: &
phase, constituent, instance, homog, sourceOffset, damageOffset
phase = phaseAt(ipc,ip,el)
@ -204,7 +200,7 @@ subroutine source_damage_isoDuctile_getRateAndItsTangent(localphiDot, dLocalphiD
use material, only: &
sourceState
integer(pInt), intent(in) :: &
integer, intent(in) :: &
phase, &
constituent
real(pReal), intent(in) :: &
@ -212,7 +208,7 @@ subroutine source_damage_isoDuctile_getRateAndItsTangent(localphiDot, dLocalphiD
real(pReal), intent(out) :: &
localphiDot, &
dLocalphiDot_dPhi
integer(pInt) :: &
integer :: &
sourceOffset
sourceOffset = source_damage_isoDuctile_offset(phase)
@ -231,25 +227,25 @@ function source_damage_isoDuctile_postResults(phase, constituent)
use material, only: &
sourceState
integer(pInt), intent(in) :: &
integer, intent(in) :: &
phase, &
constituent
real(pReal), dimension(sum(source_damage_isoDuctile_sizePostResult(:, &
source_damage_isoDuctile_instance(phase)))) :: &
source_damage_isoDuctile_postResults
integer(pInt) :: &
integer :: &
instance, sourceOffset, o, c
instance = source_damage_isoDuctile_instance(phase)
sourceOffset = source_damage_isoDuctile_offset(phase)
c = 0_pInt
c = 0
do o = 1_pInt,size(param(instance)%outputID)
do o = 1,size(param(instance)%outputID)
select case(param(instance)%outputID(o))
case (damage_drivingforce_ID)
source_damage_isoDuctile_postResults(c+1_pInt) = sourceState(phase)%p(sourceOffset)%state(1,constituent)
source_damage_isoDuctile_postResults(c+1) = sourceState(phase)%p(sourceOffset)%state(1,constituent)
c = c + 1
end select