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

removed some unused variables 

correct some spelling errors in material.config
adopted marc2008 interface routine
removed sequential interface routines as they are not needed any more
This commit is contained in:
Franz Roters 2009-03-05 14:18:50 +00:00
parent 9aa89bd9b5
commit a7135ed83b
10 changed files with 84 additions and 663 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
trunk/CPFEM.f90
View File

@ -171,7 +171,7 @@
real(pReal), dimension (3,3) :: ffn,ffn1,Kirchhoff_bar
real(pReal), dimension (3,3,3,3) :: H_bar, H_bar_sym
real(pReal), dimension(CPFEM_ngens) :: CPFEM_stress
real(pReal), dimension(CPFEM_ngens,CPFEM_ngens) :: CPFEM_jaco, odd_jaco
real(pReal), dimension(CPFEM_ngens,CPFEM_ngens) :: CPFEM_jaco
real(pReal) Temperature,CPFEM_dt,J_inverse
integer(pInt) CPFEM_mode ! 1: regular computation with aged results&
! 2: regular computation&
@ -291,7 +291,6 @@
real(pReal), intent(in) :: CPFEM_dt
integer(pInt) g,i,e
logical error
real(pReal) volfrac
real(pReal), dimension(3,3) :: U,R
@ -354,7 +353,7 @@
implicit none
character (len=128) msg
logical updateJaco, JacoOK, allConverged
logical updateJaco, allConverged
real(preal) dt
real(pReal), dimension(3,3) :: Fg_pert,Lp_pert, P_pert, Fp_pert, Fe_pert
real(pReal), dimension(6) :: Tstar_v
@ -784,5 +783,4 @@ Inner: do ! inner iteration: Lp
END SUBROUTINE
!
END MODULE
!##############################################################
!##############################################################

88
trunk/constitutive.f90
View File

@ -1,4 +1,4 @@
!************************************
!* Module: CONSTITUTIVE *
!************************************
@ -6,20 +6,20 @@
!* - constitutive equations *
!* - parameters definition *
!************************************
MODULE constitutive
!*** Include other modules ***
use prec
implicit none
type(p_vec), dimension(:,:,:), allocatable :: constitutive_state_old, & ! pointer array to old state variables of each grain
constitutive_state_new ! pointer array to new state variables of each grain
integer(pInt), dimension(:,:,:), allocatable :: constitutive_sizeDotState, & ! size of dotState array
constitutive_sizeState, & ! size of state array per grain
constitutive_sizePostResults ! size of postResults array per grain
integer(pInt) constitutive_maxSizeDotState,constitutive_maxSizeState,constitutive_maxSizePostResults
CONTAINS
!****************************************
!* - constitutive_init
@ -29,8 +29,8 @@ CONTAINS
!* - constitutive_dotState
!* - constitutive_postResults
!****************************************
subroutine constitutive_init()
!**************************************
!* Module initialization *
@ -41,21 +41,21 @@ subroutine constitutive_init()
use material
use constitutive_phenomenological
! use constitutive_dislocal
integer(pInt), parameter :: fileunit = 200
integer(pInt) e,i,g,myInstance
if(.not. IO_open_file(fileunit,material_configFile)) call IO_error (100) ! corrupt config file
call constitutive_phenomenological_init(fileunit) ! parse all phases of this constitution
! call constitutive_dislocal_init(fileunit)
close(fileunit)
allocate(constitutive_state_old(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems))
allocate(constitutive_state_new(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems))
allocate(constitutive_sizeDotState(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; constitutive_sizeDotState = 0_pInt
allocate(constitutive_sizeState(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; constitutive_sizeState = 0_pInt
allocate(constitutive_sizePostResults(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; constitutive_sizePostResults = 0_pInt
do e = 1,mesh_NcpElems ! loop over elements
do i = 1,FE_Nips(mesh_element(2,e)) ! loop over IPs
do g = 1,homogenization_Ngrains(mesh_element(3,e)) ! loop over grains
@ -79,12 +79,12 @@ subroutine constitutive_init()
constitutive_maxSizeDotState = maxval(constitutive_sizeDotState)
constitutive_maxSizeState = maxval(constitutive_sizeState)
constitutive_maxSizePostResults = maxval(constitutive_sizePostResults)
return
end subroutine
function constitutive_homogenizedC(ipc,ip,el)
!*********************************************************************
!* This function returns the homogenized elacticity matrix *
@ -99,22 +99,22 @@ function constitutive_homogenizedC(ipc,ip,el)
use constitutive_phenomenological
! use constitutive_dislocal
implicit none
!* Definition of variables
integer(pInt) ipc,ip,el
real(pReal), dimension(6,6) :: constitutive_homogenizedC
select case (phase_constitution(material_phase(ipc,ip,el)))
case (constitutive_phenomenological_label)
constitutive_homogenizedC = constitutive_phenomenological_homogenizedC(constitutive_state_new,ipc,ip,el)
case ('dislocal')
end select
return
end function
subroutine constitutive_microstructure(Temperature,ipc,ip,el)
!*********************************************************************
!* This function calculates from state needed variables *
@ -130,21 +130,21 @@ subroutine constitutive_microstructure(Temperature,ipc,ip,el)
use constitutive_phenomenological
! use constitutive_dislocal
implicit none
!* Definition of variables
integer(pInt) ipc,ip,el
real(pReal) Temperature
select case (phase_constitution(material_phase(ipc,ip,el)))
case (constitutive_phenomenological_label)
call constitutive_phenomenological_microstructure(Temperature,constitutive_state_new,ipc,ip,el)
case ('dislocal')
end select
end subroutine
subroutine constitutive_LpAndItsTangent(Lp,dLp_dTstar, Tstar_v,Temperature,ipc,ip,el)
!*********************************************************************
!* This subroutine contains the constitutive equation for *
@ -163,25 +163,25 @@ subroutine constitutive_LpAndItsTangent(Lp,dLp_dTstar, Tstar_v,Temperature,ipc,i
use constitutive_phenomenological
! use constitutive_dislocal
implicit none
!* Definition of variables
integer(pInt) ipc,ip,el
real(pReal) Temperature
real(pReal), dimension(6) :: Tstar_v
real(pReal), dimension(3,3) :: Lp
real(pReal), dimension(9,9) :: dLp_dTstar
select case (phase_constitution(material_phase(ipc,ip,el)))
case (constitutive_phenomenological_label)
call constitutive_phenomenological_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperature,constitutive_state_new,ipc,ip,el)
case ('dislocal')
end select
return
end subroutine
function constitutive_dotState(Tstar_v,Temperature,ipc,ip,el)
!*********************************************************************
!* This subroutine contains the constitutive equation for *
@ -200,23 +200,23 @@ function constitutive_dotState(Tstar_v,Temperature,ipc,ip,el)
use constitutive_phenomenological
! use constitutive_dislocal
implicit none
!* Definition of variables
integer(pInt) ipc,ip,el
real(pReal) Temperature
real(pReal), dimension(6) :: Tstar_v
real(pReal), dimension(constitutive_sizeDotState(ipc,ip,el)) :: constitutive_dotState
select case (phase_constitution(material_phase(ipc,ip,el)))
case (constitutive_phenomenological_label)
constitutive_dotState = constitutive_phenomenological_dotState(Tstar_v,Temperature,constitutive_state_new,ipc,ip,el)
case ('dislocal')
end select
return
end function
pure function constitutive_postResults(Tstar_v,Temperature,dt,ipc,ip,el)
!*********************************************************************
!* return array of constitutive results *
@ -232,22 +232,22 @@ pure function constitutive_postResults(Tstar_v,Temperature,dt,ipc,ip,el)
use constitutive_phenomenological
! use constitutive_dislocal
implicit none
!* Definition of variables
integer(pInt), intent(in) :: ipc,ip,el
real(pReal), intent(in) :: dt,Temperature
real(pReal), dimension(6), intent(in) :: Tstar_v
real(pReal), dimension(constitutive_sizePostResults(ipc,ip,el)) :: constitutive_postResults
constitutive_postResults = 0.0_pReal
select case (phase_constitution(material_phase(ipc,ip,el)))
case (constitutive_phenomenological_label)
constitutive_postResults = constitutive_phenomenological_postResults(Tstar_v,Temperature,dt,constitutive_state_new,ipc,ip,el)
case ('dislocal')
end select
return
end function
END MODULE
END MODULE

4
trunk/constitutive_phenomenological.f90
View File

@ -61,7 +61,7 @@ subroutine constitutive_phenomenological_init(file)
integer(pInt), intent(in) :: file
integer(pInt), parameter :: maxNchunks = 7
integer(pInt), dimension(1+2*maxNchunks) :: positions
integer(pInt) section, maxNinstance, i,j,k,l, output
integer(pInt) section, maxNinstance, i,j,k, output
character(len=64) tag
character(len=1024) line
@ -415,7 +415,7 @@ pure function constitutive_phenomenological_postResults(Tstar_v,Temperature,dt,s
real(pReal), dimension(6), intent(in) :: Tstar_v
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: state
integer(pInt) matID,o,i,c,n
real(pReal) tau_slip, active_rate
real(pReal) tau_slip
real(pReal), dimension(constitutive_phenomenological_sizePostResults(phase_constitutionInstance(material_phase(ipc,ip,el)))) :: &
constitutive_phenomenological_postResults

8
trunk/material.config
View File

@ -19,14 +19,14 @@ Ngrains 4
#####################
[Aluminum_CubeSX]
(constitutent) phase 1 texture 2 fraction 1.0
(constituent) phase 1 texture 2 fraction 1.0
[Copper_rCubeSX]
(constitutent) phase 2 texture 3 fraction 1.0
(constituent) phase 2 texture 3 fraction 1.0
[DPsteel]
(constitutent) phase 3 texture 1 fraction 0.8
(constitutent) phase 4 texture 1 fraction 0.2
(constituent) phase 3 texture 1 fraction 0.8
(constituent) phase 4 texture 1 fraction 0.2
#####################

8
trunk/material.f90
View File

@ -421,8 +421,8 @@ subroutine material_populateGrains()
real(pReal), dimension (:,:), allocatable :: orientationOfGrain
real(pReal), dimension (3) :: orientation
real(pReal), dimension (3,3) :: symOrientation
integer(pInt) t,n,e,i,g,j,m,homog,micro,sgn
integer(pInt) phaseID,textureID,dGrains,myNgrains,myNorientations,NorientationsOfConstituent, &
integer(pInt) t,e,i,g,j,m,homog,micro,sgn
integer(pInt) phaseID,textureID,dGrains,myNgrains,myNorientations, &
grain,constituentGrain,symExtension
real(pReal) extreme,rnd
@ -516,7 +516,7 @@ subroutine material_populateGrains()
texture_Gauss( 4,t,textureID))
enddo
constituentGrain = constituentGrain + int(myNorientations*texture_Gauss(5,t,textureID))
write (6,*) 'now at constitutent grain',constituentGrain
write (6,*) 'now at constituent grain',constituentGrain
enddo
do t = 1,texture_Nfiber(textureID) ! loop over fiber components
@ -528,7 +528,7 @@ subroutine material_populateGrains()
texture_Fiber( 5,t,textureID))
enddo
constituentGrain = constituentGrain + int(myNorientations*texture_fiber(6,t,textureID))
write (6,*) 'now at constitutent grain',constituentGrain
write (6,*) 'now at constituent grain',constituentGrain
enddo
write (6,*) 'looping',constituentGrain+1,myNorientations

2
trunk/mesh.f90
View File

@ -1454,7 +1454,7 @@ matchFace: do j = 1,FE_NfaceNodes(-neighbor,t) ! count over nodes on matc
use math
implicit none
integer(pInt) e,f,t,i,j,k,n
integer(pInt) e,f,t,i,j,n
integer(pInt), parameter :: Ntriangles = FE_NipFaceNodes-2 ! each interface is made up of this many triangles
real(pReal), dimension (3,FE_NipFaceNodes) :: nPos ! coordinates of nodes on IP face
real(pReal), dimension(3,Ntriangles,FE_NipFaceNodes) :: normal

296
trunk/mpie_cpfem_marc2007r1_sequential.f90
View File

@ -1,296 +0,0 @@
!********************************************************************
! Material subroutine for MSC.Marc Version 0.1
!
! written by F. Roters, P. Eisenlohr, L. Hantcherli, W.A. Counts
! MPI fuer Eisenforschung, Duesseldorf
!
! last modified: 27.11.2008
!********************************************************************
! Usage:
! - choose material as hypela2
! - set statevariable 2 to index of material
! - set statevariable 3 to index of texture
! - choose output of user variables if desired
! - make sure the file "mattex.mpie" exists in the working
! directory
! - use nonsymmetric option for solver (e.g. direct
! profile or multifrontal sparse, the latter seems
! to be faster!)
!********************************************************************
! Marc subroutines used:
! - hypela2
! - plotv
! - quit
!********************************************************************
! Marc common blocks included:
! - concom: lovl, ncycle, inc, incsub
! - creeps: timinc
!********************************************************************
!
include "prec.f90" ! uses nothing else
include "debug.f90" ! uses prec
include "math.f90" ! uses prec
include "IO.f90" ! uses prec, debug, math
include "FEsolving.f90" ! uses prec, IO
include "mesh.f90" ! uses prec, IO, math, FEsolving
include "lattice.f90" ! uses prec, math
include "constitutive.f90" ! uses prec, IO, math, lattice, mesh, debug
! include "crystallite.f90" ! uses prec, debug, constitutive, mesh, math, IO
include "CPFEM_sequential.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,&
nshear,disp,dispt,coord,ffn,frotn,strechn,eigvn,ffn1,&
frotn1,strechn1,eigvn1,ncrd,itel,ndeg,ndm,&
nnode,jtype,lclass,ifr,ifu)
!********************************************************************
! This is the Marc material routine
!********************************************************************
!
! ************* user subroutine for defining material behavior **************
!
!
! CAUTION : Due to calculation of the Deformation gradients, Stretch Tensors and
! Rotation tensors at previous and current states, the analysis can be
! computationally expensive. Please use the user subroutine -> hypela
! if these kinematic quantities are not needed in the constitutive model
!
!
! IMPORTANT NOTES :
!
! (1) F,R,U are only available for continuum and membrane elements (not for
! shells and beams).
!
! (2) For total Lagrangian formulation use the -> 'Elasticity,1' card(=
! total Lagrange with large disp) in the parameter section of input deck.
! For updated Lagrangian formulation use the -> 'Plasticity,3' card(=
! update+finite+large disp+constant d) in the parameter section of
! input deck.
!
!
! d stress strain law to be formed
! g change in stress due to temperature effects
! e total elastic strain
! de increment of strain
! s stress - should be updated by user
! t state variables (comes in at t=n, must be updated
! to have state variables at t=n+1)
! dt increment of state variables
! ngens size of stress - strain law
! n element number
! nn integration point number
! kcus(1) layer number
! kcus(2) internal layer number
! matus(1) user material identification number
! matus(2) internal material identification number
! ndi number of direct components
! nshear number of shear components
! disp incremental displacements
! dispt displacements at t=n (at assembly, lovl=4) and
! displacements at t=n+1 (at stress recovery, lovl=6)
! coord coordinates
! ncrd number of coordinates
! ndeg number of degrees of freedom
! itel dimension of F and R, either 2 or 3
! nnode number of nodes per element
! jtype element type
! lclass element class
! ifr set to 1 if R has been calculated
! ifu set to 1 if strech has been calculated
!
! at t=n :
!
! ffn deformation gradient
! frotn rotation tensor
! strechn square of principal stretch ratios, lambda(i)
! eigvn(i,j) i principal direction components for j eigenvalues
!
! at t=n+1 :
!
! ffn1 deformation gradient
! frotn1 rotation tensor
! strechn1 square of principal stretch ratios, lambda(i)
! eigvn1(i,j) i principal direction components for j eigenvalues
!
! The following operation obtains U (stretch tensor) at t=n+1 :
!
! call scla(un1,0.d0,itel,itel,1)
! do 3 k=1,3
! do 2 i=1,3
! do 1 j=1,3
! un1(i,j)=un1(i,j)+dsqrt(strechn1(k))*eigvn1(i,k)*eigvn1(j,k)
!1 continue
!2 continue
!3 continue
!
use prec, only: pReal,pInt, ijaco
use FEsolving
use CPFEM, only: CPFEM_general
use math, only: invnrmMandel
implicit real(pReal) (a-h,o-z)
integer(pInt) computationMode
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)
! 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
! concom is needed for inc, subinc, ncycle, lovl
! include 'concom'
common/marc_concom/ &
iacous, iasmbl, iautth, ibear, icompl, iconj, icreep, ideva(50), idyn, idynt,&
ielas, ielcma, ielect, iform, ifour, iharm, ihcps, iheat, iheatt, ihresp,&
ijoule, ilem, ilnmom, iloren, inc, incext, incsub, ipass, iplres, ipois,&
ipoist, irpflo, ismall, ismalt, isoil, ispect, ispnow, istore, iswep, ithcrp,&
itherm, iupblg, iupdat, jacflg, jel, jparks, largst, lfond, loadup, loaduq,&
lodcor, lovl, lsub, magnet, ncycle, newtnt, newton, noshr, linear, ivscpl,&
icrpim, iradrt, ipshft, itshr, iangin, iupmdr, iconjf, jincfl, jpermg, jhour,&
isolvr, jritz, jtable, jshell, jdoubl, jform, jcentr, imini, kautth, iautof,&
ibukty, iassum, icnstd, icnstt, kmakmas, imethvp,iradrte,iradrtp, iupdate,iupdatp,&
ncycnt, marmen ,idynme, ihavca, ispf, kmini, imixed, largtt, kdoela, iautofg,&
ipshftp,idntrc, ipore, jtablm, jtablc, isnecma,itrnspo,imsdif, jtrnspo,mcnear,&
imech, imecht, ielcmat, ielectt,magnett, imsdift,noplas, jtabls, jactch, jtablth,&
kgmsto ,jpzo, ifricsh, iremkin,iremfor, ishearp,jspf, machining, jlshell,icompsol,&
iupblgfo,jcondir,nstcrp, nactive,ipassref, nstspnt,ibeart,icheckmpc, noline, icuring,&
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
! creeps is needed for timinc (time increment)
! include 'creeps'
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
if (inc == 0) then
cycleCounter = 4
else
if (theCycle > ncycle .or. theInc /= inc) cycleCounter = 0 ! reset counter for each cutback or new inc
if (theCycle /= ncycle .or. theLovl /= lovl) then
cycleCounter = cycleCounter+1 ! ping pong
endif
endif
if (cptim > theTime .or. theInc /= inc) then ! reached convergence
lastIncConverged = .true.
outdatedByNewInc = .true.
endif
if (mod(cycleCounter,2) == 0) computationMode = 2 ! compute
if (mod(cycleCounter,2) /= 0) computationMode = 4 ! recycle
if (computationMode == 4 .and. ncycle == 0 .and. .not. lastIncConverged) &
computationMode = 6 ! recycle but restore known good consistent tangent
if (computationMode == 4 .and. lastIncConverged) then
computationMode = 5 ! recycle and record former consistent tangent
lastIncConverged = .false.
endif
if (computationMode == 2 .and. outdatedByNewInc) then
computationMode = 1 ! compute and age former results
outdatedByNewInc = .false.
endif
theTime = cptim ! record current starting time
theInc = inc ! record current increment number
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,2_pInt*ijaco)==0,d,ngens)
! Mandel: 11, 22, 33, SQRT(2)*12, SQRT(2)*23, SQRT(2)*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)
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)))
return
END SUBROUTINE
!
SUBROUTINE plotv(v,s,sp,etot,eplas,ecreep,t,m,nn,layer,ndi,nshear,jpltcd)
!********************************************************************
! This routine sets user defined output variables for Marc
!********************************************************************
!
! select a variable contour plotting (user subroutine).
!
! v variable
! s (idss) stress array
! sp stresses in preferred direction
! etot total strain (generalized)
! eplas total plastic strain
! ecreep total creep strain
! t current temperature
! m element number
! nn integration point number
! layer layer number
! ndi (3) number of direct stress components
! nshear (3) number of shear stress components
!
!********************************************************************
use prec, only: pReal,pInt
use CPFEM, only: CPFEM_results, CPFEM_Nresults
use constitutive, only: constitutive_maxNresults
use mesh, only: mesh_FEasCP
implicit none
!
real(pReal) s(*),etot(*),eplas(*),ecreep(*),sp(*)
real(pReal) v, t(*)
integer(pInt) m, nn, layer, ndi, nshear, jpltcd
!
! assign result variable
v=CPFEM_results(mod(jpltcd-1_pInt, CPFEM_Nresults+constitutive_maxNresults)+1_pInt,&
(jpltcd-1_pInt)/(CPFEM_Nresults+constitutive_maxNresults)+1_pInt,&
nn, mesh_FEasCP('elem', m))
return
END SUBROUTINE
!
!
! subroutine utimestep(timestep,timestepold,icall,time,timeloadcase)
!********************************************************************
! This routine modifies the addaptive time step of Marc
!********************************************************************
! use prec, only: pReal,pInt
! use CPFEM, only : CPFEM_timefactor_max
! implicit none
!
! real(pReal) timestep, timestepold, time,timeloadcase
! integer(pInt) icall
!
! user subroutine for modifying the time step in auto step
!
! timestep : the current time step as suggested by marc
! to be modified in this routine
! timestepold : the current time step before it was modified by marc
! icall : =1 for setting the initial time step
! =2 if this routine is called during an increment
! =3 if this routine is called at the beginning
! of the increment
! time : time at the start of the current increment
! timeloadcase: time period of the current load case
!
! it is in general not recommended to increase the time step
! during the increment.
! this routine is called right after the time step has (possibly)
! been updated by marc.
!
! user coding
! reduce timestep during increment in case mpie_timefactor is too large
! if(icall==2_pInt) then
! if(mpie_timefactor_max>1.25_pReal) then
! timestep=min(timestep,timestepold*0.8_pReal)
! end if
! return
! modify timestep at beginning of new increment
! else if(icall==3_pInt) then
! if(mpie_timefactor_max<=0.8_pReal) then
! timestep=min(timestep,timestepold*1.25_pReal)
! else if (mpie_timefactor_max<=1.0_pReal) then
! timestep=min(timestep,timestepold/mpie_timefactor_max)
! else if (mpie_timefactor_max<=1.25_pReal) then
! timestep=min(timestep,timestepold*1.01_pReal)
! else
! timestep=min(timestep,timestepold*0.8_pReal)
! end if
! end if
! return
! end

46
trunk/mpie_cpfem_marc2008r1.f90
View File

@ -8,10 +8,9 @@
!********************************************************************
! Usage:
! - choose material as hypela2
! - set statevariable 2 to index of material
! - set statevariable 3 to index of texture
! - choose output of user variables if desired
! - make sure the file "mattex.mpie" exists in the working
! - set statevariable 2 to index of homogenization
! - set statevariable 3 to index of microstructure
! - make sure the file "material.config" exists in the working
! directory
! - use nonsymmetric option for solver (e.g. direct
! profile or multifrontal sparse, the latter seems
@ -35,15 +34,16 @@
include "FEsolving.f90" ! uses prec, IO
include "mesh.f90" ! uses prec, IO, math, FEsolving
include "lattice.f90" ! uses prec, math
include "material.f90" ! uses prec, math, IO, mesh
include "constitutive_phenomenological.f90" ! uses prec, math, IO, lattice, material, debug
include "constitutive.f90" ! uses prec, IO, math, lattice, mesh, debug
include "CPFEM.f90" ! uses prec, math, mesh, constitutive, FEsolving, debug, lattice, IO, crystallite
logical, parameter :: parallelExecution = .true.
SUBROUTINE hypela2(d,g,e,de,s,t,dt,ngens,n,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)
!********************************************************************
! This is the Marc material routine
!********************************************************************
@ -143,16 +143,17 @@
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)
!
! 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
! concom is needed for inc, subinc, ncycle, lovl
include "concom2008r1"
! creeps is needed for timinc (time increment)
include "creeps2008r1"
!
! Beware of changes in newer Marc versions
include "concom2008r1" ! concom is needed for inc, subinc, ncycle, lovl
include "creeps2008r1" ! creeps is needed for timinc (time increment)
integer(pInt) computationMode,i
!
! write(6,'(3(3(f10.3,x),/))') ffn1(:,1),ffn1(:,2),ffn1(:,3)
if (inc == 0) then
cycleCounter = 4
else
@ -160,6 +161,7 @@
if (theCycle /= ncycle .or. theLovl /= lovl) then
cycleCounter = cycleCounter+1 ! ping pong
outdatedFFN1 = .false.
write (6,*) n(1),nn,'cycleCounter',cycleCounter
call debug_info() ! output of debugging/performance statistics of former
debug_cutbackDistribution = 0_pInt ! initialize debugging data
debug_InnerLoopDistribution = 0_pInt
@ -169,11 +171,12 @@
if (cptim > theTime .or. theInc /= inc) then ! reached convergence
lastIncConverged = .true.
outdatedByNewInc = .true.
write (6,*) n(1),nn,'lastIncConverged + outdated'
endif
if (mod(cycleCounter,2) /= 0) computationMode = 4 ! recycle
if (mod(cycleCounter,4) == 2) computationMode = 3 ! collect
if (mod(cycleCounter,4) == 0) computationMode = 2 ! compute
if (mod(cycleCounter,2) /= 0) computationMode = 4 ! recycle in odd cycles
if (mod(cycleCounter,4) == 2) computationMode = 3 ! collect in 2,6,10,...
if (mod(cycleCounter,4) == 0) computationMode = 2 ! compute in 0,4,8,...
if (computationMode == 4 .and. ncycle == 0 .and. .not. lastIncConverged) &
computationMode = 6 ! recycle but restore known good consistent tangent
if (computationMode == 4 .and. lastIncConverged) then
@ -197,6 +200,7 @@
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)
if(symmetricSolver) d(1:ngens,1:ngens) = 0.5_pReal*(d(1:ngens,1:ngens)+transpose(d(1:ngens,1:ngens)))
return
END SUBROUTINE
@ -225,7 +229,7 @@
!********************************************************************
use prec, only: pReal,pInt
use CPFEM, only: CPFEM_results, CPFEM_Nresults
use constitutive, only: constitutive_maxNresults
use constitutive, only: constitutive_maxSizePostResults
use mesh, only: mesh_FEasCP
implicit none
!
@ -234,9 +238,9 @@
integer(pInt) m, nn, layer, ndi, nshear, jpltcd
!
! assign result variable
v=CPFEM_results(mod(jpltcd-1_pInt, CPFEM_Nresults+constitutive_maxSizePostResults)+1_pInt,&
(jpltcd-1_pInt)/(CPFEM_Nresults+constitutive_maxSizePostResults)+1_pInt,&
nn, mesh_FEasCP('elem', m))
v = CPFEM_results(mod(jpltcd-1_pInt, CPFEM_Nresults+constitutive_maxSizePostResults)+1_pInt,&
(jpltcd-1_pInt)/(CPFEM_Nresults+constitutive_maxSizePostResults)+1_pInt,&
nn, mesh_FEasCP('elem', m))
return
END SUBROUTINE
!

286
trunk/mpie_cpfem_marc2008r1_sequential.f90
View File

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

1
trunk/todo.txt
View File

@ -2,3 +2,4 @@ Things to be implemented into the code
# adopt CPFEM_GIA8.f90 to new scheme, use "select case" to switch between homogenization schemes in CPFEM.f90
# make OpenMP parallelization work again
# adopt constitutive_dislo