Merge branch 'development' into misc-improvements
This commit is contained in:
Martin Diehl
commit
cf1cfb0cd9
138
src/CPFEM.f90
138
src/CPFEM.f90
|
|
@ -35,29 +35,15 @@ module CPFEM
|
|||
CPFEM_dcsdE !< Cauchy stress tangent
|
||||
real(pReal), dimension (:,:,:,:), allocatable, private :: &
|
||||
CPFEM_dcsdE_knownGood !< known good tangent
|
||||
integer(pInt), public :: &
|
||||
cycleCounter = 0_pInt, & !< needs description
|
||||
theInc = -1_pInt, & !< needs description
|
||||
lastLovl = 0_pInt !< lovl in previous call to marc hypela2
|
||||
real(pReal), public :: &
|
||||
theTime = 0.0_pReal, & !< needs description
|
||||
theDelta = 0.0_pReal
|
||||
logical, public :: &
|
||||
outdatedFFN1 = .false., & !< needs description
|
||||
lastIncConverged = .false., & !< needs description
|
||||
outdatedByNewInc = .false. !< needs description
|
||||
|
||||
logical, public, protected :: &
|
||||
CPFEM_init_done = .false. !< remember whether init has been done already
|
||||
logical, private :: &
|
||||
CPFEM_calc_done = .false. !< remember whether first ip has already calced the results
|
||||
integer(pInt), public :: &
|
||||
cycleCounter = 0_pInt !< needs description
|
||||
|
||||
integer(pInt), parameter, public :: &
|
||||
CPFEM_COLLECT = 2_pInt**0_pInt, &
|
||||
CPFEM_CALCRESULTS = 2_pInt**1_pInt, &
|
||||
CPFEM_AGERESULTS = 2_pInt**2_pInt, &
|
||||
CPFEM_BACKUPJACOBIAN = 2_pInt**3_pInt, &
|
||||
CPFEM_RESTOREJACOBIAN = 2_pInt**4_pInt
|
||||
CPFEM_CALCRESULTS = 2_pInt**0_pInt, &
|
||||
CPFEM_AGERESULTS = 2_pInt**1_pInt, &
|
||||
CPFEM_BACKUPJACOBIAN = 2_pInt**2_pInt, &
|
||||
CPFEM_RESTOREJACOBIAN = 2_pInt**3_pInt
|
||||
|
||||
public :: &
|
||||
CPFEM_general, &
|
||||
|
|
@ -68,14 +54,10 @@ contains
|
|||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief call (thread safe) all module initializations
|
||||
!> @brief call all module initializations
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
subroutine CPFEM_initAll(el,ip)
|
||||
subroutine CPFEM_initAll
|
||||
|
||||
integer(pInt), intent(in) :: el, & !< FE el number
|
||||
ip !< FE integration point number
|
||||
|
||||
CPFEM_init_done = .true.
|
||||
call DAMASK_interface_init
|
||||
call prec_init
|
||||
call IO_init
|
||||
|
|
@ -88,7 +70,7 @@ subroutine CPFEM_initAll(el,ip)
|
|||
call YAML_init
|
||||
call HDF5_utilities_init
|
||||
call results_init(.false.)
|
||||
call discretization_marc_init(ip, el)
|
||||
call discretization_marc_init
|
||||
call lattice_init
|
||||
call material_init(.false.)
|
||||
call constitutive_init
|
||||
|
|
@ -124,7 +106,7 @@ end subroutine CPFEM_init
|
|||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief perform initialization at first call, update variables and call the actual material model
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyStress, jacobian)
|
||||
subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyStress, jacobian)
|
||||
|
||||
integer(pInt), intent(in) :: elFE, & !< FE element number
|
||||
ip !< integration point number
|
||||
|
|
@ -133,24 +115,21 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt
|
|||
ffn1 !< deformation gradient for t=t1
|
||||
integer(pInt), intent(in) :: mode !< computation mode 1: regular computation plus aging of results
|
||||
real(pReal), intent(in) :: temperature_inp !< temperature
|
||||
logical, intent(in) :: parallelExecution !< flag indicating parallel computation of requested IPs
|
||||
real(pReal), dimension(6), intent(out) :: cauchyStress !< stress as 6 vector
|
||||
real(pReal), dimension(6,6), intent(out) :: jacobian !< jacobian as 66 tensor (Consistent tangent dcs/dE)
|
||||
|
||||
real(pReal) J_inverse, & ! inverse of Jacobian
|
||||
rnd
|
||||
real(pReal), dimension (3,3) :: Kirchhoff, & ! Piola-Kirchhoff stress
|
||||
cauchyStress33 ! stress vector
|
||||
real(pReal), dimension (3,3) :: Kirchhoff ! Piola-Kirchhoff stress
|
||||
real(pReal), dimension (3,3,3,3) :: H_sym, &
|
||||
H, &
|
||||
jacobian3333 ! jacobian in Matrix notation
|
||||
H
|
||||
|
||||
integer(pInt) elCP, & ! crystal plasticity element number
|
||||
i, j, k, l, m, n, ph, homog, mySource
|
||||
logical updateJaco ! flag indicating if Jacobian has to be updated
|
||||
|
||||
real(pReal), parameter :: ODD_STRESS = 1e15_pReal, & !< return value for stress in case of ping pong dummy cycle
|
||||
ODD_JACOBIAN = 1e50_pReal !< return value for jacobian in case of ping pong dummy cycle
|
||||
real(pReal), parameter :: ODD_STRESS = 1e15_pReal, & !< return value for stress if terminallyIll
|
||||
ODD_JACOBIAN = 1e50_pReal !< return value for jacobian if terminallyIll
|
||||
|
||||
elCP = mesh_FEM2DAMASK_elem(elFE)
|
||||
|
||||
|
|
@ -159,9 +138,6 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt
|
|||
write(6,'(/,a)') '#############################################'
|
||||
write(6,'(a1,a22,1x,i8,a13)') '#','element', elCP, '#'
|
||||
write(6,'(a1,a22,1x,i8,a13)') '#','ip', ip, '#'
|
||||
write(6,'(a1,a22,1x,f15.7,a6)') '#','theTime', theTime, '#'
|
||||
write(6,'(a1,a22,1x,f15.7,a6)') '#','theDelta', theDelta, '#'
|
||||
write(6,'(a1,a22,1x,i8,a13)') '#','theInc', theInc, '#'
|
||||
write(6,'(a1,a22,1x,i8,a13)') '#','cycleCounter', cycleCounter, '#'
|
||||
write(6,'(a1,a22,1x,i8,a13)') '#','computationMode',mode, '#'
|
||||
if (terminallyIll) &
|
||||
|
|
@ -174,13 +150,8 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt
|
|||
if (iand(mode, CPFEM_RESTOREJACOBIAN) /= 0_pInt) &
|
||||
CPFEM_dcsde = CPFEM_dcsde_knownGood
|
||||
|
||||
!*** age results
|
||||
if (iand(mode, CPFEM_AGERESULTS) /= 0_pInt) call CPFEM_forward
|
||||
|
||||
|
||||
!*** collection of FEM input with returning of randomize odd stress and jacobian
|
||||
!* If no parallel execution is required, there is no need to collect FEM input
|
||||
if (.not. parallelExecution) then
|
||||
chosenThermal1: select case (thermal_type(material_homogenizationAt(elCP)))
|
||||
case (THERMAL_conduction_ID) chosenThermal1
|
||||
temperature(material_homogenizationAt(elCP))%p(thermalMapping(material_homogenizationAt(elCP))%p(ip,elCP)) = &
|
||||
|
|
@ -189,64 +160,22 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt
|
|||
materialpoint_F0(1:3,1:3,ip,elCP) = ffn
|
||||
materialpoint_F(1:3,1:3,ip,elCP) = ffn1
|
||||
|
||||
elseif (iand(mode, CPFEM_COLLECT) /= 0_pInt) then
|
||||
call random_number(rnd)
|
||||
if (rnd < 0.5_pReal) rnd = rnd - 1.0_pReal
|
||||
CPFEM_cs(1:6,ip,elCP) = rnd * ODD_STRESS
|
||||
CPFEM_dcsde(1:6,1:6,ip,elCP) = ODD_JACOBIAN * math_identity2nd(6)
|
||||
chosenThermal2: select case (thermal_type(material_homogenizationAt(elCP)))
|
||||
case (THERMAL_conduction_ID) chosenThermal2
|
||||
temperature(material_homogenizationAt(elCP))%p(thermalMapping(material_homogenizationAt(elCP))%p(ip,elCP)) = &
|
||||
temperature_inp
|
||||
end select chosenThermal2
|
||||
materialpoint_F0(1:3,1:3,ip,elCP) = ffn
|
||||
materialpoint_F(1:3,1:3,ip,elCP) = ffn1
|
||||
CPFEM_calc_done = .false.
|
||||
endif
|
||||
|
||||
|
||||
!*** calculation of stress and jacobian
|
||||
if (iand(mode, CPFEM_CALCRESULTS) /= 0_pInt) then
|
||||
|
||||
!*** deformation gradient outdated or any actual deformation gradient differs more than relevantStrain from the stored one
|
||||
validCalculation: if (terminallyIll &
|
||||
.or. outdatedFFN1 &
|
||||
.or. any(abs(ffn1 - materialpoint_F(1:3,1:3,ip,elCP)) > defgradTolerance)) then
|
||||
if (any(abs(ffn1 - materialpoint_F(1:3,1:3,ip,elCP)) > defgradTolerance)) then
|
||||
if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0_pInt) then
|
||||
write(6,'(a,1x,i8,1x,i2)') '<< CPFEM >> OUTDATED at elFE ip',elFE,ip
|
||||
write(6,'(a,/,3(12x,3(f10.6,1x),/))') '<< CPFEM >> FFN1 old:',&
|
||||
transpose(materialpoint_F(1:3,1:3,ip,elCP))
|
||||
write(6,'(a,/,3(12x,3(f10.6,1x),/))') '<< CPFEM >> FFN1 now:',transpose(ffn1)
|
||||
endif
|
||||
outdatedFFN1 = .true.
|
||||
endif
|
||||
validCalculation: if (terminallyIll) then
|
||||
call random_number(rnd)
|
||||
if (rnd < 0.5_pReal) rnd = rnd - 1.0_pReal
|
||||
CPFEM_cs(1:6,ip,elCP) = ODD_STRESS * rnd
|
||||
CPFEM_dcsde(1:6,1:6,ip,elCP) = ODD_JACOBIAN * math_identity2nd(6)
|
||||
|
||||
!*** deformation gradient is not outdated
|
||||
else validCalculation
|
||||
updateJaco = mod(cycleCounter,iJacoStiffness) == 0
|
||||
!* no parallel computation, so we use just one single elFE and ip for computation
|
||||
if (.not. parallelExecution) then
|
||||
FEsolving_execElem = elCP
|
||||
FEsolving_execIP = ip
|
||||
if (iand(debug_level(debug_CPFEM), debug_levelExtensive) /= 0_pInt) &
|
||||
write(6,'(a,i8,1x,i2)') '<< CPFEM >> calculation for elFE ip ',elFE,ip
|
||||
call materialpoint_stressAndItsTangent(updateJaco, dt)
|
||||
FEsolving_execElem = elCP
|
||||
FEsolving_execIP = ip
|
||||
if (iand(debug_level(debug_CPFEM), debug_levelExtensive) /= 0_pInt) &
|
||||
write(6,'(a,i8,1x,i2)') '<< CPFEM >> calculation for elFE ip ',elFE,ip
|
||||
call materialpoint_stressAndItsTangent(updateJaco, dt)
|
||||
|
||||
!* parallel computation and calulation not yet done
|
||||
elseif (.not. CPFEM_calc_done) then
|
||||
if (iand(debug_level(debug_CPFEM), debug_levelExtensive) /= 0_pInt) &
|
||||
write(6,'(a,i8,a,i8)') '<< CPFEM >> calculation for elements ',FEsolving_execElem(1),&
|
||||
' to ',FEsolving_execElem(2)
|
||||
call materialpoint_stressAndItsTangent(updateJaco, dt)
|
||||
CPFEM_calc_done = .true.
|
||||
endif
|
||||
|
||||
!* map stress and stiffness (or return odd values if terminally ill)
|
||||
terminalIllness: if (terminallyIll) then
|
||||
|
||||
call random_number(rnd)
|
||||
|
|
@ -256,7 +185,7 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt
|
|||
|
||||
else terminalIllness
|
||||
|
||||
! translate from P to CS
|
||||
! translate from P to sigma
|
||||
Kirchhoff = matmul(materialpoint_P(1:3,1:3,ip,elCP), transpose(materialpoint_F(1:3,1:3,ip,elCP)))
|
||||
J_inverse = 1.0_pReal / math_det33(materialpoint_F(1:3,1:3,ip,elCP))
|
||||
CPFEM_cs(1:6,ip,elCP) = math_sym33to6(J_inverse * Kirchhoff,weighted=.false.)
|
||||
|
|
@ -280,7 +209,6 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt
|
|||
endif terminalIllness
|
||||
endif validCalculation
|
||||
|
||||
!* report stress and stiffness
|
||||
if ((iand(debug_level(debug_CPFEM), debug_levelExtensive) /= 0_pInt) &
|
||||
.and. ((debug_e == elCP .and. debug_i == ip) &
|
||||
.or. .not. iand(debug_level(debug_CPFEM), debug_levelSelective) /= 0_pInt)) then
|
||||
|
|
@ -293,35 +221,11 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt
|
|||
|
||||
endif
|
||||
|
||||
!*** warn if stiffness close to zero
|
||||
if (all(abs(CPFEM_dcsdE(1:6,1:6,ip,elCP)) < 1e-10_pReal)) call IO_warning(601,elCP,ip)
|
||||
|
||||
!*** copy to output if using commercial FEM solver
|
||||
cauchyStress = CPFEM_cs (1:6, ip,elCP)
|
||||
jacobian = CPFEM_dcsdE(1:6,1:6,ip,elCP)
|
||||
|
||||
|
||||
!*** remember extreme values of stress ...
|
||||
cauchyStress33 = math_6toSym33(CPFEM_cs(1:6,ip,elCP),weighted=.false.)
|
||||
if (maxval(cauchyStress33) > debug_stressMax) then
|
||||
debug_stressMaxLocation = [elCP, ip]
|
||||
debug_stressMax = maxval(cauchyStress33)
|
||||
endif
|
||||
if (minval(cauchyStress33) < debug_stressMin) then
|
||||
debug_stressMinLocation = [elCP, ip]
|
||||
debug_stressMin = minval(cauchyStress33)
|
||||
endif
|
||||
!*** ... and Jacobian
|
||||
jacobian3333 = math_66toSym3333(CPFEM_dcsdE(1:6,1:6,ip,elCP),weighted=.false.)
|
||||
if (maxval(jacobian3333) > debug_jacobianMax) then
|
||||
debug_jacobianMaxLocation = [elCP, ip]
|
||||
debug_jacobianMax = maxval(jacobian3333)
|
||||
endif
|
||||
if (minval(jacobian3333) < debug_jacobianMin) then
|
||||
debug_jacobianMinLocation = [elCP, ip]
|
||||
debug_jacobianMin = minval(jacobian3333)
|
||||
endif
|
||||
|
||||
end subroutine CPFEM_general
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -43,8 +43,6 @@ module DAMASK_interface
|
|||
logical, protected, public :: symmetricSolver
|
||||
character(len=*), parameter, public :: INPUTFILEEXTENSION = '.dat'
|
||||
|
||||
logical, dimension(:,:), public, allocatable :: &
|
||||
calcMode !< calculate or collect (ping pong scheme)
|
||||
|
||||
public :: &
|
||||
DAMASK_interface_init, &
|
||||
|
|
@ -175,62 +173,62 @@ 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
|
||||
use DAMASK_interface
|
||||
use numerics
|
||||
use FEsolving
|
||||
use debug
|
||||
use discretization_marc
|
||||
use homogenization
|
||||
use CPFEM
|
||||
use prec
|
||||
use DAMASK_interface
|
||||
use numerics
|
||||
use FEsolving
|
||||
use debug
|
||||
use discretization_marc
|
||||
use homogenization
|
||||
use CPFEM
|
||||
|
||||
implicit none
|
||||
!$ include "omp_lib.h" ! the openMP function library
|
||||
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
|
||||
nshear, & !< number of shear components
|
||||
ncrd, & !< number of coordinates
|
||||
itel, & !< dimension of F and R, either 2 or 3
|
||||
ndeg, & !< number of degrees of freedom
|
||||
ndm, & !< not specified in MSC.Marc 2012 Manual D
|
||||
nnode, & !< number of nodes per element
|
||||
jtype, & !< element type
|
||||
ifr, & !< set to 1 if R has been calculated
|
||||
ifu !< set to 1 if stretch has been calculated
|
||||
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
|
||||
lclass !< (1) element class, (2) 0: displacement, 1: low order Herrmann, 2: high order Herrmann
|
||||
real(pReal), dimension(*), intent(in) :: & ! has dimension(1) according to MSC.Marc 2012 Manual D, but according to example hypela2.f dimension(*)
|
||||
e, & !< total elastic strain
|
||||
de, & !< increment of strain
|
||||
dt !< increment of state variables
|
||||
real(pReal), dimension(itel), intent(in) :: & ! according to MSC.Marc 2012 Manual D
|
||||
strechn, & !< square of principal stretch ratios, lambda(i) at t=n
|
||||
strechn1 !< square of principal stretch ratios, lambda(i) at t=n+1
|
||||
real(pReal), dimension(3,3), intent(in) :: & ! has dimension(itel,*) according to MSC.Marc 2012 Manual D, but we alway assume dimension(3,3)
|
||||
ffn, & !< deformation gradient at t=n
|
||||
ffn1 !< deformation gradient at t=n+1
|
||||
real(pReal), dimension(itel,*), intent(in) :: & ! according to MSC.Marc 2012 Manual D
|
||||
frotn, & !< rotation tensor at t=n
|
||||
eigvn, & !< i principal direction components for j eigenvalues at t=n
|
||||
frotn1, & !< rotation tensor at t=n+1
|
||||
eigvn1 !< i principal direction components for j eigenvalues at t=n+1
|
||||
real(pReal), dimension(ndeg,*), intent(in) :: & ! according to MSC.Marc 2012 Manual D
|
||||
disp, & !< incremental displacements
|
||||
dispt !< displacements at t=n (at assembly, lovl=4) and displacements at t=n+1 (at stress recovery, lovl=6)
|
||||
real(pReal), dimension(ncrd,*), intent(in) :: & ! according to MSC.Marc 2012 Manual D
|
||||
coord !< coordinates
|
||||
real(pReal), dimension(*), intent(inout) :: & ! according to MSC.Marc 2012 Manual D
|
||||
t !< state variables (comes in at t=n, must be updated to have state variables at t=n+1)
|
||||
real(pReal), dimension(ndi+nshear), intent(out) :: & ! has dimension(*) according to MSC.Marc 2012 Manual D, but we need to loop over it
|
||||
s, & !< stress - should be updated by user
|
||||
g !< change in stress due to temperature effects
|
||||
real(pReal), dimension(ngens,ngens), intent(out) :: & ! according to MSC.Marc 2012 Manual D, but according to example hypela2.f dimension(ngens,*)
|
||||
d !< stress-strain law to be formed
|
||||
implicit none
|
||||
include "omp_lib.h" ! the openMP function library
|
||||
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
|
||||
nshear, & !< number of shear components
|
||||
ncrd, & !< number of coordinates
|
||||
itel, & !< dimension of F and R, either 2 or 3
|
||||
ndeg, & !< number of degrees of freedom
|
||||
ndm, & !< not specified in MSC.Marc 2012 Manual D
|
||||
nnode, & !< number of nodes per element
|
||||
jtype, & !< element type
|
||||
ifr, & !< set to 1 if R has been calculated
|
||||
ifu !< set to 1 if stretch has been calculated
|
||||
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
|
||||
lclass !< (1) element class, (2) 0: displacement, 1: low order Herrmann, 2: high order Herrmann
|
||||
real(pReal), dimension(*), intent(in) :: & ! has dimension(1) according to MSC.Marc 2012 Manual D, but according to example hypela2.f dimension(*)
|
||||
e, & !< total elastic strain
|
||||
de, & !< increment of strain
|
||||
dt !< increment of state variables
|
||||
real(pReal), dimension(itel), intent(in) :: & ! according to MSC.Marc 2012 Manual D
|
||||
strechn, & !< square of principal stretch ratios, lambda(i) at t=n
|
||||
strechn1 !< square of principal stretch ratios, lambda(i) at t=n+1
|
||||
real(pReal), dimension(3,3), intent(in) :: & ! has dimension(itel,*) according to MSC.Marc 2012 Manual D, but we alway assume dimension(3,3)
|
||||
ffn, & !< deformation gradient at t=n
|
||||
ffn1 !< deformation gradient at t=n+1
|
||||
real(pReal), dimension(itel,*), intent(in) :: & ! according to MSC.Marc 2012 Manual D
|
||||
frotn, & !< rotation tensor at t=n
|
||||
eigvn, & !< i principal direction components for j eigenvalues at t=n
|
||||
frotn1, & !< rotation tensor at t=n+1
|
||||
eigvn1 !< i principal direction components for j eigenvalues at t=n+1
|
||||
real(pReal), dimension(ndeg,*), intent(in) :: & ! according to MSC.Marc 2012 Manual D
|
||||
disp, & !< incremental displacements
|
||||
dispt !< displacements at t=n (at assembly, lovl=4) and displacements at t=n+1 (at stress recovery, lovl=6)
|
||||
real(pReal), dimension(ncrd,*), intent(in) :: & ! according to MSC.Marc 2012 Manual D
|
||||
coord !< coordinates
|
||||
real(pReal), dimension(*), intent(inout) :: & ! according to MSC.Marc 2012 Manual D
|
||||
t !< state variables (comes in at t=n, must be updated to have state variables at t=n+1)
|
||||
real(pReal), dimension(ndi+nshear), intent(out) :: & ! has dimension(*) according to MSC.Marc 2012 Manual D, but we need to loop over it
|
||||
s, & !< stress - should be updated by user
|
||||
g !< change in stress due to temperature effects
|
||||
real(pReal), dimension(ngens,ngens), intent(out) :: & ! according to MSC.Marc 2012 Manual D, but according to example hypela2.f dimension(ngens,*)
|
||||
d !< stress-strain law to be formed
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
! Marc common blocks are in fixed format so they have to be reformated to free format (f90)
|
||||
|
|
@ -239,137 +237,108 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, &
|
|||
#include QUOTE(PASTE(./marc/include/concom,Marc4DAMASK)) ! concom is needed for inc, lovl
|
||||
#include QUOTE(PASTE(./marc/include/creeps,Marc4DAMASK)) ! creeps is needed for timinc (time increment)
|
||||
|
||||
logical :: cutBack
|
||||
real(pReal), dimension(6) :: stress
|
||||
real(pReal), dimension(6,6) :: ddsdde
|
||||
integer :: computationMode, i, cp_en, node, CPnodeID
|
||||
!$ integer(4) :: defaultNumThreadsInt !< default value set by Marc
|
||||
logical :: cutBack
|
||||
real(pReal), dimension(6) :: stress
|
||||
real(pReal), dimension(6,6) :: ddsdde
|
||||
integer :: computationMode, i, cp_en, node, CPnodeID
|
||||
integer(4) :: defaultNumThreadsInt !< default value set by Marc
|
||||
|
||||
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
|
||||
write(6,'(a,i1)') ' Shear stress: ', nshear
|
||||
write(6,'(a,i2)') ' DoF: ', ndeg
|
||||
write(6,'(a,i2)') ' Coordinates: ', ncrd
|
||||
write(6,'(a,i12)') ' Nodes: ', nnode
|
||||
write(6,'(a,i1)') ' Deformation gradient: ', itel
|
||||
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' Deformation gradient at t=n:', &
|
||||
transpose(ffn)
|
||||
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' Deformation gradient at t=n+1:', &
|
||||
transpose(ffn1)
|
||||
endif
|
||||
integer(pInt), save :: &
|
||||
theInc = -1_pInt, & !< needs description
|
||||
lastLovl = 0_pInt !< lovl in previous call to marc hypela2
|
||||
real(pReal), save :: &
|
||||
theTime = 0.0_pReal, & !< needs description
|
||||
theDelta = 0.0_pReal
|
||||
logical, save :: &
|
||||
lastIncConverged = .false., & !< needs description
|
||||
outdatedByNewInc = .false., & !< needs description
|
||||
CPFEM_init_done = .false. !< remember whether init has been done already
|
||||
|
||||
!$ defaultNumThreadsInt = omp_get_num_threads() ! remember number of threads set by Marc
|
||||
!$ call omp_set_num_threads(DAMASK_NumThreadsInt) ! set number of threads for parallel execution set by DAMASK_NUM_THREADS
|
||||
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
|
||||
write(6,'(a,i1)') ' Shear stress: ', nshear
|
||||
write(6,'(a,i2)') ' DoF: ', ndeg
|
||||
write(6,'(a,i2)') ' Coordinates: ', ncrd
|
||||
write(6,'(a,i12)') ' Nodes: ', nnode
|
||||
write(6,'(a,i1)') ' Deformation gradient: ', itel
|
||||
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' Deformation gradient at t=n:', &
|
||||
transpose(ffn)
|
||||
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' Deformation gradient at t=n+1:', &
|
||||
transpose(ffn1)
|
||||
endif
|
||||
|
||||
if (.not. CPFEM_init_done) call CPFEM_initAll(m(1),nn)
|
||||
defaultNumThreadsInt = omp_get_num_threads() ! remember number of threads set by Marc
|
||||
call omp_set_num_threads(1) ! no openMP
|
||||
|
||||
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
|
||||
elseif (lovl == 6) then ! stress requested by marc
|
||||
cp_en = mesh_FEM2DAMASK_elem(m(1))
|
||||
if (cptim > theTime .or. inc /= theInc) then ! reached "convergence"
|
||||
terminallyIll = .false.
|
||||
cycleCounter = -1 ! first calc step increments this to cycle = 0
|
||||
if (inc == 0) then ! >> start of analysis <<
|
||||
lastIncConverged = .false. ! no Jacobian backup
|
||||
outdatedByNewInc = .false. ! no aging of state
|
||||
calcMode = .false. ! pretend last step was collection
|
||||
lastLovl = lovl ! pretend that this is NOT the first after a lovl change
|
||||
write(6,'(a,i6,1x,i2)') '<< HYPELA2 >> start of analysis..! ',m(1),nn
|
||||
flush(6)
|
||||
else if (inc - theInc > 1) then ! >> restart of broken analysis <<
|
||||
lastIncConverged = .false. ! no Jacobian backup
|
||||
outdatedByNewInc = .false. ! no aging of state
|
||||
calcMode = .true. ! pretend last step was calculation
|
||||
write(6,'(a,i6,1x,i2)') '<< HYPELA2 >> restart of analysis..! ',m(1),nn
|
||||
flush(6)
|
||||
else ! >> just the next inc <<
|
||||
lastIncConverged = .true. ! request Jacobian backup
|
||||
outdatedByNewInc = .true. ! request aging of state
|
||||
calcMode = .true. ! assure last step was calculation
|
||||
write(6,'(a,i6,1x,i2)') '<< HYPELA2 >> new increment..! ',m(1),nn
|
||||
flush(6)
|
||||
endif
|
||||
else if ( timinc < theDelta ) then ! >> cutBack <<
|
||||
lastIncConverged = .false. ! no Jacobian backup
|
||||
outdatedByNewInc = .false. ! no aging of state
|
||||
terminallyIll = .false.
|
||||
cycleCounter = -1 ! first calc step increments this to cycle = 0
|
||||
calcMode = .true. ! pretend last step was calculation
|
||||
write(6,'(a,i6,1x,i2)') '<< HYPELA2 >> cutback detected..! ',m(1),nn
|
||||
flush(6)
|
||||
endif ! convergence treatment end
|
||||
if (.not. CPFEM_init_done) then
|
||||
CPFEM_init_done = .true.
|
||||
call CPFEM_initAll
|
||||
endif
|
||||
|
||||
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
|
||||
elseif (lovl == 6) then ! stress requested by marc
|
||||
computationMode = CPFEM_CALCRESULTS
|
||||
cp_en = mesh_FEM2DAMASK_elem(m(1))
|
||||
if (cptim > theTime .or. inc /= theInc) then ! reached "convergence"
|
||||
terminallyIll = .false.
|
||||
cycleCounter = -1 ! first calc step increments this to cycle = 0
|
||||
if (inc == 0) then ! >> start of analysis <<
|
||||
lastIncConverged = .false.
|
||||
outdatedByNewInc = .false.
|
||||
lastLovl = lovl ! pretend that this is NOT the first after a lovl change
|
||||
write(6,'(a,i6,1x,i2)') '<< HYPELA2 >> start of analysis..! ',m(1),nn
|
||||
else if (inc - theInc > 1) then ! >> restart of broken analysis <<
|
||||
lastIncConverged = .false.
|
||||
outdatedByNewInc = .false.
|
||||
write(6,'(a,i6,1x,i2)') '<< HYPELA2 >> restart of analysis..! ',m(1),nn
|
||||
else ! >> just the next inc <<
|
||||
lastIncConverged = .true.
|
||||
outdatedByNewInc = .true.
|
||||
write(6,'(a,i6,1x,i2)') '<< HYPELA2 >> new increment..! ',m(1),nn
|
||||
endif
|
||||
else if ( timinc < theDelta ) then ! >> cutBack <<
|
||||
lastIncConverged = .false.
|
||||
outdatedByNewInc = .false.
|
||||
terminallyIll = .false.
|
||||
cycleCounter = -1 ! first calc step increments this to cycle = 0
|
||||
write(6,'(a,i6,1x,i2)') '<< HYPELA2 >> cutback detected..! ',m(1),nn
|
||||
endif ! convergence treatment end
|
||||
flush(6)
|
||||
|
||||
if (usePingPong) then
|
||||
calcMode(nn,cp_en) = .not. calcMode(nn,cp_en) ! ping pong (calc <--> collect)
|
||||
if (calcMode(nn,cp_en)) then ! now --- CALC ---
|
||||
computationMode = CPFEM_CALCRESULTS
|
||||
if (lastLovl /= lovl) then ! first after ping pong
|
||||
call debug_reset() ! resets debugging
|
||||
outdatedFFN1 = .false.
|
||||
cycleCounter = cycleCounter + 1
|
||||
!mesh_cellnode = mesh_build_cellnodes() ! update cell node coordinates
|
||||
!call mesh_build_ipCoordinates() ! update ip coordinates
|
||||
endif
|
||||
if (outdatedByNewInc) then
|
||||
computationMode = ior(computationMode,CPFEM_AGERESULTS) ! calc and age results
|
||||
outdatedByNewInc = .false. ! reset flag
|
||||
endif
|
||||
else ! now --- COLLECT ---
|
||||
computationMode = CPFEM_COLLECT ! plain collect
|
||||
if (lastLovl /= lovl .and. & .not. terminallyIll) &
|
||||
call debug_info() ! first after ping pong reports (meaningful) debugging
|
||||
if (lastIncConverged) then
|
||||
computationMode = ior(computationMode,CPFEM_BACKUPJACOBIAN) ! collect and backup Jacobian after convergence
|
||||
lastIncConverged = .false. ! reset flag
|
||||
endif
|
||||
!do node = 1,theMesh%elem%nNodes
|
||||
!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
|
||||
if (lastLovl /= lovl) then
|
||||
cycleCounter = cycleCounter + 1
|
||||
!mesh_cellnode = mesh_build_cellnodes() ! update cell node coordinates
|
||||
!call mesh_build_ipCoordinates() ! update ip coordinates
|
||||
endif
|
||||
if (outdatedByNewInc) then
|
||||
computationMode = ior(computationMode,CPFEM_AGERESULTS)
|
||||
outdatedByNewInc = .false.
|
||||
endif
|
||||
if (lastIncConverged) then
|
||||
computationMode = ior(computationMode,CPFEM_BACKUPJACOBIAN)
|
||||
lastIncConverged = .false.
|
||||
endif
|
||||
|
||||
else ! --- PLAIN MODE ---
|
||||
computationMode = CPFEM_CALCRESULTS ! always calc
|
||||
if (lastLovl /= lovl) then
|
||||
if (.not. terminallyIll) &
|
||||
call debug_info() ! first reports (meaningful) debugging
|
||||
call debug_reset() ! and resets debugging
|
||||
outdatedFFN1 = .false.
|
||||
cycleCounter = cycleCounter + 1
|
||||
!mesh_cellnode = mesh_build_cellnodes() ! update cell node coordinates
|
||||
!call mesh_build_ipCoordinates() ! update ip coordinates
|
||||
endif
|
||||
if (outdatedByNewInc) then
|
||||
computationMode = ior(computationMode,CPFEM_AGERESULTS)
|
||||
outdatedByNewInc = .false. ! reset flag
|
||||
endif
|
||||
if (lastIncConverged) then
|
||||
computationMode = ior(computationMode,CPFEM_BACKUPJACOBIAN) ! backup Jacobian after convergence
|
||||
lastIncConverged = .false. ! reset flag
|
||||
endif
|
||||
endif
|
||||
theTime = cptim
|
||||
theDelta = timinc
|
||||
theInc = inc
|
||||
|
||||
theTime = cptim ! record current starting time
|
||||
theDelta = timinc ! record current time increment
|
||||
theInc = inc ! record current increment number
|
||||
endif
|
||||
lastLovl = lovl
|
||||
|
||||
endif
|
||||
lastLovl = lovl ! record lovl
|
||||
call CPFEM_general(computationMode,ffn,ffn1,t(1),timinc,m(1),nn,stress,ddsdde)
|
||||
|
||||
call CPFEM_general(computationMode,usePingPong,ffn,ffn1,t(1),timinc,m(1),nn,stress,ddsdde)
|
||||
d = ddsdde(1:ngens,1:ngens)
|
||||
s = stress(1:ndi+nshear)
|
||||
g = 0.0_pReal
|
||||
if(symmetricSolver) d = 0.5_pReal*(d+transpose(d))
|
||||
|
||||
d = ddsdde(1:ngens,1:ngens)
|
||||
s = stress(1:ndi+nshear)
|
||||
g = 0.0_pReal
|
||||
if(symmetricSolver) d = 0.5_pReal*(d+transpose(d))
|
||||
|
||||
!$ call omp_set_num_threads(defaultNumThreadsInt) ! reset number of threads to stored default value
|
||||
call omp_set_num_threads(defaultNumThreadsInt) ! reset number of threads to stored default value
|
||||
|
||||
end subroutine hypela2
|
||||
|
||||
|
|
|
|||
12
src/IO.f90
12
src/IO.f90
|
|
@ -11,7 +11,6 @@ module IO
|
|||
implicit none
|
||||
private
|
||||
character(len=*), parameter, public :: &
|
||||
IO_EOF = '#EOF#', & !< end of file string
|
||||
IO_WHITESPACE = achar(44)//achar(32)//achar(9)//achar(10)//achar(13) !< whitespace characters
|
||||
character, parameter, public :: &
|
||||
IO_EOL = new_line('DAMASK'), & !< end of line character
|
||||
|
|
@ -106,7 +105,7 @@ function IO_readlines(fileName) result(fileContent)
|
|||
endif
|
||||
startPos = endPos + 2 ! jump to next line start
|
||||
|
||||
fileContent(l) = line
|
||||
fileContent(l) = trim(line)//''
|
||||
l = l + 1
|
||||
enddo
|
||||
|
||||
|
|
@ -114,7 +113,8 @@ end function IO_readlines
|
|||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief reads an entire ASCII file into a string
|
||||
!> @brief read ASCII file into a string
|
||||
!> @details ensures that the string ends with a new line (expected UNIX behavior)
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
function IO_read(fileName) result(fileContent)
|
||||
|
||||
|
|
@ -130,10 +130,14 @@ function IO_read(fileName) result(fileContent)
|
|||
status='old', position='rewind', action='read',iostat=myStat)
|
||||
if(myStat /= 0) call IO_error(100,ext_msg=trim(fileName))
|
||||
allocate(character(len=fileLength)::fileContent)
|
||||
if(fileLength==0) return
|
||||
|
||||
read(fileUnit,iostat=myStat) fileContent
|
||||
if(myStat > 0) call IO_error(102,ext_msg=trim(fileName)) ! <0 for ifort (https://software.intel.com/en-us/comment/1960081)
|
||||
if(myStat /= 0) call IO_error(102,ext_msg=trim(fileName))
|
||||
close(fileUnit)
|
||||
|
||||
if(fileContent(fileLength:fileLength) /= IO_EOL) fileContent = fileContent//IO_EOL ! ensure EOL@EOF
|
||||
|
||||
end function IO_read
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -245,7 +245,7 @@ subroutine crystallite_init
|
|||
enddo
|
||||
!$OMP END PARALLEL DO
|
||||
|
||||
if(any(plasticState%nonlocal) .and. .not. usePingPong) call IO_error(601) ! exit if nonlocal but no ping-pong ToDo: Why not check earlier? or in nonlocal?
|
||||
!if(any(plasticState%nonlocal) .and. .not. usePingPong) call IO_error(601)
|
||||
|
||||
crystallite_partionedFp0 = crystallite_Fp0
|
||||
crystallite_partionedFi0 = crystallite_Fi0
|
||||
|
|
@ -276,9 +276,6 @@ subroutine crystallite_init
|
|||
write(6,'(a42,1x,i10)') 'max # of constituents/integration point: ', cMax
|
||||
flush(6)
|
||||
endif
|
||||
|
||||
call debug_info
|
||||
call debug_reset
|
||||
#endif
|
||||
|
||||
end subroutine crystallite_init
|
||||
|
|
|
|||
|
|
@ -49,26 +49,11 @@ module debug
|
|||
debug_i = 1, &
|
||||
debug_g = 1
|
||||
|
||||
integer, dimension(2), public :: &
|
||||
debug_stressMaxLocation = 0, &
|
||||
debug_stressMinLocation = 0, &
|
||||
debug_jacobianMaxLocation = 0, &
|
||||
debug_jacobianMinLocation = 0
|
||||
|
||||
|
||||
real(pReal), public :: &
|
||||
debug_stressMax = -huge(1.0_pReal), &
|
||||
debug_stressMin = huge(1.0_pReal), &
|
||||
debug_jacobianMax = -huge(1.0_pReal), &
|
||||
debug_jacobianMin = huge(1.0_pReal)
|
||||
|
||||
#ifdef PETSc
|
||||
character(len=1024), parameter, public :: &
|
||||
PETSCDEBUG = ' -snes_view -snes_monitor '
|
||||
#endif
|
||||
public :: debug_init, &
|
||||
debug_reset, &
|
||||
debug_info
|
||||
public :: debug_init
|
||||
|
||||
contains
|
||||
|
||||
|
|
@ -230,42 +215,4 @@ subroutine debug_init
|
|||
|
||||
end subroutine debug_init
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief resets all debug values
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
subroutine debug_reset
|
||||
|
||||
debug_stressMaxLocation = 0
|
||||
debug_stressMinLocation = 0
|
||||
debug_jacobianMaxLocation = 0
|
||||
debug_jacobianMinLocation = 0
|
||||
debug_stressMax = -huge(1.0_pReal)
|
||||
debug_stressMin = huge(1.0_pReal)
|
||||
debug_jacobianMax = -huge(1.0_pReal)
|
||||
debug_jacobianMin = huge(1.0_pReal)
|
||||
|
||||
end subroutine debug_reset
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief writes debug statements to standard out
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
subroutine debug_info
|
||||
|
||||
!$OMP CRITICAL (write2out)
|
||||
debugOutputCPFEM: if (iand(debug_level(debug_CPFEM),debug_LEVELBASIC) /= 0 &
|
||||
.and. any(debug_stressMinLocation /= 0) &
|
||||
.and. any(debug_stressMaxLocation /= 0) ) then
|
||||
write(6,'(2/,a,/)') ' Extreme values of returned stress and Jacobian'
|
||||
write(6,'(a39)') ' value el ip'
|
||||
write(6,'(a14,1x,e12.3,1x,i8,1x,i4)') ' stress min :', debug_stressMin, debug_stressMinLocation
|
||||
write(6,'(a14,1x,e12.3,1x,i8,1x,i4)') ' max :', debug_stressMax, debug_stressMaxLocation
|
||||
write(6,'(a14,1x,e12.3,1x,i8,1x,i4)') ' Jacobian min :', debug_jacobianMin, debug_jacobianMinLocation
|
||||
write(6,'(a14,1x,e12.3,1x,i8,1x,i4,/)') ' max :', debug_jacobianMax, debug_jacobianMaxLocation
|
||||
endif debugOutputCPFEM
|
||||
!$OMP END CRITICAL (write2out)
|
||||
|
||||
end subroutine debug_info
|
||||
|
||||
end module debug
|
||||
|
|
|
|||
|
|
@ -140,7 +140,8 @@ program DAMASK_grid
|
|||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
! reading information from load case file and to sanity checks
|
||||
fileContent = IO_read_ASCII(trim(loadCaseFile))
|
||||
fileContent = IO_readlines(trim(loadCaseFile))
|
||||
if(size(fileContent) == 0) call IO_error(307,ext_msg='No load case specified')
|
||||
|
||||
allocate (loadCases(0)) ! array of load cases
|
||||
do currentLoadCase = 1, size(fileContent)
|
||||
|
|
|
|||
|
|
@ -249,7 +249,8 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
|
|||
NiterationHomog = 0
|
||||
|
||||
cutBackLooping: do while (.not. terminallyIll .and. &
|
||||
any(subStep(:,FEsolving_execELem(1):FEsolving_execElem(2)) > num%subStepMinHomog))
|
||||
any(subStep(FEsolving_execIP(1):FEsolving_execIP(2),&
|
||||
FEsolving_execElem(1):FEsolving_execElem(2)) > num%subStepMinHomog))
|
||||
|
||||
!$OMP PARALLEL DO PRIVATE(myNgrains)
|
||||
elementLooping1: do e = FEsolving_execElem(1),FEsolving_execElem(2)
|
||||
|
|
|
|||
|
|
@ -45,9 +45,7 @@ contains
|
|||
!> @brief initializes the mesh by calling all necessary private routines the mesh module
|
||||
!! Order and routines strongly depend on type of solver
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
subroutine discretization_marc_init(ip,el)
|
||||
|
||||
integer, intent(in) :: el, ip
|
||||
subroutine discretization_marc_init
|
||||
|
||||
real(pReal), dimension(:,:), allocatable :: &
|
||||
node0_elem, & !< node x,y,z coordinates (initially!)
|
||||
|
|
@ -70,7 +68,7 @@ subroutine discretization_marc_init(ip,el)
|
|||
real(pReal), dimension(:,:,:,:),allocatable :: &
|
||||
unscaledNormals
|
||||
|
||||
write(6,'(/,a)') ' <<<+- mesh init -+>>>'; flush(6)
|
||||
write(6,'(/,a)') ' <<<+- discretization_marc init -+>>>'; flush(6)
|
||||
|
||||
mesh_unitlength = numerics_unitlength ! set physical extent of a length unit in mesh
|
||||
|
||||
|
|
@ -83,10 +81,6 @@ subroutine discretization_marc_init(ip,el)
|
|||
FEsolving_execElem = [1,nElems]
|
||||
FEsolving_execIP = [1,elem%nIPs]
|
||||
|
||||
allocate(calcMode(elem%nIPs,nElems),source=.false.) ! pretend to have collected what first call is asking (F = I)
|
||||
calcMode(ip,mesh_FEM2DAMASK_elem(el)) = .true. ! first ip,el needs to be already pingponged to "calc"
|
||||
|
||||
|
||||
allocate(cellNodeDefinition(elem%nNodes-1))
|
||||
allocate(connectivity_cell(elem%NcellNodesPerCell,elem%nIPs,nElems))
|
||||
call buildCells(connectivity_cell,cellNodeDefinition,&
|
||||
|
|
|
|||
|
|
@ -110,7 +110,7 @@ subroutine discretization_mesh_init(restart)
|
|||
call MPI_Bcast(mesh_boundaries,mesh_Nboundaries,MPI_INTEGER,0,PETSC_COMM_WORLD,ierr)
|
||||
|
||||
if (worldrank == 0) then
|
||||
fileContent = IO_read_ASCII(geometryFile)
|
||||
fileContent = IO_readlines(geometryFile)
|
||||
l = 0
|
||||
do
|
||||
l = l + 1
|
||||
|
|
|
|||
|
|
@ -28,8 +28,6 @@ module numerics
|
|||
numerics_unitlength = 1.0_pReal, & !< determines the physical length of one computational length unit
|
||||
charLength = 1.0_pReal, & !< characteristic length scale for gradient problems
|
||||
residualStiffness = 1.0e-6_pReal !< non-zero residual damage
|
||||
logical, protected, public :: &
|
||||
usePingPong = .true.
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
! field parameters:
|
||||
|
|
@ -133,8 +131,6 @@ subroutine numerics_init
|
|||
defgradTolerance = IO_floatValue(line,chunkPos,2)
|
||||
case ('ijacostiffness')
|
||||
iJacoStiffness = IO_intValue(line,chunkPos,2)
|
||||
case ('usepingpong')
|
||||
usepingpong = IO_intValue(line,chunkPos,2) > 0
|
||||
case ('unitlength')
|
||||
numerics_unitlength = IO_floatValue(line,chunkPos,2)
|
||||
|
||||
|
|
@ -221,7 +217,6 @@ subroutine numerics_init
|
|||
! writing parameters to output
|
||||
write(6,'(a24,1x,es8.1)') ' defgradTolerance: ',defgradTolerance
|
||||
write(6,'(a24,1x,i8)') ' iJacoStiffness: ',iJacoStiffness
|
||||
write(6,'(a24,1x,L8)') ' use ping pong scheme: ',usepingpong
|
||||
write(6,'(a24,1x,es8.1,/)')' unitlength: ',numerics_unitlength
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
|
|
|
|||
Loading…
Reference in New Issue