[cig-commits] r5372 -
short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d
willic3 at geodynamics.org
willic3 at geodynamics.org
Tue Nov 28 19:35:55 PST 2006
Author: willic3
Date: 2006-11-28 19:35:54 -0800 (Tue, 28 Nov 2006)
New Revision: 5372
Modified:
short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/mat_8.f
Log:
Changed old ESF-based Maxwell model to correspond to a new material
model.
Modified: short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/mat_8.f
===================================================================
--- short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/mat_8.f 2006-11-29 03:33:58 UTC (rev 5371)
+++ short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/mat_8.f 2006-11-29 03:35:54 UTC (rev 5372)
@@ -30,13 +30,14 @@
c... Program segment to define material model 8:
c
c Model number: 8
-c Model name: ??
-c Number material properties: ??
-c Number state variables: ??
-c Tangent matrix varies with state: ??
-c Material properties: ??
-c ??
-c ??
+c Model name: IsotropicLinearMaxwellViscoelasticESF
+c Number material properties: 4
+c Number state variables: 3
+c Tangent matrix varies with state: True
+c Material properties: Density
+c Young's modulus
+c Poisson's ratio
+c Viscosity
c
subroutine mat_prt_8(prop,nprop,matnum,idout,idsk,kw,kp,
& ierr,errstrng)
@@ -45,7 +46,6 @@
c
c Error codes:
c 4: Write error
-c 101: Attempt to use undefined material model
c
include "implicit.inc"
c
@@ -61,20 +61,54 @@
c
c... local constants
c
- character labelp(3)*15,modelname*24
+ character labelp(4)*15,modelname*41
data labelp/"Density",
& "Young's modulus",
- & "Poisson's ratio"/
- parameter(modelname="???????")
+ & "Poisson's ratio",
+ & "Viscosity"/
+ parameter(modelname="Isotropic Linear Maxwell Viscoelastic ESF")
integer mattype
parameter(mattype=8)
c
-c... return error code, as this material is not yet defined
+c... local variables
c
- ierr=101
- errstrng="mat_prt_8"
+ integer i
c
+cdebug write(6,*) "Hello from mat_prt_8_f!"
+c
+c... output plot results
+c
+ if(idsk.eq.ione) then
+ write(kp,"(3i7)",err=10) matnum,mattype,nprop
+ write(kp,"(1pe15.8,20(2x,1pe15.8))",err=10) (prop(i),i=1,nprop)
+ else if(idsk.eq.itwo) then
+ write(kp,err=10) matnum,mattype,nprop
+ write(kp,err=10) prop
+ end if
+c
+c... output ascii results, if desired
+c
+ if(idout.gt.izero) then
+ write(kw,700,err=10) matnum,modelname,nprop
+ do i=1,nprop
+ write(kw,710,err=10) labelp(i),prop(i)
+ end do
+ end if
+c
return
+c
+c... error writing to output file
+c
+ 10 continue
+ ierr=4
+ errstrng="mat_prt_8"
+ return
+c
+ 700 format(/,5x,"Material number: ",i7,/,5x,
+ & "Material type: ",a37,/,5x,
+ & "Number of properties: ",i7,/)
+ 710 format(15x,a15,3x,1pe15.8)
+c
end
c
c
@@ -100,10 +134,30 @@
character errstrng*(*)
double precision dmat(nddmat),prop(nprop)
c
-c... return error code, as this material is not yet defined
+c... local variables
c
- ierr=101
- errstrng="elas_mat_8"
+ integer i,j
+ double precision e,pr,pr1,pr2,pr3,fac,dd,od,ss
+c
+cdebug write(6,*) "Hello from elas_mat_8_f!"
+c
+ call fill(dmat,zero,nddmat)
+ e=prop(2)
+ pr=prop(3)
+ pr1=one-pr
+ pr2=one+pr
+ pr3=one-two*pr
+ fac=e/(pr2*pr3)
+ dd=pr1*fac
+ od=pr*fac
+ ss=half*pr3*fac
+ do i=1,3
+ dmat(iddmat(i,i))=dd
+ dmat(iddmat(i+3,i+3))=ss
+ do j=i+1,3
+ dmat(iddmat(i,j))=od
+ end do
+ end do
return
end
c
@@ -111,10 +165,20 @@
subroutine elas_strs_8(prop,nprop,state,state0,ee,scur,dmat,tmax,
& nstate,nstate0,ierr,errstrng)
c
-c... subroutine to compute stresses for the elastic solution.
+c... subroutine to compute stresses for the elastic solution. For this
+c material, there are 3 sets of state variables: total stress,
+c total strain, and viscous strain. The Maxwell time is computed,
+c even though this is the elastic solution, as an aid in determining
+c the proper time step size for the next step.
c The current total strain is contained in ee and the computed
-c total strerss should be copied to scur.
+c total stress should be copied to scur.
c
+c state(1:6) = Cauchy stress
+c state(7:12) = linear strain
+c state(13:18) = viscous strain
+c
+c The state0 array contains initial stresses.
+c
include "implicit.inc"
c
c... parameter definitions
@@ -131,10 +195,24 @@
double precision dmat(nddmat),tmax
character errstrng*(*)
c
-c... return error code, as this material is not yet defined
+c... local variables
c
- ierr=101
- errstrng="elas_strs_8"
+ double precision e,pr,vis,rmu
+c
+cdebug write(6,*) "Hello from elas_strs_8_f!"
+c
+ call dcopy(nstr,ee,ione,state(7),ione)
+ call dcopy(nstr,state0,ione,state,ione)
+ call dspmv("u",nstr,one,dmat,state(7),ione,one,state,ione)
+ call dcopy(nstr,state,ione,scur,ione)
+c
+c... compute Maxwell time for current stress state
+c
+ e=prop(2)
+ pr=prop(3)
+ vis=prop(4)
+ rmu=half*e/(one+pr)
+ tmax=vis/rmu
return
end
c
@@ -149,6 +227,9 @@
c computed. Thus, for some time-dependent materials, the material
c matrix will only be an approximation of the material matrix for
c the current iteration.
+c As this is a linear viscoelastic material, the tangent matrix is
+c actually independent of the current stresses and strains, so they
+c are not used.
c Note that only the upper triangle is used (or available), as
c dmat is assumed to always be symmetric.
c
@@ -174,16 +255,35 @@
include "rgiter_dim.inc"
include "ntimdat_dim.inc"
c
+c... local variables
+c
+ double precision e,pr,vis,f1,f2,rmu
+c
c... included variable definitions
c
include "rtimdat_def.inc"
include "rgiter_def.inc"
include "ntimdat_def.inc"
c
-c... return error code, as this material is not yet defined
+cdebug write(6,*) "Hello from td_matinit_8_f!"
c
- ierr=101
- errstrng="td_matinit_8"
+ call fill(dmat,zero,nddmat)
+ e=prop(2)
+ pr=prop(3)
+ vis=prop(4)
+ rmu=half*e/(one+pr)
+ tmax=vis/rmu
+ f1=third*e/(one-two*pr)
+ f2=third/((one+pr)/e+half*alfap*deltp/vis)
+ dmat(iddmat(1,1))=f1+two*f2
+ dmat(iddmat(2,2))=dmat(iddmat(1,1))
+ dmat(iddmat(3,3))=dmat(iddmat(1,1))
+ dmat(iddmat(1,2))=f1-f2
+ dmat(iddmat(1,3))=dmat(iddmat(1,2))
+ dmat(iddmat(2,3))=dmat(iddmat(1,2))
+ dmat(iddmat(4,4))=half*three*f2
+ dmat(iddmat(5,5))=dmat(iddmat(4,4))
+ dmat(iddmat(6,6))=dmat(iddmat(4,4))
return
end
c
@@ -192,10 +292,8 @@
& rtimdat,rgiter,ntimdat,iddmat,tmax,nstate,nstate0,nprop,matchg,
& ierr,errstrng)
c
-c... subroutine to compute the current stress for the time-dependent
-c solution.
-c Note that only the upper triangle is used (or available), as
-c dmat is assumed to always be symmetric.
+c... subroutine to compute the current values for stress, total strain,
+c and viscous strain increment.
c
include "implicit.inc"
c
@@ -220,17 +318,68 @@
include "rgiter_dim.inc"
include "ntimdat_dim.inc"
c
+c... local constants
+c
+ double precision diag(6)
+ data diag/one,one,one,zero,zero,zero/
+ double precision eng(6)
+ data eng/one,one,one,half,half,half/
+c
+c... local variables
+c
+ integer i
+ double precision e,pr,vis,rmu,ae,tf
+ double precision emeantdt,smeant,smeantdt,smean0
+ double precision gamtau,f1,f2,f3,f4
+ double precision epptdt,sdevt,sdev0,sdevtdt,sdevtau
+c
+cdebug integer idb
+c
c... included variable definitions
c
include "rtimdat_def.inc"
include "rgiter_def.inc"
include "ntimdat_def.inc"
c
-c... return error code, as this material is not yet defined
+cdebug write(6,*) "Hello from td_strs_8!"
c
- ierr=101
- errstrng="td_strs_8"
+ ierr=0
c
+c... define material properties
+c
+ e=prop(2)
+ pr=prop(3)
+ vis=prop(4)
+ rmu=half*e/(one+pr)
+ ae=(one+pr)/e
+ tf=deltp*(one-alfap)
+ gamtau=half/vis
+ f1=one/(ae+alfap*deltp*gamtau)
+ f2=gamtau*tf
+ f3=e/(one-two*pr)
+ f4=deltp*gamtau
+ tmax=vis/rmu
+c
+c... define stress and strain quantities
+c
+ emeantdt=third*(ee(1)+ee(2)+ee(3))
+ smeantdt=emeantdt*f3
+ smeant=third*(state(1)+state(2)+state(3))
+ smean0=third*(state0(1)+state0(2)+state0(3))
+c
+c... compute new stresses and store stress and strain values in dstate
+ do i=1,nstr
+ epptdt=eng(i)*ee(i)-diag(i)*emeantdt-state(i+12)
+ sdevt=state(i)-diag(i)*smeant
+ sdev0=state0(i)-diag(i)*smean0
+ sdevtdt=f1*(epptdt-f2*sdevt+ae*sdev0)
+ sdevtau=(one-alfap)*sdevt+alfap*sdevtdt
+ dstate(i)=sdevtdt+diag(i)*(smeantdt+smean0)
+ scur(i)=dstate(i)
+ dstate(i+6)=ee(i)
+ dstate(i+12)=f4*sdevtau
+ end do
+c
return
end
c
@@ -240,9 +389,7 @@
& ierr,errstrng)
c
c... subroutine to compute the current stress and updated material
-c matrix for the time-dependent solution. Since this is a purely
-c elastic material, the material matrix should not change unless the
-c material properties have changed for a time step.
+c matrix for the time-dependent solution.
c Note that only the upper triangle is used (or available), as
c dmat is assumed to always be symmetric.
c
@@ -275,10 +422,14 @@
include "rgiter_def.inc"
include "ntimdat_def.inc"
c
-c... return error code, as this material is not yet defined
+cdebug write(6,*) "Hello from td_strs_mat_8_f!"
c
- ierr=101
- errstrng="td_strs_mat_8"
+ call td_matinit_8(state,dstate,state0,dmat,prop,rtimdat,rgiter,
+ & ntimdat,iddmat,tmax,nstate,nstate0,nprop,matchg,ierr,errstrng)
+ if(ierr.ne.izero) return
+ call td_strs_8(state,dstate,state0,ee,scur,dmat,prop,rtimdat,
+ & rgiter,ntimdat,iddmat,tmax,nstate,nstate0,nprop,matchg,ierr,
+ & errstrng)
c
return
end
@@ -312,6 +463,8 @@
c
double precision ptmp(10)
c
+cdebug write(6,*) "Hello from prestr_mat_8_f!"
+c
call dcopy(nprop,prop,ione,ptmp,ione)
if(ipauto.eq.ione) then
ptmp(2)=tyoungs
@@ -320,8 +473,8 @@
call elas_mat_8(dmat,ptmp,iddmat,nprop,ierr,errstrng)
return
end
-c
-c
+c
+c
subroutine get_state_8(state,dstate,sout,nstate)
c
c... routine to transfer state variables into sout
@@ -338,6 +491,10 @@
integer nstate
double precision state(nstate),dstate(nstate),sout(3*nstatesmax)
c
+cdebug write(6,*) "Hello from get_state_8_f!"
+c
+ call dcopy(nstate,state,ione,sout,ione)
+ call dcopy(nstate,dstate,ione,sout(nstatesmax+ione),ione)
return
end
c
@@ -365,16 +522,26 @@
c
c... local data
c
+ double precision sub
+ data sub/-1.0d0/
c
-c... local variables
+cdebug integer idb
c
+cdebug write(6,*) "Hello from update_state_8_f!"
c
+cdebug write(6,*) "state-begin:",(state(idb),idb=1,nstate)
+cdebug write(6,*) "dstate-begin:",(dstate(idb),idb=1,nstate)
+ call daxpy(nstate-6,sub,state,ione,dstate,ione)
+ call daxpy(nstate,one,dstate,ione,state,ione)
+cdebug write(6,*) "state-end:",(state(idb),idb=1,nstate)
+cdebug write(6,*) "dstate-end:",(dstate(idb),idb=1,nstate)
+c
return
end
c
-c
+c
c version
-c $Id: mat_8.f,v 1.5 2005/04/01 23:10:17 willic3 Exp $
+c $Id: mat_8.f,v 1.8 2005/04/01 23:07:43 willic3 Exp $
c Generated automatically by Fortran77Mill on Tue May 18 14:18:50 2004
More information about the cig-commits
mailing list