From 31d3dae2790c1b8c24c191a2c5ebc49ed01d4498 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Mon, 23 Mar 2009 16:08:17 +0000 Subject: [PATCH] FE (Marc) interface is now general. Requires to execute "makeMe" in order to derive the specific interface routines for each Marc version (e.g., 2007r1 and 2008r1) from mpie_cpfem_marc.f90. Pls do not commit explicit mpie_cpfem_marcXXX.f90 anymore and ignore them for SVN control. --- trunk/makeMe | 26 ++ ...fem_marc2007r1.f90 => mpie_cpfem_marc.f90} | 17 +- trunk/mpie_cpfem_marc2008r1.f90 | 301 ------------------ 3 files changed, 37 insertions(+), 307 deletions(-) create mode 100755 trunk/makeMe rename trunk/{mpie_cpfem_marc2007r1.f90 => mpie_cpfem_marc.f90} (96%) delete mode 100644 trunk/mpie_cpfem_marc2008r1.f90 diff --git a/trunk/makeMe b/trunk/makeMe new file mode 100755 index 000000000..c002ce540 --- /dev/null +++ b/trunk/makeMe @@ -0,0 +1,26 @@ +#!/usr/bin/env python + +import os,sys + +architectures = { + 'marc': { + 'parent': 'mpie_cpfem_marc.f90', + 'versions' : ['%%MARCVERSION%%','2007r1','2008r1'], + }, + } + +for arch in architectures: + try: + parent = architectures[arch]['parent'] + parentFile = open(parent) + parentContent = parentFile.readlines() + parentFile.close() + except IOError: + print 'unable to open',parent + continue + + for version in architectures[arch]['versions'][1:]: + childFile = open(os.path.splitext(parent)[0]+version+os.path.splitext(parent)[1],'w') + for line in parentContent: + childFile.write(line.replace(architectures[arch]['versions'][0],version)) + childFile.close() diff --git a/trunk/mpie_cpfem_marc2007r1.f90 b/trunk/mpie_cpfem_marc.f90 similarity index 96% rename from trunk/mpie_cpfem_marc2007r1.f90 rename to trunk/mpie_cpfem_marc.f90 index 0813140d8..6ef61f1aa 100644 --- a/trunk/mpie_cpfem_marc2007r1.f90 +++ b/trunk/mpie_cpfem_marc.f90 @@ -1,10 +1,14 @@ !******************************************************************** -! Material subroutine for MSC.Marc Version 0.1 +! Material subroutine for MSC.Marc +! +! written by P. Eisenlohr, +! F. Roters, +! L. Hantcherli, +! W.A. Counts +! D.D. Tjahjanto ! -! written by F. Roters, P. Eisenlohr, L. Hantcherli, W.A. Counts ! MPI fuer Eisenforschung, Duesseldorf ! -! last modified: 22.11.2008 !******************************************************************** ! Usage: ! - choose material as hypela2 @@ -15,7 +19,8 @@ ! - 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 +! - in case of ddm (domain decomposition)a SYMMETRIC +! solver has to be used, i.e uncheck "non-symmetric" !******************************************************************** ! Marc subroutines used: ! - hypela2 @@ -148,8 +153,8 @@ ! Marc common blocks are in fixed format so they have to be reformated to free format (f90) ! Beware of changes in newer Marc versions - include "concom2007r1" ! concom is needed for inc, subinc, ncycle, lovl - include "creeps2007r1" ! creeps is needed for timinc (time increment) + include "concom%%MARCVERSION%%" ! concom is needed for inc, subinc, ncycle, lovl + include "creeps%%MARCVERSION%%" ! creeps is needed for timinc (time increment) integer(pInt) computationMode,i diff --git a/trunk/mpie_cpfem_marc2008r1.f90 b/trunk/mpie_cpfem_marc2008r1.f90 deleted file mode 100644 index 62f7e07c2..000000000 --- a/trunk/mpie_cpfem_marc2008r1.f90 +++ /dev/null @@ -1,301 +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: 22.11.2008 -!******************************************************************** -! Usage: -! - choose material as hypela2 -! - 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 -! 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 "material.f90" ! uses prec, math, IO, mesh - include "lattice.f90" ! uses prec, math - include "constitutive_phenomenological.f90" ! uses prec, math, IO, lattice, material, debug - include "constitutive_j2.f90" ! uses prec, math, IO, lattice, material, debug - include "constitutive_dislobased.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 - - 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 - use debug -! - 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) - -! Marc common blocks are in fixed format so they have to be reformated to free format (f90) -! 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 - 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 - 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 - debug_OuterLoopDistribution = 0_pInt - debug_cumLpTicks = 0 - debug_cumDotStateTicks = 0 - debug_cumLpCalls = 0_pInt - debug_cumDotStateCalls = 0_pInt - endif - endif - 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 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 - 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,4_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_maxSizePostResults - 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_maxSizePostResults)+1_pInt,& - (jpltcd-1_pInt)/(CPFEM_Nresults+constitutive_maxSizePostResults)+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 \ No newline at end of file