From a7135ed83b8774294ab53da6e9fa382e6a336854 Mon Sep 17 00:00:00 2001 From: Franz Roters Date: Thu, 5 Mar 2009 14:18:50 +0000 Subject: [PATCH] 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 --- trunk/CPFEM.f90 | 8 +- trunk/constitutive.f90 | 88 +++--- trunk/constitutive_phenomenological.f90 | 4 +- trunk/material.config | 8 +- trunk/material.f90 | 8 +- trunk/mesh.f90 | 2 +- trunk/mpie_cpfem_marc2007r1_sequential.f90 | 296 --------------------- trunk/mpie_cpfem_marc2008r1.f90 | 46 ++-- trunk/mpie_cpfem_marc2008r1_sequential.f90 | 286 -------------------- trunk/todo.txt | 1 + 10 files changed, 84 insertions(+), 663 deletions(-) delete mode 100644 trunk/mpie_cpfem_marc2007r1_sequential.f90 delete mode 100644 trunk/mpie_cpfem_marc2008r1_sequential.f90 diff --git a/trunk/CPFEM.f90 b/trunk/CPFEM.f90 index 70d70a132..1bb458b84 100644 --- a/trunk/CPFEM.f90 +++ b/trunk/CPFEM.f90 @@ -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 -!############################################################## - +!############################################################## \ No newline at end of file diff --git a/trunk/constitutive.f90 b/trunk/constitutive.f90 index 168b054b3..8e121fd6b 100644 --- a/trunk/constitutive.f90 +++ b/trunk/constitutive.f90 @@ -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 \ No newline at end of file diff --git a/trunk/constitutive_phenomenological.f90 b/trunk/constitutive_phenomenological.f90 index 79d1c3918..9db9e31f4 100644 --- a/trunk/constitutive_phenomenological.f90 +++ b/trunk/constitutive_phenomenological.f90 @@ -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 diff --git a/trunk/material.config b/trunk/material.config index 0c83019d3..24831b741 100644 --- a/trunk/material.config +++ b/trunk/material.config @@ -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 ##################### diff --git a/trunk/material.f90 b/trunk/material.f90 index 00bf09405..39012dd21 100644 --- a/trunk/material.f90 +++ b/trunk/material.f90 @@ -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 diff --git a/trunk/mesh.f90 b/trunk/mesh.f90 index d1d66092b..4bd782a0b 100644 --- a/trunk/mesh.f90 +++ b/trunk/mesh.f90 @@ -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 diff --git a/trunk/mpie_cpfem_marc2007r1_sequential.f90 b/trunk/mpie_cpfem_marc2007r1_sequential.f90 deleted file mode 100644 index fdfd0b867..000000000 --- a/trunk/mpie_cpfem_marc2007r1_sequential.f90 +++ /dev/null @@ -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 diff --git a/trunk/mpie_cpfem_marc2008r1.f90 b/trunk/mpie_cpfem_marc2008r1.f90 index d0c67025d..74eeeb033 100644 --- a/trunk/mpie_cpfem_marc2008r1.f90 +++ b/trunk/mpie_cpfem_marc2008r1.f90 @@ -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 ! diff --git a/trunk/mpie_cpfem_marc2008r1_sequential.f90 b/trunk/mpie_cpfem_marc2008r1_sequential.f90 deleted file mode 100644 index f8f5279cd..000000000 --- a/trunk/mpie_cpfem_marc2008r1_sequential.f90 +++ /dev/null @@ -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 diff --git a/trunk/todo.txt b/trunk/todo.txt index f000bca94..d7558e20c 100644 --- a/trunk/todo.txt +++ b/trunk/todo.txt @@ -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