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

changed Line Ending to Unix style (LF)

This commit is contained in:
Philip Eisenlohr 2009-01-19 19:10:58 +00:00
parent abb2e3ef30
commit cbd7c279d4
14 changed files with 3157 additions and 3175 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

77
trunk/CPFEM_GIA8.f90
View File

@ -78,7 +78,8 @@
enddo enddo
! !
! *** Output to MARC output file *** ! *** Output to MARC output file ***
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,*) write(6,*)
write(6,*) 'CPFEM Initialization' write(6,*) 'CPFEM Initialization'
write(6,*) write(6,*)
@ -99,7 +100,8 @@
write(6,*) 'GIA_bNorm: ', shape(GIA_bNorm) write(6,*) 'GIA_bNorm: ', shape(GIA_bNorm)
write(6,*) write(6,*)
call flush(6) call flush(6)
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
return return
! !
END SUBROUTINE END SUBROUTINE
@ -158,28 +160,28 @@
endif endif
! !
cp_en = mesh_FEasCP('elem',CPFEM_en) cp_en = mesh_FEasCP('elem',CPFEM_en)
if (cp_en == 1 .and. CPFEM_in == 1) then if (cp_en == 1 .and. CPFEM_in == 1) then
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,'(a6,x,i4,x,a4,x,i4,x,a10,x,f8.4,x,a10,x,i2,x,a10,x,i2,x,a10,x,i2,x,a10,x,i2)') & write(6,'(a6,x,i4,x,a4,x,i4,x,a10,x,f8.4,x,a10,x,i2,x,a10,x,i2,x,a10,x,i2,x,a10,x,i2)') &
'elem',cp_en,'IP',CPFEM_in,& 'elem',cp_en,'IP',CPFEM_in,&
'theTime',theTime,'theInc',theInc,'theCycle',theCycle,'theLovl',theLovl,& 'theTime',theTime,'theInc',theInc,'theCycle',theCycle,'theLovl',theLovl,&
'mode',CPFEM_mode 'mode',CPFEM_mode
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
endif endif
! !
select case (CPFEM_mode) select case (CPFEM_mode)
case (2,1) ! regular computation (with aging of results) case (2,1) ! regular computation (with aging of results)
if (.not. CPFEM_calc_done) then ! puuh, me needs doing all the work... if (.not. CPFEM_calc_done) then ! puuh, me needs doing all the work...
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write (6,*) 'puuh me needs doing all the work', cp_en write (6,*) 'puuh me needs doing all the work', cp_en
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
if (CPFEM_mode == 1) then ! age results at start of new increment if (CPFEM_mode == 1) then ! age results at start of new increment
CPFEM_Fp_old = CPFEM_Fp_new CPFEM_Fp_old = CPFEM_Fp_new
constitutive_state_old = constitutive_state_new constitutive_state_old = constitutive_state_new
GIA_rVect_old = GIA_rVect_new GIA_rVect_old = GIA_rVect_new
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write (6,*) '#### aged results' write (6,*) '#### aged results'
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
endif endif
debug_cutbackDistribution = 0_pInt ! initialize debugging data debug_cutbackDistribution = 0_pInt ! initialize debugging data
debug_InnerLoopDistribution = 0_pInt debug_InnerLoopDistribution = 0_pInt
@ -195,9 +197,9 @@
CPFEM_calc_done = .true. ! now calc is done CPFEM_calc_done = .true. ! now calc is done
endif endif
! translate from P and dP/dF to CS and dCS/dE ! translate from P and dP/dF to CS and dCS/dE
!!$OMP CRITICAL (evilmatmul) !!$OMP CRITICAL (evilmatmul)
Kirchhoff_bar = math_mul33x33(CPFEM_PK1_bar(:,:,CPFEM_in, cp_en),transpose(CPFEM_ffn1_bar(:,:,CPFEM_in, cp_en))) Kirchhoff_bar = math_mul33x33(CPFEM_PK1_bar(:,:,CPFEM_in, cp_en),transpose(CPFEM_ffn1_bar(:,:,CPFEM_in, cp_en)))
!!$OMP END CRITICAL (evilmatmul) !!$OMP END CRITICAL (evilmatmul)
J_inverse = 1.0_pReal/math_det3x3(CPFEM_ffn1_bar(:,:,CPFEM_in, cp_en)) J_inverse = 1.0_pReal/math_det3x3(CPFEM_ffn1_bar(:,:,CPFEM_in, cp_en))
CPFEM_stress_bar(1:CPFEM_ngens,CPFEM_in,cp_en) = math_Mandel33to6(J_inverse*Kirchhoff_bar) CPFEM_stress_bar(1:CPFEM_ngens,CPFEM_in,cp_en) = math_Mandel33to6(J_inverse*Kirchhoff_bar)
! !
@ -208,8 +210,8 @@
math_I3(j,l)*CPFEM_ffn1_bar(i,m,CPFEM_in,cp_en)*CPFEM_PK1_bar(k,m,CPFEM_in,cp_en)) + & math_I3(j,l)*CPFEM_ffn1_bar(i,m,CPFEM_in,cp_en)*CPFEM_PK1_bar(k,m,CPFEM_in,cp_en)) + &
0.5_pReal*(math_I3(i,k)*Kirchhoff_bar(j,l) + math_I3(j,l)*Kirchhoff_bar(i,k) + & 0.5_pReal*(math_I3(i,k)*Kirchhoff_bar(j,l) + math_I3(j,l)*Kirchhoff_bar(i,k) + &
math_I3(i,l)*Kirchhoff_bar(j,k) + math_I3(j,k)*Kirchhoff_bar(i,l)) math_I3(i,l)*Kirchhoff_bar(j,k) + math_I3(j,k)*Kirchhoff_bar(i,l))
forall(i=1:3,j=1:3,k=1:3,l=1:3) & forall(i=1:3,j=1:3,k=1:3,l=1:3) &
H_bar_sym(i,j,k,l)= 0.25_pReal*(H_bar(i,j,k,l) +H_bar(j,i,k,l) +H_bar(i,j,l,k) +H_bar(j,i,l,k)) H_bar_sym(i,j,k,l)= 0.25_pReal*(H_bar(i,j,k,l) +H_bar(j,i,k,l) +H_bar(i,j,l,k) +H_bar(j,i,l,k))
CPFEM_jaco_bar(1:CPFEM_ngens,1:CPFEM_ngens,CPFEM_in,cp_en) = math_Mandel3333to66(J_inverse*H_bar) CPFEM_jaco_bar(1:CPFEM_ngens,1:CPFEM_ngens,CPFEM_in,cp_en) = math_Mandel3333to66(J_inverse*H_bar)
! !
case (3) ! collect and return odd result case (3) ! collect and return odd result
@ -219,7 +221,6 @@
CPFEM_stress_bar(1:CPFEM_ngens,CPFEM_in,cp_en) = CPFEM_odd_stress CPFEM_stress_bar(1:CPFEM_ngens,CPFEM_in,cp_en) = CPFEM_odd_stress
CPFEM_jaco_bar(1:CPFEM_ngens,1:CPFEM_ngens,CPFEM_in,cp_en) = CPFEM_odd_jacobian*math_identity2nd(CPFEM_ngens) CPFEM_jaco_bar(1:CPFEM_ngens,1:CPFEM_ngens,CPFEM_in,cp_en) = CPFEM_odd_jacobian*math_identity2nd(CPFEM_ngens)
CPFEM_calc_done = .false. CPFEM_calc_done = .false.
case (4) ! do nothing since we can recycle the former results (MARC specialty) case (4) ! do nothing since we can recycle the former results (MARC specialty)
case (5) ! record consistent tangent at beginning of new increment case (5) ! record consistent tangent at beginning of new increment
CPFEM_jaco_knownGood = CPFEM_jaco_bar CPFEM_jaco_knownGood = CPFEM_jaco_bar
@ -314,15 +315,15 @@
! -------------- grain loop ----------------- ! -------------- grain loop -----------------
do grain = 1,texture_Ngrains(mesh_element(4,cp_en)) do grain = 1,texture_Ngrains(mesh_element(4,cp_en))
call SingleCrystallite(msg,PK1(:,:,grain),dPdF(:,:,:,:,grain),& call SingleCrystallite(msg,PK1(:,:,grain),dPdF(:,:,:,:,grain),&
CPFEM_results(CPFEM_Nresults+1:CPFEM_Nresults+constitutive_Nresults(grain,CPFEM_in,cp_en),& CPFEM_results(CPFEM_Nresults+1:CPFEM_Nresults+constitutive_Nresults(grain,CPFEM_in,cp_en),&
grain,CPFEM_in,cp_en),& grain,CPFEM_in,cp_en),&
Fp1(:,:,grain),Fe1(:,:,grain),state1(:,grain),& ! output up to here Fp1(:,:,grain),Fe1(:,:,grain),state1(:,grain),& ! output up to here
dTime,cp_en,CPFEM_in,grain,.true.,& dTime,cp_en,CPFEM_in,grain,.true.,&
CPFEM_Temperature(CPFEM_in,cp_en),F1(:,:,grain),F0(:,:,grain),Fp0(:,:,grain),state0(:,grain)) CPFEM_Temperature(CPFEM_in,cp_en),F1(:,:,grain),F0(:,:,grain),Fp0(:,:,grain),state0(:,grain))
if (msg /= 'ok') then ! solution not reached --> exit NRIteration if (msg /= 'ok') then ! solution not reached --> exit NRIteration
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,*) 'GIA: grain loop failed to converge @ EL:',cp_en,' IP:',CPFEM_in write(6,*) 'GIA: grain loop failed to converge @ EL:',cp_en,' IP:',CPFEM_in
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
NRconvergent = .false. NRconvergent = .false.
exit NRiteration exit NRiteration
endif endif
@ -379,10 +380,10 @@
enddo enddo
resNorm = sqrt(resNorm) resNorm = sqrt(resNorm)
! !
if (debugger) then if (debugger) then
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,'(x,a,i3,a,i3,a,i3,a,e10.4)')'EL:',cp_en,' IP:',CPFEM_in,' Iter:',NRiter,' RNorm:',resNorm write(6,'(x,a,i3,a,i3,a,i3,a,e10.4)')'EL:',cp_en,' IP:',CPFEM_in,' Iter:',NRiter,' RNorm:',resNorm
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
if (NRiter == 1_pInt) resMax = resNorm if (NRiter == 1_pInt) resMax = resNorm
if ((resNorm < resToler*resMax) .or. (resNorm < resAbsol)) then ! resNorm < tolerance ===> convergent if ((resNorm < resToler*resMax) .or. (resNorm < resAbsol)) then ! resNorm < tolerance ===> convergent
NRconvergent = .true. NRconvergent = .true.
@ -422,9 +423,9 @@
dvardres = 0.0_pReal dvardres = 0.0_pReal
call math_invert(36,dresdvar,dvardres,dummy,failed) call math_invert(36,dresdvar,dvardres,dummy,failed)
if (failed) then if (failed) then
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,*) 'GIA: failed to invert the Jacobian @ EL:',cp_en,' IP:',CPFEM_in write(6,*) 'GIA: failed to invert the Jacobian @ EL:',cp_en,' IP:',CPFEM_in
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
NRconvergent = .false. NRconvergent = .false.
exit NRiteration exit NRiteration
endif endif
@ -452,9 +453,9 @@
! !
! return to the general subroutine when convergence is not reached ! return to the general subroutine when convergence is not reached
if (.not. NRconvergent) then if (.not. NRconvergent) then
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,'(x,a)') 'GIA: convergence is not reached @ EL:',cp_en,' IP:',CPFEM_in write(6,'(x,a)') 'GIA: convergence is not reached @ EL:',cp_en,' IP:',CPFEM_in
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
call IO_error(600) call IO_error(600)
return return
endif endif
@ -473,13 +474,13 @@
! update results plotted in MENTAT ! update results plotted in MENTAT
call math_pDecomposition(Fe1(:,:,grain),U,R,error) ! polar decomposition call math_pDecomposition(Fe1(:,:,grain),U,R,error) ! polar decomposition
if (error) then if (error) then
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,*) Fe1(:,:,grain) write(6,*) Fe1(:,:,grain)
write(6,*) 'polar decomposition' write(6,*) 'polar decomposition'
write(6,*) 'Grain: ',grain write(6,*) 'Grain: ',grain
write(6,*) 'Integration point: ',CPFEM_in write(6,*) 'Integration point: ',CPFEM_in
write(6,*) 'Element: ',mesh_element(1,cp_en) write(6,*) 'Element: ',mesh_element(1,cp_en)
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
call IO_error(650) call IO_error(650)
return return
endif endif
@ -503,15 +504,15 @@
call GIA_RelaxedDeformation(F1,F1_bar,rx) call GIA_RelaxedDeformation(F1,F1_bar,rx)
do grain = 1,8 do grain = 1,8
call SingleCrystallite(msg,PK1(:,:,grain),dPdF(:,:,:,:,grain),& call SingleCrystallite(msg,PK1(:,:,grain),dPdF(:,:,:,:,grain),&
CPFEM_results(CPFEM_Nresults+1:CPFEM_Nresults+constitutive_Nresults(grain,CPFEM_in,cp_en),& CPFEM_results(CPFEM_Nresults+1:CPFEM_Nresults+constitutive_Nresults(grain,CPFEM_in,cp_en),&
grain,CPFEM_in,cp_en),& grain,CPFEM_in,cp_en),&
Fp1(:,:,grain),Fe1(:,:,grain),state1(:,grain),& ! output up to here Fp1(:,:,grain),Fe1(:,:,grain),state1(:,grain),& ! output up to here
dTime,cp_en,CPFEM_in,grain,.true.,& dTime,cp_en,CPFEM_in,grain,.true.,&
CPFEM_Temperature(CPFEM_in,cp_en),F1(:,:,grain),F0(:,:,grain),Fp0(:,:,grain),state0(:,grain)) CPFEM_Temperature(CPFEM_in,cp_en),F1(:,:,grain),F0(:,:,grain),Fp0(:,:,grain),state0(:,grain))
if (msg /= 'ok') then ! solution not reached --> exit NRIteration if (msg /= 'ok') then ! solution not reached --> exit NRIteration
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,*) 'GIA: perturbation grain loop failed to converge within allowable step-size' write(6,*) 'GIA: perturbation grain loop failed to converge within allowable step-size'
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
NRconvergent = .false. NRconvergent = .false.
exit NRPerturbation exit NRPerturbation
endif endif
@ -563,11 +564,11 @@
enddo enddo
resNorm = sqrt(resNorm) resNorm = sqrt(resNorm)
! !
! if (debugger) then ! if (debugger) then
!!$OMP CRITICAL (write2out) !!$OMP CRITICAL (write2out)
! write(6,'(x,a,i3,a,i3,a,i3,a,i3,a,e10.4)')'EL = ',cp_en,':IP = ',CPFEM_in,':pert = ',3*(ip-1)+jp,':Iter = ',NRiter,':RNorm = ',resNorm ! write(6,'(x,a,i3,a,i3,a,i3,a,i3,a,e10.4)')'EL = ',cp_en,':IP = ',CPFEM_in,':pert = ',3*(ip-1)+jp,':Iter = ',NRiter,':RNorm = ',resNorm
!!$OMP END CRITICAL (write2out) !!$OMP END CRITICAL (write2out)
! endif ! endif
if (NRiter == 1_pInt) resMax = resNorm if (NRiter == 1_pInt) resMax = resNorm
if ((resNorm < resToler*resMax) .or. (resNorm < resAbsol)) then ! resNorm < tolerance ===> convergent if ((resNorm < resToler*resMax) .or. (resNorm < resAbsol)) then ! resNorm < tolerance ===> convergent
NRconvergent = .true. NRconvergent = .true.
@ -607,9 +608,9 @@
dvardres = 0.0_pReal dvardres = 0.0_pReal
call math_invert(36,dresdvar,dvardres,dummy,failed) call math_invert(36,dresdvar,dvardres,dummy,failed)
if (failed) then if (failed) then
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,*) 'GIA: perturbation failed to invert the Jacobian' write(6,*) 'GIA: perturbation failed to invert the Jacobian'
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
NRconvergent = .false. NRconvergent = .false.
exit NRPerturbation exit NRPerturbation
endif endif

1992
trunk/CPFEM_Taylor.f90
View File

File diff suppressed because it is too large Load Diff

1405
trunk/CPFEM_Taylor_sequential.f90
View File

File diff suppressed because it is too large Load Diff

79
trunk/FEsolving.f90
View File

@ -10,44 +10,45 @@
integer(pInt) theInc,theCycle,theLovl integer(pInt) theInc,theCycle,theLovl
real(pReal) theTime real(pReal) theTime
logical :: lastIncConverged = .false.,outdatedByNewInc = .false., outdatedFFN1 = .false. logical :: lastIncConverged = .false.,outdatedByNewInc = .false., outdatedFFN1 = .false.
logical :: symmetricSolver = .false. logical :: symmetricSolver = .false.
CONTAINS
CONTAINS
!***********************************************************
! determine wether a symmetric solver is used !***********************************************************
!*********************************************************** ! determine wether a symmetric solver is used
subroutine FE_get_solverSymmetry(unit) !***********************************************************
subroutine FE_get_solverSymmetry(unit)
use prec, only: pInt
use IO use prec, only: pInt
implicit none use IO
implicit none
integer(pInt) unit
integer(pInt), dimension (133) :: pos integer(pInt) unit
character*300 line integer(pInt), dimension (133) :: pos
character*300 line
610 FORMAT(A300)
610 FORMAT(A300)
rewind(unit)
do rewind(unit)
read (unit,610,END=630) line do
pos = IO_stringPos(line,1) read (unit,610,END=630) line
if( IO_lc(IO_stringValue(line,pos,1)) == 'solver' ) then pos = IO_stringPos(line,1)
read (unit,610,END=630) line ! Garbage line if( IO_lc(IO_stringValue(line,pos,1)) == 'solver' ) then
pos = IO_stringPos(line,2) ! limit to 64 nodes max (plus ID, type) read (unit,610,END=630) line ! Garbage line
if(IO_intValue(line,pos,2) /= 1_pInt) then pos = IO_stringPos(line,2) ! limit to 64 nodes max (plus ID, type)
symmetricSolver = .true. if(IO_intValue(line,pos,2) /= 1_pInt) then
!$OMP CRITICAL (write2out) symmetricSolver = .true.
write (6,*) !$OMP CRITICAL (write2out)
write (6,*) 'Symmetric solver detected. d-Matrix will be symmetrized!' write (6,*)
!$OMP END CRITICAL (write2out) write (6,*) 'Symmetric solver detected. d-Matrix will be symmetrized!'
endif !$OMP END CRITICAL (write2out)
endif endif
enddo endif
enddo
630 return
630 return
end subroutine
end subroutine
END MODULE FEsolving END MODULE FEsolving

27
trunk/IO.f90
View File

@ -256,8 +256,8 @@
character(len=*), intent(in) :: line character(len=*), intent(in) :: line
character(len=*), parameter :: sep=achar(32)//achar(9)//achar(10)//achar(13) ! whitespaces character(len=*), parameter :: sep=achar(32)//achar(9)//achar(10)//achar(13) ! whitespaces
integer(pInt), intent(in) :: N integer(pInt), intent(in) :: N
integer(pInt) part integer(pInt) part
integer(pInt) IO_stringPos(1+N*2) integer(pInt) IO_stringPos(1+N*2)
IO_stringPos = -1 IO_stringPos = -1
@ -324,7 +324,7 @@
character(len=*), intent(in) :: line character(len=*), intent(in) :: line
integer(pInt), intent(in) :: positions(*),pos integer(pInt), intent(in) :: positions(*),pos
real(pReal) IO_floatValue real(pReal) IO_floatValue
if (positions(1) >= pos) then if (positions(1) >= pos) then
read(UNIT=line(positions(pos*2):positions(pos*2+1)),ERR=100,FMT=*) IO_floatValue read(UNIT=line(positions(pos*2):positions(pos*2+1)),ERR=100,FMT=*) IO_floatValue
@ -346,7 +346,7 @@
character(len=*), intent(in) :: line character(len=*), intent(in) :: line
integer(pInt), intent(in) :: ends(*),pos integer(pInt), intent(in) :: ends(*),pos
real(pReal) IO_fixedFloatValue real(pReal) IO_fixedFloatValue
read(UNIT=line(ends(pos-1)+1:ends(pos)),ERR=100,FMT=*) IO_fixedFloatValue read(UNIT=line(ends(pos-1)+1:ends(pos)),ERR=100,FMT=*) IO_fixedFloatValue
return return
@ -365,9 +365,9 @@
implicit none implicit none
character(len=*), intent(in) :: line character(len=*), intent(in) :: line
integer(pInt), intent(in) :: ends(*),pos integer(pInt), intent(in) :: ends(*),pos
integer(pInt) pos_exp,expon integer(pInt) pos_exp,expon
real(pReal) IO_fixedNoEFloatValue,base real(pReal) IO_fixedNoEFloatValue,base
pos_exp = scan(line(ends(pos)+1:ends(pos+1)),'+-',back=.true.) pos_exp = scan(line(ends(pos)+1:ends(pos+1)),'+-',back=.true.)
if (pos_exp > 1) then if (pos_exp > 1) then
@ -395,7 +395,7 @@
character(len=*), intent(in) :: line character(len=*), intent(in) :: line
integer(pInt), intent(in) :: positions(*),pos integer(pInt), intent(in) :: positions(*),pos
integer(pInt) IO_intValue integer(pInt) IO_intValue
if (positions(1) >= pos) then if (positions(1) >= pos) then
read(UNIT=line(positions(pos*2):positions(pos*2+1)),ERR=100,FMT=*) IO_intValue read(UNIT=line(positions(pos*2):positions(pos*2+1)),ERR=100,FMT=*) IO_intValue
@ -417,7 +417,7 @@
character(len=*), intent(in) :: line character(len=*), intent(in) :: line
integer(pInt), intent(in) :: ends(*),pos integer(pInt), intent(in) :: ends(*),pos
integer(pInt) IO_fixedIntValue integer(pInt) IO_fixedIntValue
read(UNIT=line(ends(pos)+1:ends(pos+1)),ERR=100,FMT=*) IO_fixedIntValue read(UNIT=line(ends(pos)+1:ends(pos+1)),ERR=100,FMT=*) IO_fixedIntValue
return return
@ -595,19 +595,18 @@
case default case default
msg='Unknown error number' msg='Unknown error number'
end select end select
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,*) 'MPIE Material Routine Ver. 0.0 by the coding team' write(6,*) 'MPIE Material Routine Ver. 0.0 by the coding team'
write(6,*) write(6,*)
write(6,*) msg write(6,*) msg
write(6,*) write(6,*)
call debug_info() call debug_info()
call flush(6) call flush(6)
call quit(9000+ID) call quit(9000+ID)
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
! ABAQUS returns in some cases ! ABAQUS returns in some cases
return return

372
trunk/concom2008r1
View File

@ -1,186 +1,186 @@
! reformated to free format ! reformated to free format
!*********************************************************************** !***********************************************************************
! !
! File: concom.cmn ! File: concom.cmn
! !
! MSC.Marc include file ! MSC.Marc include file
! !
integer(pInt) & integer(pInt) &
iacous, iasmbl, iautth, ibear, icompl, iconj, icreep, ideva, idyn, idynt,& iacous, iasmbl, iautth, ibear, icompl, iconj, icreep, ideva, idyn, idynt,&
ielas, ielcma, ielect, iform, ifour, iharm, ihcps, iheat, iheatt, ihresp,& ielas, ielcma, ielect, iform, ifour, iharm, ihcps, iheat, iheatt, ihresp,&
ijoule, ilem, ilnmom, iloren, inc, incext, incsub, ipass, iplres, ipois,& ijoule, ilem, ilnmom, iloren, inc, incext, incsub, ipass, iplres, ipois,&
ipoist, irpflo, ismall, ismalt, isoil, ispect, ispnow, istore, iswep, ithcrp,& ipoist, irpflo, ismall, ismalt, isoil, ispect, ispnow, istore, iswep, ithcrp,&
itherm, iupblg, iupdat, jacflg, jel, jparks, largst, lfond, loadup, loaduq,& itherm, iupblg, iupdat, jacflg, jel, jparks, largst, lfond, loadup, loaduq,&
lodcor, lovl, lsub, magnet, ncycle, newtnt, newton, noshr, linear, ivscpl,& lodcor, lovl, lsub, magnet, ncycle, newtnt, newton, noshr, linear, ivscpl,&
icrpim, iradrt, ipshft, itshr, iangin, iupmdr, iconjf, jincfl, jpermg, jhour,& icrpim, iradrt, ipshft, itshr, iangin, iupmdr, iconjf, jincfl, jpermg, jhour,&
isolvr, jritz, jtable, jshell, jdoubl, jform, jcentr, imini, kautth, iautof,& isolvr, jritz, jtable, jshell, jdoubl, jform, jcentr, imini, kautth, iautof,&
ibukty, iassum, icnstd, icnstt, kmakmas, imethvp,iradrte,iradrtp, iupdate,iupdatp,& ibukty, iassum, icnstd, icnstt, kmakmas, imethvp,iradrte,iradrtp, iupdate,iupdatp,&
ncycnt, marmen ,idynme, ihavca, ispf, kmini, imixed, largtt, kdoela, iautofg,& ncycnt, marmen ,idynme, ihavca, ispf, kmini, imixed, largtt, kdoela, iautofg,&
ipshftp,idntrc, ipore, jtablm, jtablc, isnecma,itrnspo,imsdif, jtrnspo,mcnear,& ipshftp,idntrc, ipore, jtablm, jtablc, isnecma,itrnspo,imsdif, jtrnspo,mcnear,&
imech, imecht, ielcmat, ielectt,magnett, imsdift,noplas, jtabls, jactch, jtablth,& imech, imecht, ielcmat, ielectt,magnett, imsdift,noplas, jtabls, jactch, jtablth,&
kgmsto ,jpzo, ifricsh, iremkin,iremfor, ishearp,jspf, machining, jlshell,icompsol,& kgmsto ,jpzo, ifricsh, iremkin,iremfor, ishearp,jspf, machining, jlshell,icompsol,&
iupblgfo,jcondir,nstcrp, nactive,ipassref, nstspnt,ibeart,icheckmpc, noline, icuring,& iupblgfo,jcondir,nstcrp, nactive,ipassref, nstspnt,ibeart,icheckmpc, noline, icuring,&
ishrink,ioffsflg,isetoff, ioffsetm,iharmt, inc_incdat,iautspc,ibrake, icbush ,istream_input,& ishrink,ioffsflg,isetoff, ioffsetm,iharmt, inc_incdat,iautspc,ibrake, icbush ,istream_input,&
iprsinp,ivlsinp,ifirst_time,ipin_m,jgnstr_glb,imarc_return,iqvcinp,nqvceid,istpnx,imicro1,& iprsinp,ivlsinp,ifirst_time,ipin_m,jgnstr_glb,imarc_return,iqvcinp,nqvceid,istpnx,imicro1,&
iaxisymm,jbreakglue,iglstif,jfastasm,iwear, iwearcf, imixmeth,ielcmadyn,idinout,igena_meth iaxisymm,jbreakglue,iglstif,jfastasm,iwear, iwearcf, imixmeth,ielcmadyn,idinout,igena_meth
integer(pInt) num_concom integer(pInt) num_concom
parameter(num_concom=219) parameter(num_concom=219)
common/marc_concom/& common/marc_concom/&
iacous, iasmbl, iautth, ibear, icompl, iconj, icreep, ideva(50), idyn, idynt,& iacous, iasmbl, iautth, ibear, icompl, iconj, icreep, ideva(50), idyn, idynt,&
ielas, ielcma, ielect, iform, ifour, iharm, ihcps, iheat, iheatt, ihresp,& ielas, ielcma, ielect, iform, ifour, iharm, ihcps, iheat, iheatt, ihresp,&
ijoule, ilem, ilnmom, iloren, inc, incext, incsub, ipass, iplres, ipois,& ijoule, ilem, ilnmom, iloren, inc, incext, incsub, ipass, iplres, ipois,&
ipoist, irpflo, ismall, ismalt, isoil, ispect, ispnow, istore, iswep, ithcrp,& ipoist, irpflo, ismall, ismalt, isoil, ispect, ispnow, istore, iswep, ithcrp,&
itherm, iupblg, iupdat, jacflg, jel, jparks, largst, lfond, loadup, loaduq,& itherm, iupblg, iupdat, jacflg, jel, jparks, largst, lfond, loadup, loaduq,&
lodcor, lovl, lsub, magnet, ncycle, newtnt, newton, noshr, linear, ivscpl,& lodcor, lovl, lsub, magnet, ncycle, newtnt, newton, noshr, linear, ivscpl,&
icrpim, iradrt, ipshft, itshr, iangin, iupmdr, iconjf, jincfl, jpermg, jhour,& icrpim, iradrt, ipshft, itshr, iangin, iupmdr, iconjf, jincfl, jpermg, jhour,&
isolvr, jritz, jtable, jshell, jdoubl, jform, jcentr, imini, kautth, iautof,& isolvr, jritz, jtable, jshell, jdoubl, jform, jcentr, imini, kautth, iautof,&
ibukty, iassum, icnstd, icnstt, kmakmas, imethvp,iradrte,iradrtp, iupdate,iupdatp,& ibukty, iassum, icnstd, icnstt, kmakmas, imethvp,iradrte,iradrtp, iupdate,iupdatp,&
ncycnt, marmen ,idynme, ihavca, ispf, kmini, imixed, largtt, kdoela, iautofg,& ncycnt, marmen ,idynme, ihavca, ispf, kmini, imixed, largtt, kdoela, iautofg,&
ipshftp,idntrc, ipore, jtablm, jtablc, isnecma,itrnspo,imsdif, jtrnspo,mcnear,& ipshftp,idntrc, ipore, jtablm, jtablc, isnecma,itrnspo,imsdif, jtrnspo,mcnear,&
imech, imecht, ielcmat, ielectt,magnett, imsdift,noplas, jtabls, jactch, jtablth,& imech, imecht, ielcmat, ielectt,magnett, imsdift,noplas, jtabls, jactch, jtablth,&
kgmsto ,jpzo, ifricsh, iremkin,iremfor, ishearp,jspf, machining, jlshell,icompsol,& kgmsto ,jpzo, ifricsh, iremkin,iremfor, ishearp,jspf, machining, jlshell,icompsol,&
iupblgfo,jcondir,nstcrp, nactive,ipassref, nstspnt,ibeart,icheckmpc, noline, icuring,& iupblgfo,jcondir,nstcrp, nactive,ipassref, nstspnt,ibeart,icheckmpc, noline, icuring,&
ishrink,ioffsflg,isetoff, ioffsetm,iharmt, inc_incdat,iautspc,ibrake, icbush ,istream_input,& ishrink,ioffsflg,isetoff, ioffsetm,iharmt, inc_incdat,iautspc,ibrake, icbush ,istream_input,&
iprsinp,ivlsinp,ifirst_time,ipin_m,jgnstr_glb,imarc_return,iqvcinp,nqvceid,istpnx,imicro1,& iprsinp,ivlsinp,ifirst_time,ipin_m,jgnstr_glb,imarc_return,iqvcinp,nqvceid,istpnx,imicro1,&
iaxisymm,jbreakglue,iglstif,jfastasm,iwear, iwearcf, imixmeth, ielcmadyn,idinout,igena_meth iaxisymm,jbreakglue,iglstif,jfastasm,iwear, iwearcf, imixmeth, ielcmadyn,idinout,igena_meth
! !
! comments of variables: ! comments of variables:
! !
! ideva(50) - debug print out flag ! ideva(50) - debug print out flag
! 1 print element stiffness matrices, mass matrix ! 1 print element stiffness matrices, mass matrix
! 2 output matrices used in tying ! 2 output matrices used in tying
! 3 force the solution of a nonpositive definite matrix ! 3 force the solution of a nonpositive definite matrix
! 4 print info of connections to each node ! 4 print info of connections to each node
! 5 info of gap convergence, internal heat generated, contact ! 5 info of gap convergence, internal heat generated, contact
! touching and separation ! touching and separation
! 6 nodal value array during rezoning ! 6 nodal value array during rezoning
! 7 tying info in CONRAD GAP option, fluid element numbers in ! 7 tying info in CONRAD GAP option, fluid element numbers in
! CHANNEL option ! CHANNEL option
! 8 output incremental displacements in local coord. system ! 8 output incremental displacements in local coord. system
! 9 latent heat output ! 9 latent heat output
! 10 stress-strain in local coord. system ! 10 stress-strain in local coord. system
! 11 additional info on interlaminar stress ! 11 additional info on interlaminar stress
! 12 output right hand side and solution vector ! 12 output right hand side and solution vector
! 13 info of CPU resources used and memory available on NT ! 13 info of CPU resources used and memory available on NT
! 14 info of mesh adaption process, 2D outline information ! 14 info of mesh adaption process, 2D outline information
! info of penetration checking for remeshing ! info of penetration checking for remeshing
! save .fem files after afmesh3d meshing ! save .fem files after afmesh3d meshing
! 15 surface energy balance flag ! 15 surface energy balance flag
! 16 print info regarding pyrolysis ! 16 print info regarding pyrolysis
! 17 print info of "streamline topology" ! 17 print info of "streamline topology"
! 18 print mesh data changes after remeshing ! 18 print mesh data changes after remeshing
! 19 print material flow stress data read in from *.mat file ! 19 print material flow stress data read in from *.mat file
! if unit flag is on, print out flow stress after conversion ! if unit flag is on, print out flow stress after conversion
! 20 print information on table input ! 20 print information on table input
! 21 print out information regarding kinematic boundary conditions ! 21 print out information regarding kinematic boundary conditions
! 22 print out information regarding dist loads, point loads, film ! 22 print out information regarding dist loads, point loads, film
! and foundations ! and foundations
! 23 print out information about automatic domain decomposition ! 23 print out information about automatic domain decomposition
! 24 print out iteration information in SuperForm status report file ! 24 print out iteration information in SuperForm status report file
! 25 print out information for ablation ! 25 print out information for ablation
! 26 print out information for films - Table input ! 26 print out information for films - Table input
! 27 print out the tying forces ! 27 print out the tying forces
! 28 print out for CASI solver, convection, ! 28 print out for CASI solver, convection,
! 29 DDM single file debug printout ! 29 DDM single file debug printout
! 30 print out cavity debug info ! 30 print out cavity debug info
! 31 print out welding related info ! 31 print out welding related info
! 32 prints categorized DDM memory usage ! 32 prints categorized DDM memory usage
! 33 print out the cutting info regarding machining feature ! 33 print out the cutting info regarding machining feature
! 34 print out the list of quantities which can be defined via a table ! 34 print out the list of quantities which can be defined via a table
! and for each quantity the supported independent variables ! and for each quantity the supported independent variables
! 35 print out detailed coupling region info ! 35 print out detailed coupling region info
! 36 print out solver debug info level 1 (Least Detailed) ! 36 print out solver debug info level 1 (Least Detailed)
! 37 print out solver debug info level 1 (Medium Detailed) ! 37 print out solver debug info level 1 (Medium Detailed)
! 38 print out solver debug info level 1 (Very Detailed) ! 38 print out solver debug info level 1 (Very Detailed)
! 39 print detailed memory allocation info ! 39 print detailed memory allocation info
! 40 print out marc-adams debug info ! 40 print out marc-adams debug info
! 41 output rezone mapping post file for debugging ! 41 output rezone mapping post file for debugging
! 42 output post file after calling oprofos() for debugging ! 42 output post file after calling oprofos() for debugging
! 43 debug printout for vcct ! 43 debug printout for vcct
! 44 debug printout for progressive failure ! 44 debug printout for progressive failure
! 45 print out automatically generated midside node coordinates (arecrd) ! 45 print out automatically generated midside node coordinates (arecrd)
! 46 print out message about routine and location, where the ibort is raised (ibort_inc) ! 46 print out message about routine and location, where the ibort is raised (ibort_inc)
! 47 print out summary message of element variables on a ! 47 print out summary message of element variables on a
! group-basis after all the automatic changes have been ! group-basis after all the automatic changes have been
! made (em_ellibp) ! made (em_ellibp)
! 48 Automatically generate check results based on max and min vals. ! 48 Automatically generate check results based on max and min vals.
! These vals are stored in the checkr file, which is inserted ! These vals are stored in the checkr file, which is inserted
! into the *dat file by the generate_check_results script from /marc/tools ! into the *dat file by the generate_check_results script from /marc/tools
! 49 Automatically generate check results based on the real calculated values ! 49 Automatically generate check results based on the real calculated values
! at the sppecified check result locations. ! at the sppecified check result locations.
! These vals are stored in the checkr file, which is inserted ! These vals are stored in the checkr file, which is inserted
! into the *dat file by the update_check_results script from /marc/tools ! into the *dat file by the update_check_results script from /marc/tools
! !
! !
! jactch = 1 or 2 if elements are activated or deactivated ! jactch = 1 or 2 if elements are activated or deactivated
! = 3 if elements are adaptively remeshed or rezoned ! = 3 if elements are adaptively remeshed or rezoned
! = 0 normally / reset to 0 when assembly is done ! = 0 normally / reset to 0 when assembly is done
! ifricsh = 0 call to fricsh in otest not needed ! ifricsh = 0 call to fricsh in otest not needed
! = 1 call to fricsh (nodal friction) in otest needed ! = 1 call to fricsh (nodal friction) in otest needed
! iremkin = 0 remove deactivated kinematic boundary conditions ! iremkin = 0 remove deactivated kinematic boundary conditions
! immediately - only in new input format (this is default) ! immediately - only in new input format (this is default)
! = 1 remove deactivated kinematic boundary conditions ! = 1 remove deactivated kinematic boundary conditions
! gradually - only in new input format ! gradually - only in new input format
! iremfor = 0 remove force boundary conditions immediately - ! iremfor = 0 remove force boundary conditions immediately -
! only in new input format (this is default) ! only in new input format (this is default)
! = 1 remove force boundary conditions gradually - ! = 1 remove force boundary conditions gradually -
! only in new input format (this is default) ! only in new input format (this is default)
! ishearp set to 1 if shear panel elements are present in the model ! ishearp set to 1 if shear panel elements are present in the model
! !
! jspf = 0 not in spf loadcase ! jspf = 0 not in spf loadcase
! > 0 in spf loadcase (jspf=1 during first increment) ! > 0 in spf loadcase (jspf=1 during first increment)
! machining = 1 if the metal cutting feature is used, for memory allocation purpose ! machining = 1 if the metal cutting feature is used, for memory allocation purpose
! = 0 (default) if no metal cutting feature required ! = 0 (default) if no metal cutting feature required
! !
! jlshell = 1 if there is a shell element in the mesh ! jlshell = 1 if there is a shell element in the mesh
! icompsol = 1 if there is a composite solid element in the mesh ! icompsol = 1 if there is a composite solid element in the mesh
! iupblgfo = 1 if follower force for point loads ! iupblgfo = 1 if follower force for point loads
! jcondir = 1 if contact priority option is used ! jcondir = 1 if contact priority option is used
! nstcrp = 0 (default) steady state creep flag (undocumented feature. ! nstcrp = 0 (default) steady state creep flag (undocumented feature.
! if not 0, turns off special ncycle = 0 code in radial.f) ! if not 0, turns off special ncycle = 0 code in radial.f)
! nactive = number of active passes, if =1 then it's not a coupled analysis ! nactive = number of active passes, if =1 then it's not a coupled analysis
! ipassref = reference ipass, if not in a multiphysics pass ipass=ipassref ! ipassref = reference ipass, if not in a multiphysics pass ipass=ipassref
! icheckmpc = value of mpc-check parameter option ! icheckmpc = value of mpc-check parameter option
! noline = set to 1 in osolty if no line seacrh should be done in ogetst ! noline = set to 1 in osolty if no line seacrh should be done in ogetst
! icuring = set to 1 if the curing is included for the heat transfer analysis. ! icuring = set to 1 if the curing is included for the heat transfer analysis.
! ishrink = set to 1 if shrinkage strain is included for mechancial analysis. ! ishrink = set to 1 if shrinkage strain is included for mechancial analysis.
! ioffsflg = 1 for small displacement beam/shell offsets ! ioffsflg = 1 for small displacement beam/shell offsets
! = 2 for large displacement beam/shell offsets ! = 2 for large displacement beam/shell offsets
! isetoff = 0 - do not apply beam/shell offsets ! isetoff = 0 - do not apply beam/shell offsets
! = 1 - apply beam/shell offsets ! = 1 - apply beam/shell offsets
! ioffsetm = min. value of offset flag ! ioffsetm = min. value of offset flag
! inc_incdat = flag to record increment number of a new loadcase in incdat.f ! inc_incdat = flag to record increment number of a new loadcase in incdat.f
! iautspc = flag for AutoSPC option ! iautspc = flag for AutoSPC option
! ibrake = brake squeal in this increment ! ibrake = brake squeal in this increment
! icbush = set to 1 if cbush elements present in model ! icbush = set to 1 if cbush elements present in model
! istream_input = set to 1 for streaming input calling Marc as library ! istream_input = set to 1 for streaming input calling Marc as library
! iprsinp = set to 1 if pressure input, introduced so other variables ! iprsinp = set to 1 if pressure input, introduced so other variables
! such as h could be a function of pressure ! such as h could be a function of pressure
! ivlsinp = set to 1 if velocity input, introduced so other variables ! ivlsinp = set to 1 if velocity input, introduced so other variables
! such as h could be a function of velocity ! such as h could be a function of velocity
! ipin_m = # of beam element with PIN flag ! ipin_m = # of beam element with PIN flag
! jgnstr_glb = global control over pre or fast integrated composite shells ! jgnstr_glb = global control over pre or fast integrated composite shells
! imarc_return = Marc return flag for streaming input control ! imarc_return = Marc return flag for streaming input control
! iqvcimp = if non-zero, then the number of QVECT boundary conditions ! iqvcimp = if non-zero, then the number of QVECT boundary conditions
! nqvceid = number of QVECT boundary conditions, where emisivity/absorbtion id entered ! nqvceid = number of QVECT boundary conditions, where emisivity/absorbtion id entered
! istpnx = 1 if to stop at end of increment ! istpnx = 1 if to stop at end of increment
! imicro1 = 1 if micro1 interface is used ! imicro1 = 1 if micro1 interface is used
! iaxisymm = set to 1 if axisymmetric analysis ! iaxisymm = set to 1 if axisymmetric analysis
! jbreakglue = set to 1 if breaking glued option is used ! jbreakglue = set to 1 if breaking glued option is used
! iglstif = 1 if ddm and global stiffness matrix formed (sgi solver 6 or solver9) ! iglstif = 1 if ddm and global stiffness matrix formed (sgi solver 6 or solver9)
! jfastasm = 1 do fast assembly using SuperForm code ! jfastasm = 1 do fast assembly using SuperForm code
! iwear = set to 1 if wear model, set to 2 if wear model and coordinates updated ! iwear = set to 1 if wear model, set to 2 if wear model and coordinates updated
! iwearcf = set to 1 to store nodal coefficient of friction for wear calculation ! iwearcf = set to 1 to store nodal coefficient of friction for wear calculation
! imixmeth = set=1 then use nonlinear mixture material - allocate memory ! imixmeth = set=1 then use nonlinear mixture material - allocate memory
! ielcmadyn = flag for magnetodynamics ! ielcmadyn = flag for magnetodynamics
! 0 - electromagnetics using newmark beta ! 0 - electromagnetics using newmark beta
! 1 - transient magnetics using backward euler ! 1 - transient magnetics using backward euler
! idinout = flag to control if inside out elements should be deactivated ! idinout = flag to control if inside out elements should be deactivated
! igena_meth = 0 - generalized alpha parameters depend on whether or not contact ! igena_meth = 0 - generalized alpha parameters depend on whether or not contact
! is flagged (dynamic,7) ! is flagged (dynamic,7)
! 10 - generalized alpha parameters are optimized for a contact ! 10 - generalized alpha parameters are optimized for a contact
! analysis (dynamic,8) ! analysis (dynamic,8)
! 11 - generalized alpha parameters are optimized for an analysis ! 11 - generalized alpha parameters are optimized for an analysis
! without contact (dynamic,8) ! without contact (dynamic,8)
! !
!*********************************************************************** !***********************************************************************

56
trunk/creeps2008r1
View File

@ -1,28 +1,28 @@
! reformated to free format ! reformated to free format
!*********************************************************************** !***********************************************************************
! !
! File: creeps.cmn ! File: creeps.cmn
! !
! MSC.Marc include file ! MSC.Marc include file
! !
real(pReal) cptim,timinc,timinc_p,timinc_s,timincm,timinc_a,timinc_b,creept real(pReal) cptim,timinc,timinc_p,timinc_s,timincm,timinc_a,timinc_b,creept
integer(pInt) icptim,icfte,icfst,icfeq,icftm,icetem,mcreep,jcreep,icpa,icftmp,icfstr,& integer(pInt) icptim,icfte,icfst,icfeq,icftm,icetem,mcreep,jcreep,icpa,icftmp,icfstr,&
icfqcp,icfcpm,icrppr,icrcha,icpb,iicpmt,iicpa icfqcp,icfcpm,icrppr,icrcha,icpb,iicpmt,iicpa
real(pReal) time_beg_lcase,time_beg_inc,fractol,time_beg_pst real(pReal) time_beg_lcase,time_beg_inc,fractol,time_beg_pst
! !
integer num_creepsr,num_creepsi,num_creeps2r integer num_creepsr,num_creepsi,num_creeps2r
parameter(num_creepsr=40) parameter(num_creepsr=40)
parameter(num_creepsi=18) parameter(num_creepsi=18)
parameter(num_creeps2r=4) parameter(num_creeps2r=4)
common/marc_creeps/cptim,timinc,timinc_p,timinc_s,timincm,timinc_a,timinc_b,creept(33),icptim,icfte,icfst,& common/marc_creeps/cptim,timinc,timinc_p,timinc_s,timincm,timinc_a,timinc_b,creept(33),icptim,icfte,icfst,&
icfeq,icftm,icetem,mcreep,jcreep,icpa,icftmp,icfstr,icfqcp,icfcpm,icrppr,icrcha,icpb,iicpmt,iicpa icfeq,icftm,icetem,mcreep,jcreep,icpa,icftmp,icfstr,icfqcp,icfcpm,icrppr,icrcha,icpb,iicpmt,iicpa
common/marc_creeps2/time_beg_lcase,time_beg_inc,fractol,time_beg_pst common/marc_creeps2/time_beg_lcase,time_beg_inc,fractol,time_beg_pst
! !
! time_beg_lcase time at the beginning of the current load case ! time_beg_lcase time at the beginning of the current load case
! time_beg_inc time at the beginning of the current increment ! time_beg_inc time at the beginning of the current increment
! fractol fraction of loadcase or increment time when we ! fractol fraction of loadcase or increment time when we
! consider it to be finished ! consider it to be finished
! time_beg_pst time corresponding to first increment to be ! time_beg_pst time corresponding to first increment to be
! read in from thermal post file for auto step ! read in from thermal post file for auto step
! !
!*********************************************************************** !***********************************************************************

1258
trunk/lattice.f90
View File

File diff suppressed because it is too large Load Diff

206
trunk/math.f90
View File

@ -221,8 +221,8 @@
use prec, only: pReal, pInt use prec, only: pReal, pInt
implicit none implicit none
integer(pInt), intent(in) :: dimen integer(pInt), intent(in) :: dimen
integer(pInt) i,j,k,l integer(pInt) i,j,k,l
real(pReal), dimension(dimen,dimen,dimen,dimen) :: math_identity4th real(pReal), dimension(dimen,dimen,dimen,dimen) :: math_identity4th
forall (i=1:dimen,j=1:dimen,k=1:dimen,l=1:dimen) math_identity4th(i,j,k,l) = & forall (i=1:dimen,j=1:dimen,k=1:dimen,l=1:dimen) math_identity4th(i,j,k,l) = &
@ -232,47 +232,47 @@
END FUNCTION END FUNCTION
!************************************************************************** !**************************************************************************
! vector product a x b ! vector product a x b
!************************************************************************** !**************************************************************************
PURE FUNCTION math_vectorproduct(A,B) PURE FUNCTION math_vectorproduct(A,B)
use prec, only: pReal, pInt use prec, only: pReal, pInt
implicit none implicit none
real(pReal), dimension(3), intent(in) :: A,B real(pReal), dimension(3), intent(in) :: A,B
real(pReal), dimension(3) :: math_vectorproduct real(pReal), dimension(3) :: math_vectorproduct
math_vectorproduct(1) = A(2)*B(3)-A(3)*B(2) math_vectorproduct(1) = A(2)*B(3)-A(3)*B(2)
math_vectorproduct(2) = A(3)*B(1)-A(1)*B(3) math_vectorproduct(2) = A(3)*B(1)-A(1)*B(3)
math_vectorproduct(3) = A(1)*B(2)-A(2)*B(1) math_vectorproduct(3) = A(1)*B(2)-A(2)*B(1)
return return
END FUNCTION END FUNCTION
!************************************************************************** !**************************************************************************
! matrix multiplication 3x3 ! matrix multiplication 3x3
!************************************************************************** !**************************************************************************
PURE FUNCTION math_mul33x33(A,B)
PURE FUNCTION math_mul33x33(A,B)
use prec, only: pReal, pInt
implicit none use prec, only: pReal, pInt
implicit none
integer(pInt) i,j
real(pReal), dimension(3,3), intent(in) :: A,B integer(pInt) i,j
real(pReal), dimension(3,3) :: math_mul33x33 real(pReal), dimension(3,3), intent(in) :: A,B
real(pReal), dimension(3,3) :: math_mul33x33
forall (i=1:3,j=1:3) math_mul33x33(i,j) = &
A(i,1)*B(1,j) + A(i,2)*B(2,j) + A(i,3)*B(3,j) forall (i=1:3,j=1:3) math_mul33x33(i,j) = &
return A(i,1)*B(1,j) + A(i,2)*B(2,j) + A(i,3)*B(3,j)
return
END FUNCTION
END FUNCTION
!************************************************************************** !**************************************************************************
! matrix multiplication 6x6 ! matrix multiplication 6x6
!************************************************************************** !**************************************************************************
@ -282,8 +282,8 @@
implicit none implicit none
integer(pInt) i,j integer(pInt) i,j
real(pReal), dimension(6,6), intent(in) :: A,B real(pReal), dimension(6,6), intent(in) :: A,B
real(pReal), dimension(6,6) :: math_mul66x66 real(pReal), dimension(6,6) :: math_mul66x66
forall (i=1:6,j=1:6) math_mul66x66(i,j) = & forall (i=1:6,j=1:6) math_mul66x66(i,j) = &
A(i,1)*B(1,j) + A(i,2)*B(2,j) + A(i,3)*B(3,j) + & A(i,1)*B(1,j) + A(i,2)*B(2,j) + A(i,3)*B(3,j) + &
@ -292,25 +292,27 @@
END FUNCTION END FUNCTION
!**************************************************************************
! matrix multiplication 6x6 !**************************************************************************
!************************************************************************** ! matrix multiplication 6x6
PURE FUNCTION math_mul66x6(A,B) !**************************************************************************
PURE FUNCTION math_mul66x6(A,B)
use prec, only: pReal, pInt
implicit none use prec, only: pReal, pInt
implicit none
integer(pInt) i
real(pReal), dimension(6,6), intent(in) :: A integer(pInt) i
real(pReal), dimension(6), intent(in) :: B real(pReal), dimension(6,6), intent(in) :: A
real(pReal), dimension(6) :: math_mul66x6 real(pReal), dimension(6), intent(in) :: B
real(pReal), dimension(6) :: math_mul66x6
forall (i=1:6) math_mul66x6(i) = &
A(i,1)*B(1) + A(i,2)*B(2) + A(i,3)*B(3) + & forall (i=1:6) math_mul66x6(i) = &
A(i,4)*B(4) + A(i,5)*B(5) + A(i,6)*B(6) A(i,1)*B(1) + A(i,2)*B(2) + A(i,3)*B(3) + &
return A(i,4)*B(4) + A(i,5)*B(5) + A(i,6)*B(6)
return
END FUNCTION
END FUNCTION
!************************************************************************** !**************************************************************************
! matrix multiplication 9x9 ! matrix multiplication 9x9
@ -321,8 +323,10 @@
implicit none implicit none
integer(pInt) i,j integer(pInt) i,j
real(pReal), dimension(9,9), intent(in) :: A,B real(pReal), dimension(9,9), intent(in) :: A,B
real(pReal), dimension(9,9) :: math_mul99x99
real(pReal), dimension(9,9) :: math_mul99x99
forall (i=1:9,j=1:9) math_mul99x99(i,j) = & forall (i=1:9,j=1:9) math_mul99x99(i,j) = &
A(i,1)*B(1,j) + A(i,2)*B(2,j) + A(i,3)*B(3,j) + & A(i,1)*B(1,j) + A(i,2)*B(2,j) + A(i,3)*B(3,j) + &
@ -610,9 +614,9 @@
use prec, only: pReal,pInt use prec, only: pReal,pInt
implicit none implicit none
integer(pInt) i,j integer(pInt) i,j
real(pReal), dimension(6,6), intent(in) :: m real(pReal), dimension(6,6), intent(in) :: m
real(pReal), dimension(6,6) :: math_symmetric6x6 real(pReal), dimension(6,6) :: math_symmetric6x6
forall (i=1:6,j=1:6) math_symmetric6x6(i,j) = 1.0_pReal/2.0_pReal * & forall (i=1:6,j=1:6) math_symmetric6x6(i,j) = 1.0_pReal/2.0_pReal * &
(m(i,j) + m(j,i)) (m(i,j) + m(j,i))
@ -629,7 +633,7 @@
implicit none implicit none
real(pReal), dimension(3,3), intent(in) :: m real(pReal), dimension(3,3), intent(in) :: m
real(pReal) math_det3x3 real(pReal) math_det3x3
math_det3x3 = m(1,1)*(m(2,2)*m(3,3)-m(2,3)*m(3,2)) & math_det3x3 = m(1,1)*(m(2,2)*m(3,3)-m(2,3)*m(3,2)) &
-m(1,2)*(m(2,1)*m(3,3)-m(2,3)*m(3,1)) & -m(1,2)*(m(2,1)*m(3,3)-m(2,3)*m(3,1)) &
@ -807,10 +811,10 @@
use prec, only: pReal, pInt use prec, only: pReal, pInt
implicit none implicit none
real(pReal), dimension (3,3), intent(in) :: R real(pReal), dimension (3,3), intent(in) :: R
real(pReal), dimension(3) :: math_RtoEuler real(pReal), dimension(3) :: math_RtoEuler
real(pReal) sqhkl, squvw, sqhk, val real(pReal) sqhkl, squvw, sqhk, val
sqhkl=sqrt(R(1,3)*R(1,3)+R(2,3)*R(2,3)+R(3,3)*R(3,3)) sqhkl=sqrt(R(1,3)*R(1,3)+R(2,3)*R(2,3)+R(3,3)*R(3,3))
squvw=sqrt(R(1,1)*R(1,1)+R(2,1)*R(2,1)+R(3,1)*R(3,1)) squvw=sqrt(R(1,1)*R(1,1)+R(2,1)*R(2,1)+R(3,1)*R(3,1))
@ -862,9 +866,9 @@
use prec, only: pReal, pInt use prec, only: pReal, pInt
implicit none implicit none
real(pReal), dimension(3), intent(in) :: axis real(pReal), dimension(3), intent(in) :: axis
real(pReal), intent(in) :: omega real(pReal), intent(in) :: omega
real(pReal), dimension(3) :: axisNrm real(pReal), dimension(3) :: axisNrm
real(pReal), dimension(3,3) :: math_RodrigToR real(pReal), dimension(3,3) :: math_RodrigToR
real(pReal) s,c real(pReal) s,c
integer(pInt) i integer(pInt) i
@ -976,16 +980,11 @@
real(pReal) noise,scatter,cosScatter real(pReal) noise,scatter,cosScatter
integer(pInt) i integer(pInt) i
if (noise==0.0) then if (noise==0.0) then
math_sampleGaussOri = center math_sampleGaussOri = center
return return
endif endif
! Helming uses different distribution with Bessel functions ! Helming uses different distribution with Bessel functions
! therefore the gauss scatter width has to be scaled differently ! therefore the gauss scatter width has to be scaled differently
scatter = 0.95_pReal * noise scatter = 0.95_pReal * noise
@ -2034,28 +2033,27 @@ math_sampleFiberOri = math_RtoEuler(math_mul33x33(pRot,math_mul33x33(fRot,oRot))
END FUNCTION END FUNCTION
!************************************************************************** !**************************************************************************
! volume of tetrahedron given by four vertices ! volume of tetrahedron given by four vertices
!************************************************************************** !**************************************************************************
PURE FUNCTION math_volTetrahedron(v1,v2,v3,v4) PURE FUNCTION math_volTetrahedron(v1,v2,v3,v4)
use prec, only: pReal
implicit none
real(pReal) math_volTetrahedron
real(pReal), dimension (3), intent(in) :: v1,v2,v3,v4
real(pReal), dimension (3,3) :: m
m(:,1) = v1-v2
m(:,2) = v2-v3
m(:,3) = v3-v4
math_volTetrahedron = math_det3x3(m)/6.0_pReal
return
END FUNCTION
END MODULE math
use prec, only: pReal
implicit none
real(pReal) math_volTetrahedron
real(pReal), dimension (3), intent(in) :: v1,v2,v3,v4
real(pReal), dimension (3,3) :: m
m(:,1) = v1-v2
m(:,2) = v2-v3
m(:,3) = v3-v4
math_volTetrahedron = math_det3x3(m)/6.0_pReal
return
END FUNCTION
END MODULE math

166
trunk/mattex.mpie
View File

@ -1,83 +1,83 @@
<materials> <materials>
[TWIP steel FeMnC] [TWIP steel FeMnC]
lattice_structure 1 lattice_structure 1
Nslip 12 Nslip 12
Ntwin 0 Ntwin 0
## Elastic constants ## Elastic constants
# Unit in [Pa] # Unit in [Pa]
C11 183.9e9 C11 183.9e9
C12 101.9e9 C12 101.9e9
C44 115.4e9 C44 115.4e9
## Parameters for phenomenological modeling ## Parameters for phenomenological modeling
# Unit in [Pa] # Unit in [Pa]
s0_slip 85.0e6 s0_slip 85.0e6
gdot0_slip 0.001 gdot0_slip 0.001
n_slip 100.0 n_slip 100.0
h0 355.0e6 h0 355.0e6
s_sat 265.0e6 s_sat 265.0e6
w0 1.0 w0 1.0
# Self and latent hardening coefficients # Self and latent hardening coefficients
hardening_coefficients 1.0 1.4 hardening_coefficients 1.0 1.4
## Parameters for dislocation-based modeling ## Parameters for dislocation-based modeling
# Burgers vector [m] # Burgers vector [m]
burgers 2.56e-10 burgers 2.56e-10
# Activation energy for dislocation glide [J/K] (0.5*G*b^3) # Activation energy for dislocation glide [J/K] (0.5*G*b^3)
Qedge 5.5e-19 Qedge 5.5e-19
# Activation energy for self diffusion [J/K] (gamma-iron) # Activation energy for self diffusion [J/K] (gamma-iron)
Qsd 4.7e-19 Qsd 4.7e-19
# Vacancy diffusion coeffficent (gamma-iron) # Vacancy diffusion coeffficent (gamma-iron)
diff0 4.0e-5 diff0 4.0e-5
# Average grain size [m] # Average grain size [m]
grain_size 2.0e-5 grain_size 2.0e-5
# Dislocation interaction coefficients # Dislocation interaction coefficients
interaction_coefficients 1.0 2.2 3.0 1.6 3.8 4.5 interaction_coefficients 1.0 2.2 3.0 1.6 3.8 4.5
# Initial dislocation density [m]² # Initial dislocation density [m]²
rho0 6.0e12 rho0 6.0e12
# Passing stress adjustment # Passing stress adjustment
c1 0.1 c1 0.1
# Jump width adjustment # Jump width adjustment
c2 2.0 c2 2.0
# Activation volume adjustment # Activation volume adjustment
c3 1.0 c3 1.0
# Average slip distance adjustment for lock formation # Average slip distance adjustment for lock formation
c4 50.0 c4 50.0
# Average slip distance adjustment when grain boundaries # Average slip distance adjustment when grain boundaries
c5 1.0 c5 1.0
# Athermal recovery adjustment # Athermal recovery adjustment
c7 8.0 c7 8.0
# Thermal recovery adjustment (plays no role for me) # Thermal recovery adjustment (plays no role for me)
c8 1.0e10 c8 1.0e10
## Parameters for mechanical twinning ## Parameters for mechanical twinning
# Average twin thickness (stacks) [m] # Average twin thickness (stacks) [m]
stack_size 5.0e-8 stack_size 5.0e-8
# Total twin volume fraction saturation # Total twin volume fraction saturation
f_sat 1.0 f_sat 1.0
# Average slip distance adjustment when twin boundaries # Average slip distance adjustment when twin boundaries
c6 c6
# Scaling potential nucleation sites # Scaling potential nucleation sites
site_scaling 1.0e-6 site_scaling 1.0e-6
# Scaling the P-K force on the twinning dislocation # Scaling the P-K force on the twinning dislocation
q1 1.0 q1 1.0
# Scaling the resolved shear stress # Scaling the resolved shear stress
q2 1.0 q2 1.0
<textures> <textures>
[cube SX] [cube SX]
symmetry no /monoclinic /orthorhombic symmetry no /monoclinic /orthorhombic
Ngrains 10 /2 /4 Ngrains 10 /2 /4
#(gauss) phi1 0.0 phi 29.21 phi2 -26.57 scatter 0.0 fraction 1.0 #(gauss) phi1 0.0 phi 29.21 phi2 -26.57 scatter 0.0 fraction 1.0
#(gauss) phi1 0.0 phi 54.74 phi2 -45.0 scatter 0.0 fraction 0.1 #(gauss) phi1 0.0 phi 54.74 phi2 -45.0 scatter 0.0 fraction 0.1
#(gauss) phi1 0.0 phi 45.0 phi2 0.0 scatter 0.0 fraction 0.1 #(gauss) phi1 0.0 phi 45.0 phi2 0.0 scatter 0.0 fraction 0.1
#(gauss) phi1 0.0 phi 0.0 phi2 0.0 scatter 0.0 fraction 0.1 #(gauss) phi1 0.0 phi 0.0 phi2 0.0 scatter 0.0 fraction 0.1
#(gauss) phi1 0.0 phi 35.26 phi2 -45.0 scatter 0.0 fraction 0.1 #(gauss) phi1 0.0 phi 35.26 phi2 -45.0 scatter 0.0 fraction 0.1
#(gauss) phi1 0.0 phi 48.19 phi2 -26.57 scatter 0.0 fraction 0.1 #(gauss) phi1 0.0 phi 48.19 phi2 -26.57 scatter 0.0 fraction 0.1
#(gauss) phi1 0.0 phi 26.57 phi2 0.0 scatter 0.0 fraction 0.1 #(gauss) phi1 0.0 phi 26.57 phi2 0.0 scatter 0.0 fraction 0.1
#(gauss) phi1 0.0 phi 42.03 phi2 -33.69 scatter 0.0 fraction 0.1 #(gauss) phi1 0.0 phi 42.03 phi2 -33.69 scatter 0.0 fraction 0.1
#(gauss) phi1 0.0 phi 40.36 phi2 -11.31 scatter 0.0 fraction 0.1 #(gauss) phi1 0.0 phi 40.36 phi2 -11.31 scatter 0.0 fraction 0.1
#(gauss) phi1 0.0 phi 15.62 phi2 -26.57 scatter 0.0 fraction 0.1 #(gauss) phi1 0.0 phi 15.62 phi2 -26.57 scatter 0.0 fraction 0.1

28
trunk/mpie_cpfem_marc2007r1.f90
View File

@ -31,7 +31,7 @@
include "debug.f90" ! uses prec include "debug.f90" ! uses prec
include "math.f90" ! uses prec include "math.f90" ! uses prec
include "IO.f90" ! uses prec, debug, math include "IO.f90" ! uses prec, debug, math
include "FEsolving.f90" ! uses prec, IO include "FEsolving.f90" ! uses prec, IO
include "mesh.f90" ! uses prec, IO, math, FEsolving include "mesh.f90" ! uses prec, IO, math, FEsolving
include "lattice.f90" ! uses prec, math include "lattice.f90" ! uses prec, math
include "constitutive.f90" ! uses prec, IO, math, lattice, mesh, debug include "constitutive.f90" ! uses prec, IO, math, lattice, mesh, debug
@ -130,16 +130,11 @@
use math, only: invnrmMandel use math, only: invnrmMandel
implicit real(pReal) (a-h,o-z) implicit real(pReal) (a-h,o-z)
integer(pInt) computationMode integer(pInt) computationMode
dimension e(*),de(*),t(*),dt(*),g(*),d(ngens,*),s(*), n(2),coord(ncrd,*),disp(ndeg,*),matus(2),dispt(ndeg,*),ffn(itel,*),& dimension e(*),de(*),t(*),dt(*),g(*),d(ngens,*),s(*), n(2),coord(ncrd,*),disp(ndeg,*),matus(2),dispt(ndeg,*),ffn(itel,*),&
frotn(itel,*),strechn(itel),eigvn(itel,*),ffn1(itel,*),frotn1(itel,*),strechn1(itel),eigvn1(itel,*),kcus(2) frotn(itel,*),strechn(itel),eigvn(itel,*),ffn1(itel,*),frotn1(itel,*),strechn1(itel),eigvn1(itel,*),kcus(2)
! Marc common blocks are in fixed format so they have to be pasted in here ! Marc common blocks are in fixed format so they have to be pasted in here
! Beware of changes in newer Marc versions -- these are from 2005r3 ! Beware of changes in newer Marc versions -- these are from 2005r3
! concom is needed for inc, subinc, ncycle, lovl ! concom is needed for inc, subinc, ncycle, lovl
@ -168,14 +163,13 @@
cptim,timinc,timinc_p,timinc_s,timincm,timinc_a,timinc_b,creept(33),icptim,icfte,icfst,& cptim,timinc,timinc_p,timinc_s,timincm,timinc_a,timinc_b,creept(33),icptim,icfte,icfst,&
icfeq,icftm,icetem,mcreep,jcreep,icpa,icftmp,icfstr,icfqcp,icfcpm,icrppr,icrcha,icpb,iicpmt,iicpa icfeq,icftm,icetem,mcreep,jcreep,icpa,icftmp,icfstr,icfqcp,icfcpm,icrppr,icrcha,icpb,iicpmt,iicpa
if (inc == 0) then if (inc == 0) then
cycleCounter = 4 cycleCounter = 4
else else
if (theCycle > ncycle .or. theInc /= inc) cycleCounter = 0 ! reset counter for each cutback or new inc if (theCycle > ncycle .or. theInc /= inc) cycleCounter = 0 ! reset counter for each cutback or new inc
if (theCycle /= ncycle .or. theLovl /= lovl) then if (theCycle /= ncycle .or. theLovl /= lovl) then
cycleCounter = cycleCounter+1 ! ping pong cycleCounter = cycleCounter+1 ! ping pong
outdatedFFN1 = .false. outdatedFFN1 = .false.
endif endif
endif endif
if (cptim > theTime .or. theInc /= inc) then ! reached convergence if (cptim > theTime .or. theInc /= inc) then ! reached convergence
@ -195,12 +189,12 @@
if (computationMode == 2 .and. outdatedByNewInc) then if (computationMode == 2 .and. outdatedByNewInc) then
computationMode = 1 ! compute and age former results computationMode = 1 ! compute and age former results
outdatedByNewInc = .false. outdatedByNewInc = .false.
endif endif
if (computationMode == 2 .and. outdatedFFN1) then
computationMode = 4 ! return odd results to force new vyvle if (computationMode == 2 .and. outdatedFFN1) then
endif computationMode = 4 ! return odd results to force new vyvle
endif
theTime = cptim ! record current starting time theTime = cptim ! record current starting time
theInc = inc ! record current increment number theInc = inc ! record current increment number
theCycle = ncycle ! record current cycle count theCycle = ncycle ! record current cycle count
@ -208,7 +202,6 @@
call CPFEM_general(computationMode,ffn,ffn1,t(1),timinc,n(1),nn,s,mod(cycleCounter-4,4_pInt*ijaco)==0,d,ngens) call CPFEM_general(computationMode,ffn,ffn1,t(1),timinc,n(1),nn,s,mod(cycleCounter-4,4_pInt*ijaco)==0,d,ngens)
! Mandel: 11, 22, 33, SQRT(2)*12, SQRT(2)*23, SQRT(2)*13 ! Mandel: 11, 22, 33, SQRT(2)*12, SQRT(2)*23, SQRT(2)*13
! Marc: 11, 22, 33, 12, 23, 13 ! Marc: 11, 22, 33, 12, 23, 13
forall(i=1:ngens) d(1:ngens,i) = invnrmMandel(i)*d(1:ngens,i)*invnrmMandel(1:ngens) forall(i=1:ngens) d(1:ngens,i) = invnrmMandel(i)*d(1:ngens,i)*invnrmMandel(1:ngens)
@ -253,7 +246,8 @@
! assign result variable ! assign result variable
v=CPFEM_results(mod(jpltcd-1_pInt, CPFEM_Nresults+constitutive_maxNresults)+1_pInt,& v=CPFEM_results(mod(jpltcd-1_pInt, CPFEM_Nresults+constitutive_maxNresults)+1_pInt,&
(jpltcd-1_pInt)/(CPFEM_Nresults+constitutive_maxNresults)+1_pInt,& (jpltcd-1_pInt)/(CPFEM_Nresults+constitutive_maxNresults)+1_pInt,&
nn, mesh_FEasCP('elem', m)) nn, mesh_FEasCP('elem', m))
return return
END SUBROUTINE END SUBROUTINE
! !

25
trunk/mpie_cpfem_marc2007r1_sequential.f90
View File

@ -27,15 +27,15 @@
! - creeps: timinc ! - creeps: timinc
!******************************************************************** !********************************************************************
! !
include "prec.f90" ! uses nothing else include "prec.f90" ! uses nothing else
include "debug.f90" ! uses prec include "debug.f90" ! uses prec
include "math.f90" ! uses prec include "math.f90" ! uses prec
include "IO.f90" ! uses prec, debug, math include "IO.f90" ! uses prec, debug, math
include "FEsolving.f90" ! uses prec, IO include "FEsolving.f90" ! uses prec, IO
include "mesh.f90" ! uses prec, IO, math, FEsolving include "mesh.f90" ! uses prec, IO, math, FEsolving
include "lattice.f90" ! uses prec, math include "lattice.f90" ! uses prec, math
include "constitutive.f90" ! uses prec, IO, math, lattice, mesh, debug include "constitutive.f90" ! uses prec, IO, math, lattice, mesh, debug
! include "crystallite.f90" ! uses prec, debug, constitutive, mesh, math, IO ! include "crystallite.f90" ! uses prec, debug, constitutive, mesh, math, IO
include "CPFEM_sequential.f90" ! uses prec, math, mesh, constitutive, FEsolving, debug, lattice, IO, crystallite include "CPFEM_sequential.f90" ! uses prec, math, mesh, constitutive, FEsolving, debug, lattice, IO, crystallite
! !
@ -130,16 +130,11 @@
use math, only: invnrmMandel use math, only: invnrmMandel
implicit real(pReal) (a-h,o-z) implicit real(pReal) (a-h,o-z)
integer(pInt) computationMode integer(pInt) computationMode
dimension e(*),de(*),t(*),dt(*),g(*),d(ngens,*),s(*), n(2),coord(ncrd,*),disp(ndeg,*),matus(2),dispt(ndeg,*),ffn(itel,*),& dimension e(*),de(*),t(*),dt(*),g(*),d(ngens,*),s(*), n(2),coord(ncrd,*),disp(ndeg,*),matus(2),dispt(ndeg,*),ffn(itel,*),&
frotn(itel,*),strechn(itel),eigvn(itel,*),ffn1(itel,*),frotn1(itel,*),strechn1(itel),eigvn1(itel,*),kcus(2) frotn(itel,*),strechn(itel),eigvn(itel,*),ffn1(itel,*),frotn1(itel,*),strechn1(itel),eigvn1(itel,*),kcus(2)
! Marc common blocks are in fixed format so they have to be pasted in here ! Marc common blocks are in fixed format so they have to be pasted in here
! Beware of changes in newer Marc versions -- these are from 2005r3 ! Beware of changes in newer Marc versions -- these are from 2005r3
! concom is needed for inc, subinc, ncycle, lovl ! concom is needed for inc, subinc, ncycle, lovl
@ -168,7 +163,6 @@
cptim,timinc,timinc_p,timinc_s,timincm,timinc_a,timinc_b,creept(33),icptim,icfte,icfst,& cptim,timinc,timinc_p,timinc_s,timincm,timinc_a,timinc_b,creept(33),icptim,icfte,icfst,&
icfeq,icftm,icetem,mcreep,jcreep,icpa,icftmp,icfstr,icfqcp,icfcpm,icrppr,icrcha,icpb,iicpmt,iicpa icfeq,icftm,icetem,mcreep,jcreep,icpa,icftmp,icfstr,icfqcp,icfcpm,icrppr,icrcha,icpb,iicpmt,iicpa
if (inc == 0) then if (inc == 0) then
cycleCounter = 4 cycleCounter = 4
else else
@ -202,7 +196,6 @@
call CPFEM_general(computationMode,ffn,ffn1,t(1),timinc,n(1),nn,s,mod(cycleCounter-4,2_pInt*ijaco)==0,d,ngens) call CPFEM_general(computationMode,ffn,ffn1,t(1),timinc,n(1),nn,s,mod(cycleCounter-4,2_pInt*ijaco)==0,d,ngens)
! Mandel: 11, 22, 33, SQRT(2)*12, SQRT(2)*23, SQRT(2)*13 ! Mandel: 11, 22, 33, SQRT(2)*12, SQRT(2)*23, SQRT(2)*13
! Marc: 11, 22, 33, 12, 23, 13 ! Marc: 11, 22, 33, 12, 23, 13
forall(i=1:ngens) d(1:ngens,i) = invnrmMandel(i)*d(1:ngens,i)*invnrmMandel(1:ngens) forall(i=1:ngens) d(1:ngens,i) = invnrmMandel(i)*d(1:ngens,i)*invnrmMandel(1:ngens)

562
trunk/mpie_cpfem_marc2008r1.f90
View File

@ -1,292 +1,290 @@
!******************************************************************** !********************************************************************
! Material subroutine for MSC.Marc Version 0.1 ! Material subroutine for MSC.Marc Version 0.1
! !
! written by F. Roters, P. Eisenlohr, L. Hantcherli, W.A. Counts ! written by F. Roters, P. Eisenlohr, L. Hantcherli, W.A. Counts
! MPI fuer Eisenforschung, Duesseldorf ! MPI fuer Eisenforschung, Duesseldorf
! !
! last modified: 22.11.2008 ! last modified: 22.11.2008
!******************************************************************** !********************************************************************
! Usage: ! Usage:
! - choose material as hypela2 ! - choose material as hypela2
! - set statevariable 2 to index of material ! - set statevariable 2 to index of material
! - set statevariable 3 to index of texture ! - set statevariable 3 to index of texture
! - choose output of user variables if desired ! - choose output of user variables if desired
! - make sure the file "mattex.mpie" exists in the working ! - make sure the file "mattex.mpie" exists in the working
! directory ! directory
! - use nonsymmetric option for solver (e.g. direct ! - use nonsymmetric option for solver (e.g. direct
! profile or multifrontal sparse, the latter seems ! profile or multifrontal sparse, the latter seems
! to be faster!) ! to be faster!)
! - in case of ddm a symmetric solver has to be used ! - in case of ddm a symmetric solver has to be used
!******************************************************************** !********************************************************************
! Marc subroutines used: ! Marc subroutines used:
! - hypela2 ! - hypela2
! - plotv ! - plotv
! - quit ! - quit
!******************************************************************** !********************************************************************
! Marc common blocks included: ! Marc common blocks included:
! - concom: lovl, ncycle, inc, incsub ! - concom: lovl, ncycle, inc, incsub
! - creeps: timinc ! - creeps: timinc
!******************************************************************** !********************************************************************
! !
include "prec.f90" ! uses nothing else include "prec.f90" ! uses nothing else
include "debug.f90" ! uses prec include "debug.f90" ! uses prec
include "math.f90" ! uses prec include "math.f90" ! uses prec
include "IO.f90" ! uses prec, debug, math include "IO.f90" ! uses prec, debug, math
include "FEsolving.f90" ! uses prec, IO include "FEsolving.f90" ! uses prec, IO
include "mesh.f90" ! uses prec, IO, math, FEsolving include "mesh.f90" ! uses prec, IO, math, FEsolving
include "lattice.f90" ! uses prec, math include "lattice.f90" ! uses prec, math
include "constitutive.f90" ! uses prec, IO, math, lattice, mesh, debug include "constitutive.f90" ! uses prec, IO, math, lattice, mesh, debug
! include "crystallite.f90" ! uses prec, debug, constitutive, mesh, math, IO ! include "crystallite.f90" ! uses prec, debug, constitutive, mesh, math, IO
include "CPFEM.f90" ! uses prec, math, mesh, constitutive, FEsolving, debug, lattice, IO, crystallite include "CPFEM.f90" ! uses prec, math, mesh, constitutive, FEsolving, debug, lattice, IO, crystallite
! !
! !
SUBROUTINE hypela2(d,g,e,de,s,t,dt,ngens,n,nn,kcus,matus,ndi,& SUBROUTINE hypela2(d,g,e,de,s,t,dt,ngens,n,nn,kcus,matus,ndi,&
nshear,disp,dispt,coord,ffn,frotn,strechn,eigvn,ffn1,& nshear,disp,dispt,coord,ffn,frotn,strechn,eigvn,ffn1,&
frotn1,strechn1,eigvn1,ncrd,itel,ndeg,ndm,& frotn1,strechn1,eigvn1,ncrd,itel,ndeg,ndm,&
nnode,jtype,lclass,ifr,ifu) nnode,jtype,lclass,ifr,ifu)
!******************************************************************** !********************************************************************
! This is the Marc material routine ! This is the Marc material routine
!******************************************************************** !********************************************************************
! !
! ************* user subroutine for defining material behavior ************** ! ************* user subroutine for defining material behavior **************
! !
! !
! CAUTION : Due to calculation of the Deformation gradients, Stretch Tensors and ! CAUTION : Due to calculation of the Deformation gradients, Stretch Tensors and
! Rotation tensors at previous and current states, the analysis can be ! Rotation tensors at previous and current states, the analysis can be
! computationally expensive. Please use the user subroutine -> hypela ! computationally expensive. Please use the user subroutine -> hypela
! if these kinematic quantities are not needed in the constitutive model ! if these kinematic quantities are not needed in the constitutive model
! !
! !
! IMPORTANT NOTES : ! IMPORTANT NOTES :
! !
! (1) F,R,U are only available for continuum and membrane elements (not for ! (1) F,R,U are only available for continuum and membrane elements (not for
! shells and beams). ! shells and beams).
! !
! (2) For total Lagrangian formulation use the -> 'Elasticity,1' card(= ! (2) For total Lagrangian formulation use the -> 'Elasticity,1' card(=
! total Lagrange with large disp) in the parameter section of input deck. ! total Lagrange with large disp) in the parameter section of input deck.
! For updated Lagrangian formulation use the -> 'Plasticity,3' card(= ! For updated Lagrangian formulation use the -> 'Plasticity,3' card(=
! update+finite+large disp+constant d) in the parameter section of ! update+finite+large disp+constant d) in the parameter section of
! input deck. ! input deck.
! !
! !
! d stress strain law to be formed ! d stress strain law to be formed
! g change in stress due to temperature effects ! g change in stress due to temperature effects
! e total elastic strain ! e total elastic strain
! de increment of strain ! de increment of strain
! s stress - should be updated by user ! s stress - should be updated by user
! t state variables (comes in at t=n, must be updated ! t state variables (comes in at t=n, must be updated
! to have state variables at t=n+1) ! to have state variables at t=n+1)
! dt increment of state variables ! dt increment of state variables
! ngens size of stress - strain law ! ngens size of stress - strain law
! n element number ! n element number
! nn integration point number ! nn integration point number
! kcus(1) layer number ! kcus(1) layer number
! kcus(2) internal layer number ! kcus(2) internal layer number
! matus(1) user material identification number ! matus(1) user material identification number
! matus(2) internal material identification number ! matus(2) internal material identification number
! ndi number of direct components ! ndi number of direct components
! nshear number of shear components ! nshear number of shear components
! disp incremental displacements ! disp incremental displacements
! dispt displacements at t=n (at assembly, lovl=4) and ! dispt displacements at t=n (at assembly, lovl=4) and
! displacements at t=n+1 (at stress recovery, lovl=6) ! displacements at t=n+1 (at stress recovery, lovl=6)
! coord coordinates ! coord coordinates
! ncrd number of coordinates ! ncrd number of coordinates
! ndeg number of degrees of freedom ! ndeg number of degrees of freedom
! itel dimension of F and R, either 2 or 3 ! itel dimension of F and R, either 2 or 3
! nnode number of nodes per element ! nnode number of nodes per element
! jtype element type ! jtype element type
! lclass element class ! lclass element class
! ifr set to 1 if R has been calculated ! ifr set to 1 if R has been calculated
! ifu set to 1 if strech has been calculated ! ifu set to 1 if strech has been calculated
! !
! at t=n : ! at t=n :
! !
! ffn deformation gradient ! ffn deformation gradient
! frotn rotation tensor ! frotn rotation tensor
! strechn square of principal stretch ratios, lambda(i) ! strechn square of principal stretch ratios, lambda(i)
! eigvn(i,j) i principal direction components for j eigenvalues ! eigvn(i,j) i principal direction components for j eigenvalues
! !
! at t=n+1 : ! at t=n+1 :
! !
! ffn1 deformation gradient ! ffn1 deformation gradient
! frotn1 rotation tensor ! frotn1 rotation tensor
! strechn1 square of principal stretch ratios, lambda(i) ! strechn1 square of principal stretch ratios, lambda(i)
! eigvn1(i,j) i principal direction components for j eigenvalues ! eigvn1(i,j) i principal direction components for j eigenvalues
! !
! The following operation obtains U (stretch tensor) at t=n+1 : ! The following operation obtains U (stretch tensor) at t=n+1 :
! !
! call scla(un1,0.d0,itel,itel,1) ! call scla(un1,0.d0,itel,itel,1)
! do 3 k=1,3 ! do 3 k=1,3
! do 2 i=1,3 ! do 2 i=1,3
! do 1 j=1,3 ! do 1 j=1,3
! un1(i,j)=un1(i,j)+dsqrt(strechn1(k))*eigvn1(i,k)*eigvn1(j,k) ! un1(i,j)=un1(i,j)+dsqrt(strechn1(k))*eigvn1(i,k)*eigvn1(j,k)
!1 continue !1 continue
!2 continue !2 continue
!3 continue !3 continue
! !
use prec, only: pReal,pInt, ijaco use prec, only: pReal,pInt, ijaco
use FEsolving use FEsolving
use CPFEM, only: CPFEM_general use CPFEM, only: CPFEM_general
use math, only: invnrmMandel use math, only: invnrmMandel
! !
implicit none implicit none
!
! ** Start of generated type statements ** ! ** Start of generated type statements **
real(pReal) coord, d, de, disp, dispt, dt, e, eigvn, eigvn1, ffn, ffn1 real(pReal) coord, d, de, disp, dispt, dt, e, eigvn, eigvn1, ffn, ffn1
real(pReal) frotn, frotn1, g real(pReal) frotn, frotn1, g
integer(pInt) ifr, ifu, itel, jtype, kcus, lclass, matus, n, ncrd, ndeg integer(pInt) ifr, ifu, itel, jtype, kcus, lclass, matus, n, ncrd, ndeg
integer(pInt) ndi, ndm, ngens, nn, nnode, nshear integer(pInt) ndi, ndm, ngens, nn, nnode, nshear
real(pReal) s, strechn, strechn1, t real(pReal) s, strechn, strechn1, t
! ** End of generated type statements ** ! ** End of generated type statements **
!
dimension e(*),de(*),t(*),dt(*),g(*),d(ngens,*),s(*), n(2),coord(ncrd,*),disp(ndeg,*),matus(2),dispt(ndeg,*),ffn(itel,*),&
frotn(itel,*),strechn(itel),eigvn(itel,*),ffn1(itel,*),frotn1(itel,*),strechn1(itel),eigvn1(itel,*),kcus(2),&
lclass(2)
! !
dimension e(*),de(*),t(*),dt(*),g(*),d(ngens,*),s(*), n(2),coord(ncrd,*),disp(ndeg,*),matus(2),dispt(ndeg,*),ffn(itel,*),& ! Marc common blocks are in fixed format so they have to be reformated to free format (f90)
frotn(itel,*),strechn(itel),eigvn(itel,*),ffn1(itel,*),frotn1(itel,*),strechn1(itel),eigvn1(itel,*),kcus(2),& ! Beware of changes in newer Marc versions -- these are from 2005r3
lclass(2) ! concom is needed for inc, subinc, ncycle, lovl
! include "concom_f90"
! Marc common blocks are in fixed format so they have to be reformated to free format (f90) ! creeps is needed for timinc (time increment)
! Beware of changes in newer Marc versions -- these are from 2005r3 include "creeps_f90"
! concom is needed for inc, subinc, ncycle, lovl
include "concom_f90"
! creeps is needed for timinc (time increment)
include "creeps_f90"
!
integer(pInt) computationMode,i
! !
if (inc == 0) then integer(pInt) computationMode,i
cycleCounter = 4 !
else if (inc == 0) then
if (theCycle > ncycle .or. theInc /= inc) cycleCounter = 0 ! reset counter for each cutback or new inc cycleCounter = 4
if (theCycle /= ncycle .or. theLovl /= lovl) then else
cycleCounter = cycleCounter+1 ! ping pong if (theCycle > ncycle .or. theInc /= inc) cycleCounter = 0 ! reset counter for each cutback or new inc
outdatedFFN1 = .false. if (theCycle /= ncycle .or. theLovl /= lovl) then
endif cycleCounter = cycleCounter+1 ! ping pong
endif outdatedFFN1 = .false.
if (cptim > theTime .or. theInc /= inc) then ! reached convergence endif
lastIncConverged = .true. endif
outdatedByNewInc = .true. if (cptim > theTime .or. theInc /= inc) then ! reached convergence
endif lastIncConverged = .true.
outdatedByNewInc = .true.
if (mod(cycleCounter,2) /= 0) computationMode = 4 ! recycle endif
if (mod(cycleCounter,4) == 2) computationMode = 3 ! collect
if (mod(cycleCounter,4) == 0) computationMode = 2 ! compute if (mod(cycleCounter,2) /= 0) computationMode = 4 ! recycle
if (computationMode == 4 .and. ncycle == 0 .and. .not. lastIncConverged) & if (mod(cycleCounter,4) == 2) computationMode = 3 ! collect
computationMode = 6 ! recycle but restore known good consistent tangent if (mod(cycleCounter,4) == 0) computationMode = 2 ! compute
if (computationMode == 4 .and. lastIncConverged) then if (computationMode == 4 .and. ncycle == 0 .and. .not. lastIncConverged) &
computationMode = 5 ! recycle and record former consistent tangent computationMode = 6 ! recycle but restore known good consistent tangent
lastIncConverged = .false. if (computationMode == 4 .and. lastIncConverged) then
endif computationMode = 5 ! recycle and record former consistent tangent
if (computationMode == 2 .and. outdatedByNewInc) then lastIncConverged = .false.
computationMode = 1 ! compute and age former results endif
outdatedByNewInc = .false. if (computationMode == 2 .and. outdatedByNewInc) then
endif computationMode = 1 ! compute and age former results
if (computationMode == 2 .and. outdatedFFN1) then outdatedByNewInc = .false.
computationMode = 4 ! return odd results to force new vyvle endif
endif if (computationMode == 2 .and. outdatedFFN1) then
computationMode = 4 ! return odd results to force new vyvle
endif
theTime = cptim ! record current starting time
theInc = inc ! record current increment number theTime = cptim ! record current starting time
theCycle = ncycle ! record current cycle count theInc = inc ! record current increment number
theLovl = lovl ! record current lovl theCycle = ncycle ! record current cycle count
theLovl = lovl ! record current lovl
call CPFEM_general(computationMode,ffn,ffn1,t(1),timinc,n(1),nn,s,mod(cycleCounter-4,4_pInt*ijaco)==0,d,ngens)
call CPFEM_general(computationMode,ffn,ffn1,t(1),timinc,n(1),nn,s,mod(cycleCounter-4,4_pInt*ijaco)==0,d,ngens)
! Mandel: 11, 22, 33, SQRT(2)*12, SQRT(2)*23, SQRT(2)*13 ! Mandel: 11, 22, 33, SQRT(2)*12, SQRT(2)*23, SQRT(2)*13
! Marc: 11, 22, 33, 12, 23, 13 ! Marc: 11, 22, 33, 12, 23, 13
forall(i=1:ngens) d(1:ngens,i) = invnrmMandel(i)*d(1:ngens,i)*invnrmMandel(1:ngens) forall(i=1:ngens) d(1:ngens,i) = invnrmMandel(i)*d(1:ngens,i)*invnrmMandel(1:ngens)
s(1:ngens) = s(1:ngens)*invnrmMandel(1:ngens) s(1:ngens) = s(1:ngens)*invnrmMandel(1:ngens)
if(symmetricSolver) d(1:ngens,1:ngens) = 0.5_pReal*(d(1:ngens,1:ngens)+transpose(d(1:ngens,1:ngens))) if(symmetricSolver) d(1:ngens,1:ngens) = 0.5_pReal*(d(1:ngens,1:ngens)+transpose(d(1:ngens,1:ngens)))
return return
END SUBROUTINE END SUBROUTINE
! !
SUBROUTINE plotv(v,s,sp,etot,eplas,ecreep,t,m,nn,layer,ndi,nshear,jpltcd) SUBROUTINE plotv(v,s,sp,etot,eplas,ecreep,t,m,nn,layer,ndi,nshear,jpltcd)
!******************************************************************** !********************************************************************
! This routine sets user defined output variables for Marc ! This routine sets user defined output variables for Marc
!******************************************************************** !********************************************************************
! !
! select a variable contour plotting (user subroutine). ! select a variable contour plotting (user subroutine).
! !
! v variable ! v variable
! s (idss) stress array ! s (idss) stress array
! sp stresses in preferred direction ! sp stresses in preferred direction
! etot total strain (generalized) ! etot total strain (generalized)
! eplas total plastic strain ! eplas total plastic strain
! ecreep total creep strain ! ecreep total creep strain
! t current temperature ! t current temperature
! m element number ! m element number
! nn integration point number ! nn integration point number
! layer layer number ! layer layer number
! ndi (3) number of direct stress components ! ndi (3) number of direct stress components
! nshear (3) number of shear stress components ! nshear (3) number of shear stress components
! !
!******************************************************************** !********************************************************************
use prec, only: pReal,pInt use prec, only: pReal,pInt
use CPFEM, only: CPFEM_results, CPFEM_Nresults use CPFEM, only: CPFEM_results, CPFEM_Nresults
use constitutive, only: constitutive_maxNresults use constitutive, only: constitutive_maxNresults
use mesh, only: mesh_FEasCP use mesh, only: mesh_FEasCP
implicit none implicit none
! !
real(pReal) s(*),etot(*),eplas(*),ecreep(*),sp(*) real(pReal) s(*),etot(*),eplas(*),ecreep(*),sp(*)
real(pReal) v, t(*) real(pReal) v, t(*)
integer(pInt) m, nn, layer, ndi, nshear, jpltcd integer(pInt) m, nn, layer, ndi, nshear, jpltcd
! !
! assign result variable ! assign result variable
v=CPFEM_results(mod(jpltcd-1_pInt, CPFEM_Nresults+constitutive_maxNresults)+1_pInt,& v=CPFEM_results(mod(jpltcd-1_pInt, CPFEM_Nresults+constitutive_maxNresults)+1_pInt,&
(jpltcd-1_pInt)/(CPFEM_Nresults+constitutive_maxNresults)+1_pInt,& (jpltcd-1_pInt)/(CPFEM_Nresults+constitutive_maxNresults)+1_pInt,&
nn, mesh_FEasCP('elem', m)) nn, mesh_FEasCP('elem', m))
return return
END SUBROUTINE END SUBROUTINE
! !
! !
! subroutine utimestep(timestep,timestepold,icall,time,timeloadcase) ! subroutine utimestep(timestep,timestepold,icall,time,timeloadcase)
!******************************************************************** !********************************************************************
! This routine modifies the addaptive time step of Marc ! This routine modifies the addaptive time step of Marc
!******************************************************************** !********************************************************************
! use prec, only: pReal,pInt ! use prec, only: pReal,pInt
! use CPFEM, only : CPFEM_timefactor_max ! use CPFEM, only : CPFEM_timefactor_max
! implicit none ! implicit none
! !
! real(pReal) timestep, timestepold, time,timeloadcase ! real(pReal) timestep, timestepold, time,timeloadcase
! integer(pInt) icall ! integer(pInt) icall
! !
! user subroutine for modifying the time step in auto step ! user subroutine for modifying the time step in auto step
! !
! timestep : the current time step as suggested by marc ! timestep : the current time step as suggested by marc
! to be modified in this routine ! to be modified in this routine
! timestepold : the current time step before it was modified by marc ! timestepold : the current time step before it was modified by marc
! icall : =1 for setting the initial time step ! icall : =1 for setting the initial time step
! =2 if this routine is called during an increment ! =2 if this routine is called during an increment
! =3 if this routine is called at the beginning ! =3 if this routine is called at the beginning
! of the increment ! of the increment
! time : time at the start of the current increment ! time : time at the start of the current increment
! timeloadcase: time period of the current load case ! timeloadcase: time period of the current load case
! !
! it is in general not recommended to increase the time step ! it is in general not recommended to increase the time step
! during the increment. ! during the increment.
! this routine is called right after the time step has (possibly) ! this routine is called right after the time step has (possibly)
! been updated by marc. ! been updated by marc.
! !
! user coding ! user coding
! reduce timestep during increment in case mpie_timefactor is too large ! reduce timestep during increment in case mpie_timefactor is too large
! if(icall==2_pInt) then ! if(icall==2_pInt) then
! if(mpie_timefactor_max>1.25_pReal) then ! if(mpie_timefactor_max>1.25_pReal) then
! timestep=min(timestep,timestepold*0.8_pReal) ! timestep=min(timestep,timestepold*0.8_pReal)
! end if ! end if
! return ! return
! modify timestep at beginning of new increment ! modify timestep at beginning of new increment
! else if(icall==3_pInt) then ! else if(icall==3_pInt) then
! if(mpie_timefactor_max<=0.8_pReal) then ! if(mpie_timefactor_max<=0.8_pReal) then
! timestep=min(timestep,timestepold*1.25_pReal) ! timestep=min(timestep,timestepold*1.25_pReal)
! else if (mpie_timefactor_max<=1.0_pReal) then ! else if (mpie_timefactor_max<=1.0_pReal) then
! timestep=min(timestep,timestepold/mpie_timefactor_max) ! timestep=min(timestep,timestepold/mpie_timefactor_max)
! else if (mpie_timefactor_max<=1.25_pReal) then ! else if (mpie_timefactor_max<=1.25_pReal) then
! timestep=min(timestep,timestepold*1.01_pReal) ! timestep=min(timestep,timestepold*1.01_pReal)
! else ! else
! timestep=min(timestep,timestepold*0.8_pReal) ! timestep=min(timestep,timestepold*0.8_pReal)
! end if ! end if
! end if ! end if
! return ! return
! end ! end

79
trunk/mpie_cpfem_marc2008r1_sequential.f90
View File

@ -16,7 +16,8 @@
! - use nonsymmetric option for solver (e.g. direct ! - use nonsymmetric option for solver (e.g. direct
! profile or multifrontal sparse, the latter seems ! profile or multifrontal sparse, the latter seems
! to be faster!) ! to be faster!)
! - in case of ddm a symmetric solver has to be used ! - in case of ddm a symmetric solver has to be used
!******************************************************************** !********************************************************************
! Marc subroutines used: ! Marc subroutines used:
! - hypela2 ! - hypela2
@ -28,15 +29,15 @@
! - creeps: timinc ! - creeps: timinc
!******************************************************************** !********************************************************************
! !
include "prec.f90" ! uses nothing else include "prec.f90" ! uses nothing else
include "debug.f90" ! uses prec include "debug.f90" ! uses prec
include "math.f90" ! uses prec include "math.f90" ! uses prec
include "IO.f90" ! uses prec, debug, math include "IO.f90" ! uses prec, debug, math
include "FEsolving.f90" ! uses prec, IO include "FEsolving.f90" ! uses prec, IO
include "mesh.f90" ! uses prec, IO, math, FEsolving include "mesh.f90" ! uses prec, IO, math, FEsolving
include "lattice.f90" ! uses prec, math include "lattice.f90" ! uses prec, math
include "constitutive.f90" ! uses prec, IO, math, lattice, mesh, debug include "constitutive.f90" ! uses prec, IO, math, lattice, mesh, debug
! include "crystallite.f90" ! uses prec, debug, constitutive, mesh, math, IO ! include "crystallite.f90" ! uses prec, debug, constitutive, mesh, math, IO
include "CPFEM_sequential.f90" ! uses prec, math, mesh, constitutive, FEsolving, debug, lattice, IO, crystallite include "CPFEM_sequential.f90" ! uses prec, math, mesh, constitutive, FEsolving, debug, lattice, IO, crystallite
! !
@ -124,35 +125,34 @@
!2 continue !2 continue
!3 continue !3 continue
! !
use prec, only: pReal,pInt, ijaco
use prec, only: pReal,pInt, ijaco use FEsolving
use FEsolving use CPFEM, only: CPFEM_general
use CPFEM, only: CPFEM_general use math, only: invnrmMandel
use math, only: invnrmMandel !
! implicit none
implicit none !
! ! ** Start of generated type statements **
! ** Start of generated type statements ** real(pReal) coord, d, de, disp, dispt, dt, e, eigvn, eigvn1, ffn, ffn1
real(pReal) coord, d, de, disp, dispt, dt, e, eigvn, eigvn1, ffn, ffn1 real(pReal) frotn, frotn1, g
real(pReal) frotn, frotn1, g integer(pInt) ifr, ifu, itel, jtype, kcus, lclass, matus, n, ncrd, ndeg
integer(pInt) ifr, ifu, itel, jtype, kcus, lclass, matus, n, ncrd, ndeg integer(pInt) ndi, ndm, ngens, nn, nnode, nshear
integer(pInt) ndi, ndm, ngens, nn, nnode, nshear real(pReal) s, strechn, strechn1, t
real(pReal) s, strechn, strechn1, t ! ** End of generated type statements **
! ** End of generated type statements ** !
! dimension e(*),de(*),t(*),dt(*),g(*),d(ngens,*),s(*), n(2),coord(ncrd,*),disp(ndeg,*),matus(2),dispt(ndeg,*),ffn(itel,*),&
dimension e(*),de(*),t(*),dt(*),g(*),d(ngens,*),s(*), n(2),coord(ncrd,*),disp(ndeg,*),matus(2),dispt(ndeg,*),ffn(itel,*),& frotn(itel,*),strechn(itel),eigvn(itel,*),ffn1(itel,*),frotn1(itel,*),strechn1(itel),eigvn1(itel,*),kcus(2),&
frotn(itel,*),strechn(itel),eigvn(itel,*),ffn1(itel,*),frotn1(itel,*),strechn1(itel),eigvn1(itel,*),kcus(2),& lclass(2)
lclass(2) !
! ! Marc common blocks are in fixed format so they have to be reformated to free format (f90)
! Marc common blocks are in fixed format so they have to be reformated to free format (f90) ! Beware of changes in newer Marc versions -- these are from 2005r3
! Beware of changes in newer Marc versions -- these are from 2005r3 ! concom is needed for inc, subinc, ncycle, lovl
! concom is needed for inc, subinc, ncycle, lovl include "concom_f90"
include "concom_f90" ! creeps is needed for timinc (time increment)
! creeps is needed for timinc (time increment) include "creeps_f90"
include "creeps_f90" !
! integer(pInt) computationMode,i
integer(pInt) computationMode,i !
!
if (inc == 0) then if (inc == 0) then
cycleCounter = 4 cycleCounter = 4
else else
@ -186,7 +186,6 @@
call CPFEM_general(computationMode,ffn,ffn1,t(1),timinc,n(1),nn,s,mod(cycleCounter-4,2_pInt*ijaco)==0,d,ngens) call CPFEM_general(computationMode,ffn,ffn1,t(1),timinc,n(1),nn,s,mod(cycleCounter-4,2_pInt*ijaco)==0,d,ngens)
! Mandel: 11, 22, 33, SQRT(2)*12, SQRT(2)*23, SQRT(2)*13 ! Mandel: 11, 22, 33, SQRT(2)*12, SQRT(2)*23, SQRT(2)*13
! Marc: 11, 22, 33, 12, 23, 13 ! Marc: 11, 22, 33, 12, 23, 13
forall(i=1:ngens) d(1:ngens,i) = invnrmMandel(i)*d(1:ngens,i)*invnrmMandel(1:ngens) forall(i=1:ngens) d(1:ngens,i) = invnrmMandel(i)*d(1:ngens,i)*invnrmMandel(1:ngens)